[libc-commits] [libc] [libc][math] Add skip accurate pass option for exp*, log*, and powf functions. (PR #129831)
via libc-commits
libc-commits at lists.llvm.org
Tue Mar 4 21:27:54 PST 2025
https://github.com/lntue created https://github.com/llvm/llvm-project/pull/129831
None
>From ef73fac52209348cc5e86107de596ffb4297959f Mon Sep 17 00:00:00 2001
From: Tue Ly <lntue.h at gmail.com>
Date: Wed, 5 Mar 2025 05:25:27 +0000
Subject: [PATCH] [libc][math] Add skip accurate pass option for exp*, log*,
and powf functions.
---
libc/src/math/generic/common_constants.cpp | 24 ++--
libc/src/math/generic/exp.cpp | 22 +++
libc/src/math/generic/exp10.cpp | 19 +++
libc/src/math/generic/exp2.cpp | 18 +++
libc/src/math/generic/explogxf.cpp | 8 +-
libc/src/math/generic/explogxf.h | 34 +++--
libc/src/math/generic/log.cpp | 12 +-
libc/src/math/generic/log10.cpp | 12 +-
libc/src/math/generic/log1p.cpp | 40 +++---
libc/src/math/generic/log2.cpp | 22 +--
libc/src/math/generic/powf.cpp | 153 +++++++++++++++++++++
11 files changed, 304 insertions(+), 60 deletions(-)
diff --git a/libc/src/math/generic/common_constants.cpp b/libc/src/math/generic/common_constants.cpp
index e29c083fa474a..3088ef96e3b93 100644
--- a/libc/src/math/generic/common_constants.cpp
+++ b/libc/src/math/generic/common_constants.cpp
@@ -112,7 +112,7 @@ const double LOG_F[128] = {
// precision, and -2^-8 <= v < 2^-7.
// TODO(lntue): Add reference to how the constants are derived after the
// resulting paper is ready.
-alignas(32) const float R[128] = {
+alignas(8) const float R[128] = {
0x1p0, 0x1.fcp-1, 0x1.f8p-1, 0x1.f4p-1, 0x1.fp-1, 0x1.ecp-1, 0x1.e8p-1,
0x1.e4p-1, 0x1.ep-1, 0x1.dep-1, 0x1.dap-1, 0x1.d6p-1, 0x1.d4p-1, 0x1.dp-1,
0x1.ccp-1, 0x1.cap-1, 0x1.c6p-1, 0x1.c4p-1, 0x1.cp-1, 0x1.bep-1, 0x1.bap-1,
@@ -133,7 +133,7 @@ alignas(32) const float R[128] = {
0x1.0ap-1, 0x1.08p-1, 0x1.08p-1, 0x1.06p-1, 0x1.06p-1, 0x1.04p-1, 0x1.04p-1,
0x1.02p-1, 0x1.0p-1};
-alignas(64) const double RD[128] = {
+const double RD[128] = {
0x1p0, 0x1.fcp-1, 0x1.f8p-1, 0x1.f4p-1, 0x1.fp-1, 0x1.ecp-1, 0x1.e8p-1,
0x1.e4p-1, 0x1.ep-1, 0x1.dep-1, 0x1.dap-1, 0x1.d6p-1, 0x1.d4p-1, 0x1.dp-1,
0x1.ccp-1, 0x1.cap-1, 0x1.c6p-1, 0x1.c4p-1, 0x1.cp-1, 0x1.bep-1, 0x1.bap-1,
@@ -158,7 +158,7 @@ alignas(64) const double RD[128] = {
// available.
// Generated by Sollya with the formula: CD[i] = RD[i]*(1 + i*2^-7) - 1
// for RD[i] defined on the table above.
-alignas(64) const double CD[128] = {
+const double CD[128] = {
0.0, -0x1p-14, -0x1p-12, -0x1.2p-11, -0x1p-10, -0x1.9p-10,
-0x1.2p-9, -0x1.88p-9, -0x1p-8, -0x1.9p-11, -0x1.fp-10, -0x1.9cp-9,
-0x1p-12, -0x1.cp-10, -0x1.bp-9, -0x1.5p-11, -0x1.4p-9, 0x1p-14,
@@ -183,7 +183,7 @@ alignas(64) const double CD[128] = {
-0x1p-14, -0x1p-8,
};
-alignas(64) const double LOG_R[128] = {
+const double LOG_R[128] = {
0x0.0000000000000p0, 0x1.010157588de71p-7, 0x1.0205658935847p-6,
0x1.8492528c8cabfp-6, 0x1.0415d89e74444p-5, 0x1.466aed42de3eap-5,
0x1.894aa149fb343p-5, 0x1.ccb73cdddb2ccp-5, 0x1.08598b59e3a07p-4,
@@ -228,7 +228,7 @@ alignas(64) const double LOG_R[128] = {
0x1.5707a26bb8c66p-1, 0x1.5af405c3649ep-1, 0x1.5af405c3649ep-1,
0x1.5ee82aa24192p-1, 0x0.000000000000p0};
-alignas(64) const double LOG2_R[128] = {
+const double LOG2_R[128] = {
0x0.0000000000000p+0, 0x1.72c7ba20f7327p-7, 0x1.743ee861f3556p-6,
0x1.184b8e4c56af8p-5, 0x1.77394c9d958d5p-5, 0x1.d6ebd1f1febfep-5,
0x1.1bb32a600549dp-4, 0x1.4c560fe68af88p-4, 0x1.7d60496cfbb4cp-4,
@@ -281,7 +281,7 @@ alignas(64) const double LOG2_R[128] = {
// print("{", -c, ",", -b, "},");
// };
// We replace LOG_R[0] with log10(1.0) == 0.0
-alignas(64) const NumberPair<double> LOG_R_DD[128] = {
+alignas(16) const NumberPair<double> LOG_R_DD[128] = {
{0.0, 0.0},
{-0x1.0c76b999d2be8p-46, 0x1.010157589p-7},
{-0x1.3dc5b06e2f7d2p-45, 0x1.0205658938p-6},
@@ -417,7 +417,7 @@ alignas(64) const NumberPair<double> LOG_R_DD[128] = {
// Output range:
// [-0x1.3ffcp-15, 0x1.3e3dp-15]
// We store S2[i] = 2^16 (r(i - 2^6) - 1).
-alignas(64) const int S2[193] = {
+alignas(8) const int S2[193] = {
0x101, 0xfd, 0xf9, 0xf5, 0xf1, 0xed, 0xe9, 0xe5, 0xe1,
0xdd, 0xd9, 0xd5, 0xd1, 0xcd, 0xc9, 0xc5, 0xc1, 0xbd,
0xb9, 0xb4, 0xb0, 0xac, 0xa8, 0xa4, 0xa0, 0x9c, 0x98,
@@ -441,7 +441,7 @@ alignas(64) const int S2[193] = {
-0x1cd, -0x1d1, -0x1d5, -0x1d9, -0x1dd, -0x1e0, -0x1e4, -0x1e8, -0x1ec,
-0x1f0, -0x1f4, -0x1f8, -0x1fc};
-alignas(64) const double R2[193] = {
+const double R2[193] = {
0x1.0101p0, 0x1.00fdp0, 0x1.00f9p0, 0x1.00f5p0, 0x1.00f1p0,
0x1.00edp0, 0x1.00e9p0, 0x1.00e5p0, 0x1.00e1p0, 0x1.00ddp0,
0x1.00d9p0, 0x1.00d5p0, 0x1.00d1p0, 0x1.00cdp0, 0x1.00c9p0,
@@ -488,7 +488,7 @@ alignas(64) const double R2[193] = {
// Output range:
// [-0x1.01928p-22 , 0x1p-22]
// We store S[i] = 2^21 (r(i - 80) - 1).
-alignas(64) const int S3[161] = {
+alignas(8) const int S3[161] = {
0x50, 0x4f, 0x4e, 0x4d, 0x4c, 0x4b, 0x4a, 0x49, 0x48, 0x47, 0x46,
0x45, 0x44, 0x43, 0x42, 0x41, 0x40, 0x3f, 0x3e, 0x3d, 0x3c, 0x3b,
0x3a, 0x39, 0x38, 0x37, 0x36, 0x35, 0x34, 0x33, 0x32, 0x31, 0x30,
@@ -511,7 +511,7 @@ alignas(64) const int S3[161] = {
// Output range:
// [-0x1.0002143p-29 , 0x1p-29]
// We store S[i] = 2^28 (r(i - 65) - 1).
-alignas(64) const int S4[130] = {
+alignas(8) const int S4[130] = {
0x41, 0x40, 0x3f, 0x3e, 0x3d, 0x3c, 0x3b, 0x3a, 0x39, 0x38, 0x37,
0x36, 0x35, 0x34, 0x33, 0x32, 0x31, 0x30, 0x2f, 0x2e, 0x2d, 0x2c,
0x2b, 0x2a, 0x29, 0x28, 0x27, 0x26, 0x25, 0x24, 0x23, 0x22, 0x21,
@@ -661,7 +661,7 @@ const double EXP_M2[128] = {
// d = round(a - b - c, D, RN);
// print("{", d, ",", c, ",", b, "},");
// };
-const fputil::TripleDouble EXP2_MID1[64] = {
+alignas(16) const fputil::TripleDouble EXP2_MID1[64] = {
{0, 0, 0x1p0},
{-0x1.9085b0a3d74d5p-110, -0x1.19083535b085dp-56, 0x1.02c9a3e778061p0},
{0x1.05ff94f8d257ep-110, 0x1.d73e2a475b465p-55, 0x1.059b0d3158574p0},
@@ -739,7 +739,7 @@ const fputil::TripleDouble EXP2_MID1[64] = {
// d = round(a - b - c, D, RN);
// print("{", d, ",", c, ",", b, "},");
// };
-const fputil::TripleDouble EXP2_MID2[64] = {
+alignas(16) const fputil::TripleDouble EXP2_MID2[64] = {
{0, 0, 0x1p0},
{0x1.39726694630e3p-108, 0x1.ae8e38c59c72ap-54, 0x1.000b175effdc7p0},
{0x1.e5e06ddd31156p-112, -0x1.7b5d0d58ea8f4p-58, 0x1.00162f3904052p0},
diff --git a/libc/src/math/generic/exp.cpp b/libc/src/math/generic/exp.cpp
index 38b683aa01166..143800ca078a6 100644
--- a/libc/src/math/generic/exp.cpp
+++ b/libc/src/math/generic/exp.cpp
@@ -39,8 +39,11 @@ constexpr double LOG2_E = 0x1.71547652b82fep+0;
// Error bounds:
// Errors when using double precision.
constexpr double ERR_D = 0x1.8p-63;
+
+#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
// Errors when using double-double precision.
constexpr double ERR_DD = 0x1.0p-99;
+#endif // LIBC_MATH_HAS_SKIP_ACCURATE_PASS
// -2^-12 * log(2)
// > a = -2^-12 * log(2);
@@ -50,8 +53,11 @@ constexpr double ERR_DD = 0x1.0p-99;
// Errors < 1.5 * 2^-133
constexpr double MLOG_2_EXP2_M12_HI = -0x1.62e42ffp-13;
constexpr double MLOG_2_EXP2_M12_MID = 0x1.718432a1b0e26p-47;
+
+#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
constexpr double MLOG_2_EXP2_M12_MID_30 = 0x1.718432ap-47;
constexpr double MLOG_2_EXP2_M12_LO = 0x1.b0e2633fe0685p-79;
+#endif // LIBC_MATH_HAS_SKIP_ACCURATE_PASS
namespace {
@@ -72,6 +78,7 @@ LIBC_INLINE double poly_approx_d(double dx) {
return p;
}
+#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
// Polynomial approximation with double-double precision:
// Return exp(dx) ~ 1 + dx + dx^2 / 2 + ... + dx^6 / 720
// For |dx| < 2^-13 + 2^-30:
@@ -171,6 +178,7 @@ DoubleDouble exp_double_double(double x, double kd,
return r;
}
+#endif // LIBC_MATH_HAS_SKIP_ACCURATE_PASS
// Check for exceptional cases when
// |x| <= 2^-53 or x < log(2^-1075) or x >= 0x1.6232bdd7abcd3p+9
@@ -373,6 +381,19 @@ LLVM_LIBC_FUNCTION(double, exp, (double x)) {
double lo = fputil::multiply_add(p, mid_lo, exp_mid.lo);
+#ifdef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
+ if (LIBC_UNLIKELY(denorm)) {
+ return ziv_test_denorm</*SKIP_ZIV_TEST=*/true>(hi, exp_mid.hi, lo, ERR_D)
+ .value();
+ } else {
+ // to multiply by 2^hi, a fast way is to simply add hi to the exponent
+ // field.
+ int64_t exp_hi = static_cast<int64_t>(hi) << FPBits::FRACTION_LEN;
+ double r =
+ cpp::bit_cast<double>(exp_hi + cpp::bit_cast<int64_t>(exp_mid.hi + lo));
+ return r;
+ }
+#else
if (LIBC_UNLIKELY(denorm)) {
if (auto r = ziv_test_denorm(hi, exp_mid.hi, lo, ERR_D);
LIBC_LIKELY(r.has_value()))
@@ -413,6 +434,7 @@ LLVM_LIBC_FUNCTION(double, exp, (double x)) {
Float128 r_f128 = exp_f128(x, kd, idx1, idx2);
return static_cast<double>(r_f128);
+#endif // LIBC_MATH_HAS_SKIP_ACCURATE_PASS
}
} // namespace LIBC_NAMESPACE_DECL
diff --git a/libc/src/math/generic/exp10.cpp b/libc/src/math/generic/exp10.cpp
index 748c8a22b2368..c464979b092c3 100644
--- a/libc/src/math/generic/exp10.cpp
+++ b/libc/src/math/generic/exp10.cpp
@@ -44,15 +44,20 @@ constexpr double LOG2_10 = 0x1.a934f0979a371p+1;
// Errors < 1.5 * 2^-144
constexpr double MLOG10_2_EXP2_M12_HI = -0x1.3441350ap-14;
constexpr double MLOG10_2_EXP2_M12_MID = 0x1.0c0219dc1da99p-51;
+
+#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
constexpr double MLOG10_2_EXP2_M12_MID_32 = 0x1.0c0219dcp-51;
constexpr double MLOG10_2_EXP2_M12_LO = 0x1.da994fd20dba2p-87;
+#endif // LIBC_MATH_HAS_SKIP_ACCURATE_PASS
// Error bounds:
// Errors when using double precision.
constexpr double ERR_D = 0x1.8p-63;
+#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
// Errors when using double-double precision.
constexpr double ERR_DD = 0x1.8p-99;
+#endif // LIBC_MATH_HAS_SKIP_ACCURATE_PASS
namespace {
@@ -72,6 +77,7 @@ LIBC_INLINE double poly_approx_d(double dx) {
return p;
}
+#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
// Polynomial approximation with double-double precision. Generated by Solya
// with:
// > P = fpminimax((10^x - 1)/x, 5, [|DD...|], [-2^-14, 2^-14]);
@@ -171,6 +177,7 @@ DoubleDouble exp10_double_double(double x, double kd,
return r;
}
+#endif // LIBC_MATH_HAS_SKIP_ACCURATE_PASS
// When output is denormal.
double exp10_denorm(double x) {
@@ -199,6 +206,10 @@ double exp10_denorm(double x) {
double lo = fputil::multiply_add(p, mid_lo, exp_mid.lo);
+#ifdef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
+ return ziv_test_denorm</*SKIP_ZIV_TEST=*/true>(hi, exp_mid.hi, lo, ERR_D)
+ .value();
+#else
if (auto r = ziv_test_denorm(hi, exp_mid.hi, lo, ERR_D);
LIBC_LIKELY(r.has_value()))
return r.value();
@@ -214,6 +225,7 @@ double exp10_denorm(double x) {
Float128 r_f128 = exp10_f128(x, kd, idx1, idx2);
return static_cast<double>(r_f128);
+#endif // LIBC_MATH_HAS_SKIP_ACCURATE_PASS
}
// Check for exceptional cases when:
@@ -391,6 +403,12 @@ LLVM_LIBC_FUNCTION(double, exp10, (double x)) {
double lo = fputil::multiply_add(p, mid_lo, exp_mid.lo);
+#ifdef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
+ int64_t exp_hi = static_cast<int64_t>(hi) << FPBits::FRACTION_LEN;
+ double r =
+ cpp::bit_cast<double>(exp_hi + cpp::bit_cast<int64_t>(exp_mid.hi + lo));
+ return r;
+#else
double upper = exp_mid.hi + (lo + ERR_D);
double lower = exp_mid.hi + (lo - ERR_D);
@@ -473,6 +491,7 @@ LLVM_LIBC_FUNCTION(double, exp10, (double x)) {
Float128 r_f128 = exp10_f128(x, kd, idx1, idx2);
return static_cast<double>(r_f128);
+#endif // LIBC_MATH_HAS_SKIP_ACCURATE_PASS
}
} // namespace LIBC_NAMESPACE_DECL
diff --git a/libc/src/math/generic/exp2.cpp b/libc/src/math/generic/exp2.cpp
index 935548b538b94..2c612777c9cb5 100644
--- a/libc/src/math/generic/exp2.cpp
+++ b/libc/src/math/generic/exp2.cpp
@@ -41,8 +41,10 @@ constexpr double ERR_D = 0x1.0p-63;
constexpr double ERR_D = 0x1.8p-63;
#endif // LIBC_TARGET_CPU_HAS_FMA
+#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
// Errors when using double-double precision.
constexpr double ERR_DD = 0x1.0p-100;
+#endif // LIBC_MATH_HAS_SKIP_ACCURATE_PASS
namespace {
@@ -62,6 +64,7 @@ LIBC_INLINE double poly_approx_d(double dx) {
return p;
}
+#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
// Polynomial approximation with double-double precision. Generated by Solya
// with:
// > P = fpminimax((2^x - 1)/x, 5, [|DD...|], [-2^-13 - 2^-30, 2^-13 + 2^-30]);
@@ -147,6 +150,7 @@ DoubleDouble exp2_double_double(double x, const DoubleDouble &exp_mid) {
return r;
}
+#endif // LIBC_MATH_HAS_SKIP_ACCURATE_PASS
// When output is denormal.
double exp2_denorm(double x) {
@@ -174,6 +178,10 @@ double exp2_denorm(double x) {
double lo = fputil::multiply_add(p, mid_lo, exp_mid.lo);
+#ifdef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
+ return ziv_test_denorm</*SKIP_ZIV_TEST=*/true>(hi, exp_mid.hi, lo, ERR_D)
+ .value();
+#else
if (auto r = ziv_test_denorm(hi, exp_mid.hi, lo, ERR_D);
LIBC_LIKELY(r.has_value()))
return r.value();
@@ -189,6 +197,7 @@ double exp2_denorm(double x) {
Float128 r_f128 = exp2_f128(dx, hi, idx1, idx2);
return static_cast<double>(r_f128);
+#endif // LIBC_MATH_HAS_SKIP_ACCURATE_PASS
}
// Check for exceptional cases when:
@@ -358,6 +367,14 @@ LLVM_LIBC_FUNCTION(double, exp2, (double x)) {
double lo = fputil::multiply_add(p, mid_lo, exp_mid.lo);
+#ifdef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
+ // To multiply by 2^hi, a fast way is to simply add hi to the exponent
+ // field.
+ int64_t exp_hi = static_cast<int64_t>(hi) << FPBits::FRACTION_LEN;
+ double r =
+ cpp::bit_cast<double>(exp_hi + cpp::bit_cast<int64_t>(exp_mid.hi + lo));
+ return r;
+#else
double upper = exp_mid.hi + (lo + ERR_D);
double lower = exp_mid.hi + (lo - ERR_D);
@@ -387,6 +404,7 @@ LLVM_LIBC_FUNCTION(double, exp2, (double x)) {
Float128 r_f128 = exp2_f128(dx, hi, idx1, idx2);
return static_cast<double>(r_f128);
+#endif // LIBC_MATH_HAS_SKIP_ACCURATE_PASS
}
} // namespace LIBC_NAMESPACE_DECL
diff --git a/libc/src/math/generic/explogxf.cpp b/libc/src/math/generic/explogxf.cpp
index 9e945eca7aed6..d38efa0269693 100644
--- a/libc/src/math/generic/explogxf.cpp
+++ b/libc/src/math/generic/explogxf.cpp
@@ -12,7 +12,7 @@
namespace LIBC_NAMESPACE_DECL {
// N[Table[Log[2, 1 + x], {x, 0/64, 63/64, 1/64}], 40]
-alignas(64) const double LOG_P1_LOG2[LOG_P1_SIZE] = {
+alignas(8) const double LOG_P1_LOG2[LOG_P1_SIZE] = {
0x0.0000000000000p+0, 0x1.6e79685c2d22ap-6, 0x1.6bad3758efd87p-5,
0x1.0eb389fa29f9bp-4, 0x1.663f6fac91316p-4, 0x1.bc84240adabbap-4,
0x1.08c588cda79e4p-3, 0x1.32ae9e278ae1ap-3, 0x1.5c01a39fbd688p-3,
@@ -38,7 +38,7 @@ alignas(64) const double LOG_P1_LOG2[LOG_P1_SIZE] = {
};
// N[Table[1/(1 + x), {x, 0/64, 63/64, 1/64}], 40]
-alignas(64) const double LOG_P1_1_OVER[LOG_P1_SIZE] = {
+alignas(8) const double LOG_P1_1_OVER[LOG_P1_SIZE] = {
0x1.0000000000000p+0, 0x1.f81f81f81f820p-1, 0x1.f07c1f07c1f08p-1,
0x1.e9131abf0b767p-1, 0x1.e1e1e1e1e1e1ep-1, 0x1.dae6076b981dbp-1,
0x1.d41d41d41d41dp-1, 0x1.cd85689039b0bp-1, 0x1.c71c71c71c71cp-1,
@@ -64,11 +64,11 @@ alignas(64) const double LOG_P1_1_OVER[LOG_P1_SIZE] = {
// Taylos series expansion for Log[2, 1 + x] splitted to EVEN AND ODD numbers
// K_LOG2_ODD starts from x^3
-alignas(64) const
+alignas(8) const
double K_LOG2_ODD[4] = {0x1.ec709dc3a03fdp-2, 0x1.2776c50ef9bfep-2,
0x1.a61762a7aded9p-3, 0x1.484b13d7c02a9p-3};
-alignas(64) const
+alignas(8) const
double K_LOG2_EVEN[4] = {-0x1.71547652b82fep-1, -0x1.71547652b82fep-2,
-0x1.ec709dc3a03fdp-3, -0x1.2776c50ef9bfep-3};
diff --git a/libc/src/math/generic/explogxf.h b/libc/src/math/generic/explogxf.h
index 651524a165f03..e79aa13eb57f7 100644
--- a/libc/src/math/generic/explogxf.h
+++ b/libc/src/math/generic/explogxf.h
@@ -338,8 +338,9 @@ LIBC_INLINE static double log_eval(double x) {
// double(2^-1022 + x) - 2^-1022 = double(x).
// So if we scale x up by 2^1022, we can use
// double(1.0 + 2^1022 * x) - 1.0 to test how x is rounded in denormal range.
-LIBC_INLINE cpp::optional<double> ziv_test_denorm(int hi, double mid, double lo,
- double err) {
+template <bool SKIP_ZIV_TEST = false>
+LIBC_INLINE static cpp::optional<double>
+ziv_test_denorm(int hi, double mid, double lo, double err) {
using FPBits = typename fputil::FPBits<double>;
// Scaling factor = 1/(min normal number) = 2^1022
@@ -360,22 +361,27 @@ LIBC_INLINE cpp::optional<double> ziv_test_denorm(int hi, double mid, double lo,
scale_down = 0x3FF0'0000'0000'0000; // 1023 in the exponent field.
}
- double err_scaled =
- cpp::bit_cast<double>(exp_hi + cpp::bit_cast<int64_t>(err));
-
- double lo_u = lo_scaled + err_scaled;
- double lo_l = lo_scaled - err_scaled;
-
// By adding 1.0, the results will have similar rounding points as denormal
// outputs.
- double upper = extra_factor + (mid_hi + lo_u);
- double lower = extra_factor + (mid_hi + lo_l);
+ if constexpr (SKIP_ZIV_TEST) {
+ double r = extra_factor + (mid_hi + lo_scaled);
+ return cpp::bit_cast<double>(cpp::bit_cast<uint64_t>(r) - scale_down);
+ } else {
+ double err_scaled =
+ cpp::bit_cast<double>(exp_hi + cpp::bit_cast<int64_t>(err));
- if (LIBC_LIKELY(upper == lower)) {
- return cpp::bit_cast<double>(cpp::bit_cast<uint64_t>(upper) - scale_down);
- }
+ double lo_u = lo_scaled + err_scaled;
+ double lo_l = lo_scaled - err_scaled;
- return cpp::nullopt;
+ double upper = extra_factor + (mid_hi + lo_u);
+ double lower = extra_factor + (mid_hi + lo_l);
+
+ if (LIBC_LIKELY(upper == lower)) {
+ return cpp::bit_cast<double>(cpp::bit_cast<uint64_t>(upper) - scale_down);
+ }
+
+ return cpp::nullopt;
+ }
}
} // namespace LIBC_NAMESPACE_DECL
diff --git a/libc/src/math/generic/log.cpp b/libc/src/math/generic/log.cpp
index 4302c64c8abac..04eebab975cd5 100644
--- a/libc/src/math/generic/log.cpp
+++ b/libc/src/math/generic/log.cpp
@@ -30,6 +30,7 @@ using LIBC_NAMESPACE::operator""_u128;
namespace {
+#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
// A simple upper bound for the error of e_x * log(2) - log(r).
constexpr double HI_ERR = 0x1.0p-85;
@@ -47,7 +48,7 @@ constexpr double P_ERR = 0x1.0p-50;
constexpr Float128 LOG_2(Sign::POS, /*exponent=*/-128, /*mantissa=*/
0xb17217f7'd1cf79ab'c9e3b398'03f2f6af_u128);
-alignas(64) constexpr LogRR LOG_TABLE = {
+alignas(16) constexpr LogRR LOG_TABLE = {
// -log(r) with 128-bit precision generated by SageMath with:
// for i in range(128):
// r = 2^-8 * ceil( 2^8 * (1 - 2^(-8)) / (1 + i*2^(-7)) );
@@ -731,6 +732,7 @@ double log_accurate(int e_x, int index, double m_x) {
return static_cast<double>(r);
}
+#endif // LIBC_MATH_HAS_SKIP_ACCURATE_PASS
} // namespace
@@ -794,7 +796,7 @@ LLVM_LIBC_FUNCTION(double, log, (double x)) {
uint64_t x_m = (x_u & 0x000F'FFFF'FFFF'FFFFULL) | 0x3FF0'0000'0000'0000ULL;
double m = FPBits_t(x_m).get_val();
- double u, u_sq, err;
+ double u, u_sq;
fputil::DoubleDouble r1;
// Perform exact range reduction
@@ -819,12 +821,15 @@ LLVM_LIBC_FUNCTION(double, log, (double x)) {
double p2 = fputil::multiply_add(u, LOG_COEFFS[5], LOG_COEFFS[4]);
double p = fputil::polyeval(u_sq, lo + r1.lo, p0, p1, p2);
+#ifdef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
+ return r1.hi + p;
+#else
// Technicallly error of r1.lo is bounded by:
// hi*ulp(log(2)_lo) + C*ulp(u^2)
// To simplify the error computation a bit, we replace |hi|*ulp(log(2)_lo)
// with the upper bound: 2^11 * ulp(log(2)_lo) = 2^-85.
// Total error is bounded by ~ C * ulp(u^2) + 2^-85.
- err = fputil::multiply_add(u_sq, P_ERR, HI_ERR);
+ double err = fputil::multiply_add(u_sq, P_ERR, HI_ERR);
// Lower bound from the result
double left = r1.hi + (p - err);
@@ -836,6 +841,7 @@ LLVM_LIBC_FUNCTION(double, log, (double x)) {
return left;
return log_accurate(x_e, index, u);
+#endif // LIBC_MATH_HAS_SKIP_ACCURATE_PASS
}
} // namespace LIBC_NAMESPACE_DECL
diff --git a/libc/src/math/generic/log10.cpp b/libc/src/math/generic/log10.cpp
index 7df57ef85b81b..fd8d5a8aae938 100644
--- a/libc/src/math/generic/log10.cpp
+++ b/libc/src/math/generic/log10.cpp
@@ -33,6 +33,7 @@ namespace {
constexpr fputil::DoubleDouble LOG10_E = {0x1.95355baaafad3p-57,
0x1.bcb7b1526e50ep-2};
+#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
// A simple upper bound for the error of e_x * log(2) - log(r).
constexpr double HI_ERR = 0x1.0p-85;
@@ -50,7 +51,7 @@ constexpr double P_ERR = 0x1.0p-51;
constexpr Float128 LOG10_2(Sign::POS, /*exponent=*/-129, /*mantissa=*/
0x9a209a84'fbcff798'8f8959ac'0b7c9178_u128);
-alignas(64) constexpr LogRR LOG10_TABLE = {
+alignas(16) constexpr LogRR LOG10_TABLE = {
// -log10(r) with 128-bit precision generated by SageMath with:
//
// for i in range(128):
@@ -733,6 +734,7 @@ double log10_accurate(int e_x, int index, double m_x) {
return static_cast<double>(r);
}
+#endif // LIBC_MATH_HAS_SKIP_ACCURATE_PASS
} // namespace
@@ -796,7 +798,7 @@ LLVM_LIBC_FUNCTION(double, log10, (double x)) {
uint64_t x_m = (x_u & 0x000F'FFFF'FFFF'FFFFULL) | 0x3FF0'0000'0000'0000ULL;
double m = FPBits_t(x_m).get_val();
- double u, u_sq, err;
+ double u, u_sq;
fputil::DoubleDouble r1;
// Perform exact range reduction
@@ -827,12 +829,15 @@ LLVM_LIBC_FUNCTION(double, log10, (double x)) {
// 4*ulp( ulp(r2.hi) )
fputil::DoubleDouble r2 = fputil::quick_mult(r1, LOG10_E);
+#ifdef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
+ return r2.hi + r2.lo;
+#else
// Technicallly error of r1.lo is bounded by:
// |hi|*ulp(log(2)_lo) + C*ulp(u^2)
// To simplify the error computation a bit, we replace |hi|*ulp(log(2)_lo)
// with the upper bound: 2^11 * ulp(log(2)_lo) = 2^-85.
// Total error is bounded by ~ C * ulp(u^2) + 2^-85.
- err = fputil::multiply_add(u_sq, P_ERR, HI_ERR);
+ double err = fputil::multiply_add(u_sq, P_ERR, HI_ERR);
// Lower bound from the result
double left = r2.hi + (r2.lo - err);
@@ -897,6 +902,7 @@ LLVM_LIBC_FUNCTION(double, log10, (double x)) {
}
return log10_accurate(x_e, index, u);
+#endif // LIBC_MATH_HAS_SKIP_ACCURATE_PASS
}
} // namespace LIBC_NAMESPACE_DECL
diff --git a/libc/src/math/generic/log1p.cpp b/libc/src/math/generic/log1p.cpp
index b9c58b843a240..b1f02164b6a28 100644
--- a/libc/src/math/generic/log1p.cpp
+++ b/libc/src/math/generic/log1p.cpp
@@ -29,21 +29,6 @@ using LIBC_NAMESPACE::operator""_u128;
namespace {
-// Extra errors from P is from using x^2 to reduce evaluation latency and
-// directional rounding.
-constexpr double P_ERR = 0x1.0p-49;
-
-// log(2) with 128-bit precision generated by SageMath with:
-// def format_hex(value):
-// l = hex(value)[2:]
-// n = 8
-// x = [l[i:i + n] for i in range(0, len(l), n)]
-// return "0x" + "'".join(x) + "_u128"
-// (s, m, e) = RealField(128)(2).log().sign_mantissa_exponent();
-// print(format_hex(m));
-constexpr Float128 LOG_2(Sign::POS, /*exponent=*/-128, /*mantissa=*/
- 0xb17217f7'd1cf79ab'c9e3b398'03f2f6af_u128);
-
// R1[i] = 2^-8 * nearestint( 2^8 / (1 + i * 2^-7) )
constexpr double R1[129] = {
0x1p0, 0x1.fcp-1, 0x1.f8p-1, 0x1.f4p-1, 0x1.fp-1, 0x1.ecp-1, 0x1.eap-1,
@@ -105,7 +90,7 @@ constexpr double R1[129] = {
// print("{", -c, ",", -b, "},");
// };
// We replace LOG_R1_DD[128] with log(1.0) == 0.0
-constexpr fputil::DoubleDouble LOG_R1_DD[129] = {
+alignas(16) constexpr fputil::DoubleDouble LOG_R1_DD[129] = {
{0.0, 0.0},
{-0x1.0c76b999d2be8p-46, 0x1.010157589p-7},
{-0x1.3dc5b06e2f7d2p-45, 0x1.0205658938p-6},
@@ -248,6 +233,22 @@ constexpr double P_COEFFS[6] = {-0x1p-1,
-0x1.555874ce8ce22p-3,
0x1.24335555ddbe5p-3};
+#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
+// Extra errors from P is from using x^2 to reduce evaluation latency and
+// directional rounding.
+constexpr double P_ERR = 0x1.0p-49;
+
+// log(2) with 128-bit precision generated by SageMath with:
+// def format_hex(value):
+// l = hex(value)[2:]
+// n = 8
+// x = [l[i:i + n] for i in range(0, len(l), n)]
+// return "0x" + "'".join(x) + "_u128"
+// (s, m, e) = RealField(128)(2).log().sign_mantissa_exponent();
+// print(format_hex(m));
+constexpr Float128 LOG_2(Sign::POS, /*exponent=*/-128, /*mantissa=*/
+ 0xb17217f7'd1cf79ab'c9e3b398'03f2f6af_u128);
+
// -log(r1) with 128-bit precision generated by SageMath with:
//
// for i in range(129):
@@ -874,6 +875,7 @@ constexpr Float128 BIG_COEFFS[4]{
return static_cast<double>(r);
}
+#endif // LIBC_MATH_HAS_SKIP_ACCURATE_PASS
} // namespace
@@ -969,9 +971,11 @@ LLVM_LIBC_FUNCTION(double, log1p, (double x)) {
// <= 2^11 * 2^(-43-53) = 2^-85
double lo = fputil::multiply_add(e_x, LOG_2_LO, LOG_R1_DD[idx].lo);
+#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
// Error bound of e_x * log(2) - log(r1)
constexpr double ERR_HI[2] = {0x1.0p-85, 0.0};
double err_hi = ERR_HI[hi == 0.0];
+#endif // LIBC_MATH_HAS_SKIP_ACCURATE_PASS
// Scale x_dd by 2^(-xh_bits.get_exponent()).
int64_t s_u = static_cast<int64_t>(x_u & FPBits_t::EXP_MASK) -
@@ -1033,6 +1037,9 @@ LLVM_LIBC_FUNCTION(double, log1p, (double x)) {
double p2 = fputil::multiply_add(v_dd.hi, P_COEFFS[5], P_COEFFS[4]);
double p = fputil::polyeval(v_sq, (v_dd.lo + r1.lo) + lo, p0, p1, p2);
+#ifdef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
+ return r1.hi + p;
+#else
double err = fputil::multiply_add(v_sq, P_ERR, err_hi);
double left = r1.hi + (p - err);
@@ -1043,6 +1050,7 @@ LLVM_LIBC_FUNCTION(double, log1p, (double x)) {
return left;
return log1p_accurate(x_e, idx, v_dd);
+#endif // LIBC_MATH_HAS_SKIP_ACCURATE_PASS
}
} // namespace LIBC_NAMESPACE_DECL
diff --git a/libc/src/math/generic/log2.cpp b/libc/src/math/generic/log2.cpp
index 37ea0c8f13431..f46ff724a4f37 100644
--- a/libc/src/math/generic/log2.cpp
+++ b/libc/src/math/generic/log2.cpp
@@ -33,10 +33,7 @@ namespace {
constexpr fputil::DoubleDouble LOG2_E = {0x1.777d0ffda0d24p-56,
0x1.71547652b82fep0};
-// Extra errors from P is from using x^2 to reduce evaluation latency.
-constexpr double P_ERR = 0x1.0p-49;
-
-const fputil::DoubleDouble LOG_R1[128] = {
+alignas(16) const fputil::DoubleDouble LOG_R1[128] = {
{0.0, 0.0},
{0x1.46662d417cedp-62, 0x1.010157588de71p-7},
{0x1.27c8e8416e71fp-60, 0x1.0205658935847p-6},
@@ -167,7 +164,11 @@ const fputil::DoubleDouble LOG_R1[128] = {
{0.0, 0.0},
};
-alignas(64) constexpr LogRR LOG2_TABLE = {
+#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
+// Extra errors from P is from using x^2 to reduce evaluation latency.
+constexpr double P_ERR = 0x1.0p-49;
+
+alignas(16) constexpr LogRR LOG2_TABLE = {
// -log2(r) with 128-bit precision generated by SageMath with:
// def format_hex(value):
// l = hex(value)[2:]
@@ -853,6 +854,7 @@ double log2_accurate(int e_x, int index, double m_x) {
return static_cast<double>(r);
}
+#endif // LIBC_MATH_HAS_SKIP_ACCURATE_PASS
} // namespace
@@ -909,7 +911,7 @@ LLVM_LIBC_FUNCTION(double, log2, (double x)) {
uint64_t x_m = (x_u & 0x000F'FFFF'FFFF'FFFFULL) | 0x3FF0'0000'0000'0000ULL;
double m = FPBits_t(x_m).get_val();
- double u, u_sq, err;
+ double u, u_sq;
fputil::DoubleDouble r1;
// Perform exact range reduction
@@ -927,8 +929,6 @@ LLVM_LIBC_FUNCTION(double, log2, (double x)) {
// Error of u_sq = ulp(u^2);
u_sq = u * u;
- // Total error is bounded by ~ C * ulp(u^2).
- err = u_sq * P_ERR;
// Degree-7 minimax polynomial
double p0 = fputil::multiply_add(u, LOG_COEFFS[1], LOG_COEFFS[0]);
double p1 = fputil::multiply_add(u, LOG_COEFFS[3], LOG_COEFFS[2]);
@@ -948,6 +948,11 @@ LLVM_LIBC_FUNCTION(double, log2, (double x)) {
// Overall, if we choose sufficiently large constant C, the total error is
// bounded by (C * ulp(u^2)).
+#ifdef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
+ return r3.hi + r3.lo;
+#else
+ // Total error is bounded by ~ C * ulp(u^2).
+ double err = u_sq * P_ERR;
// Lower bound from the result
double left = r3.hi + (r3.lo - err);
// Upper bound from the result
@@ -958,6 +963,7 @@ LLVM_LIBC_FUNCTION(double, log2, (double x)) {
return left;
return log2_accurate(x_e, index, u);
+#endif // LIBC_MATH_HAS_SKIP_ACCURATE_PASS
}
} // namespace LIBC_NAMESPACE_DECL
diff --git a/libc/src/math/generic/powf.cpp b/libc/src/math/generic/powf.cpp
index c84ce0da34b10..7f4417d275702 100644
--- a/libc/src/math/generic/powf.cpp
+++ b/libc/src/math/generic/powf.cpp
@@ -32,6 +32,139 @@ using fputil::TripleDouble;
namespace {
+#ifdef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
+alignas(16) constexpr DoubleDouble LOG2_R_DD[128] = {
+ {0.0, 0.0},
+ {-0x1.177c23362928cp-25, 0x1.72c8p-7},
+ {-0x1.179e0caa9c9abp-22, 0x1.744p-6},
+ {-0x1.c6cea541f5b7p-23, 0x1.184cp-5},
+ {-0x1.66c4d4e554434p-22, 0x1.773ap-5},
+ {-0x1.70700a00fdd55p-24, 0x1.d6ecp-5},
+ {0x1.53002a4e86631p-23, 0x1.1bb3p-4},
+ {0x1.fcd15f101c142p-25, 0x1.4c56p-4},
+ {0x1.25b3eed319cedp-22, 0x1.7d6p-4},
+ {-0x1.4195120d8486fp-22, 0x1.960dp-4},
+ {0x1.45b878e27d0d9p-23, 0x1.c7b5p-4},
+ {0x1.770744593a4cbp-22, 0x1.f9c9p-4},
+ {0x1.c673032495d24p-22, 0x1.097ep-3},
+ {-0x1.1eaa65b49696ep-22, 0x1.22dbp-3},
+ {0x1.b2866f2850b22p-22, 0x1.3c6f8p-3},
+ {0x1.8ee37cd2ea9d3p-25, 0x1.494f8p-3},
+ {0x1.7e86f9c2154fbp-24, 0x1.633a8p-3},
+ {0x1.8e3cfc25f0ce6p-26, 0x1.7046p-3},
+ {0x1.57f7a64ccd537p-28, 0x1.8a898p-3},
+ {-0x1.a761c09fbd2aep-22, 0x1.97c2p-3},
+ {0x1.24bea9a2c66f3p-22, 0x1.b26p-3},
+ {-0x1.60002ccfe43f5p-25, 0x1.bfc68p-3},
+ {0x1.69f220e97f22cp-22, 0x1.dac2p-3},
+ {-0x1.6164f64c210ep-22, 0x1.e858p-3},
+ {-0x1.0c1678ae89767p-24, 0x1.01d9cp-2},
+ {-0x1.f26a05c813d57p-22, 0x1.08bdp-2},
+ {0x1.4d8fc561c8d44p-24, 0x1.169cp-2},
+ {-0x1.362ad8f7ca2dp-22, 0x1.1d984p-2},
+ {0x1.2b13cd6c4d042p-22, 0x1.249ccp-2},
+ {-0x1.1c8f11979a5dbp-22, 0x1.32cp-2},
+ {0x1.c2ab3edefe569p-23, 0x1.39de8p-2},
+ {0x1.7c3eca28e69cap-26, 0x1.4106p-2},
+ {-0x1.34c4e99e1c6c6p-24, 0x1.4f6fcp-2},
+ {-0x1.194a871b63619p-22, 0x1.56b24p-2},
+ {0x1.e3dd5c1c885aep-23, 0x1.5dfdcp-2},
+ {-0x1.6ccf3b1129b7cp-23, 0x1.6552cp-2},
+ {-0x1.2f346e2bf924bp-23, 0x1.6cb1p-2},
+ {-0x1.fa61aaa59c1d8p-23, 0x1.7b8ap-2},
+ {0x1.90c11fd32a3abp-22, 0x1.8304cp-2},
+ {0x1.57f7a64ccd537p-27, 0x1.8a898p-2},
+ {0x1.249ba76fee235p-27, 0x1.9218p-2},
+ {-0x1.aad2729b21ae5p-23, 0x1.99b08p-2},
+ {0x1.71810a5e1818p-22, 0x1.a8ff8p-2},
+ {-0x1.6172fe015e13cp-27, 0x1.b0b68p-2},
+ {0x1.5ec6c1bfbf89ap-24, 0x1.b877cp-2},
+ {0x1.678bf6cdedf51p-24, 0x1.c0438p-2},
+ {0x1.c2d45fe43895ep-22, 0x1.c819cp-2},
+ {-0x1.9ee52ed49d71dp-22, 0x1.cffbp-2},
+ {0x1.5786af187a96bp-27, 0x1.d7e6cp-2},
+ {0x1.3ab0dc56138c9p-23, 0x1.dfdd8p-2},
+ {0x1.fe538ab34efb5p-22, 0x1.e7df4p-2},
+ {-0x1.e4fee07aa4b68p-22, 0x1.efec8p-2},
+ {-0x1.172f32fe67287p-22, 0x1.f804cp-2},
+ {-0x1.9a83ff9ab9cc8p-22, 0x1.00144p-1},
+ {-0x1.68cb06cece193p-22, 0x1.042bep-1},
+ {0x1.8cd71ddf82e2p-22, 0x1.08494p-1},
+ {0x1.5e18ab2df3ae6p-22, 0x1.0c6cap-1},
+ {0x1.5dee4d9d8a273p-25, 0x1.1096p-1},
+ {0x1.fcd15f101c142p-26, 0x1.14c56p-1},
+ {-0x1.2474b0f992ba1p-23, 0x1.18faep-1},
+ {0x1.4b5a92a606047p-24, 0x1.1d368p-1},
+ {0x1.16186fcf54bbdp-22, 0x1.21786p-1},
+ {0x1.18efabeb7d722p-27, 0x1.25c0ap-1},
+ {-0x1.e5fc7d238691dp-24, 0x1.2a0f4p-1},
+ {0x1.f5809faf6283cp-22, 0x1.2e644p-1},
+ {0x1.f5809faf6283cp-22, 0x1.2e644p-1},
+ {0x1.c6e1dcd0cb449p-22, 0x1.32bfep-1},
+ {0x1.76e0e8f74b4d5p-22, 0x1.37222p-1},
+ {-0x1.cb82c89692d99p-24, 0x1.3b8b2p-1},
+ {-0x1.63161c5432aebp-22, 0x1.3ffaep-1},
+ {0x1.458104c41b901p-22, 0x1.44716p-1},
+ {0x1.458104c41b901p-22, 0x1.44716p-1},
+ {-0x1.cd9d0cde578d5p-22, 0x1.48efp-1},
+ {0x1.b9884591add87p-26, 0x1.4d738p-1},
+ {0x1.c6042978605ffp-22, 0x1.51ff2p-1},
+ {-0x1.fc4c96b37dcf6p-22, 0x1.56922p-1},
+ {-0x1.2f346e2bf924bp-24, 0x1.5b2c4p-1},
+ {-0x1.2f346e2bf924bp-24, 0x1.5b2c4p-1},
+ {0x1.c4e4fbb68a4d1p-22, 0x1.5fcdcp-1},
+ {-0x1.9d499bd9b3226p-23, 0x1.6476ep-1},
+ {-0x1.f89b355ede26fp-23, 0x1.69278p-1},
+ {-0x1.f89b355ede26fp-23, 0x1.69278p-1},
+ {0x1.53c7e319f6e92p-24, 0x1.6ddfcp-1},
+ {-0x1.b291f070528c7p-22, 0x1.729fep-1},
+ {0x1.2967a451a7b48p-25, 0x1.7767cp-1},
+ {0x1.2967a451a7b48p-25, 0x1.7767cp-1},
+ {0x1.244fcff690fcep-22, 0x1.7c37ap-1},
+ {0x1.46fd97f5dc572p-23, 0x1.810fap-1},
+ {0x1.46fd97f5dc572p-23, 0x1.810fap-1},
+ {-0x1.f3a7352663e5p-22, 0x1.85efep-1},
+ {0x1.b3cda690370b5p-23, 0x1.8ad84p-1},
+ {0x1.b3cda690370b5p-23, 0x1.8ad84p-1},
+ {0x1.3226b211bf1d9p-23, 0x1.8fc92p-1},
+ {0x1.d24b136c101eep-23, 0x1.94c28p-1},
+ {0x1.d24b136c101eep-23, 0x1.94c28p-1},
+ {0x1.7c40c7907e82ap-22, 0x1.99c48p-1},
+ {-0x1.e81781d97ee91p-22, 0x1.9ecf6p-1},
+ {-0x1.e81781d97ee91p-22, 0x1.9ecf6p-1},
+ {-0x1.6a77813f94e01p-22, 0x1.a3e3p-1},
+ {-0x1.1cfdeb43cfdp-22, 0x1.a8ffap-1},
+ {-0x1.1cfdeb43cfdp-22, 0x1.a8ffap-1},
+ {-0x1.f983f74d3138fp-23, 0x1.ae256p-1},
+ {-0x1.e278ae1a1f51fp-23, 0x1.b3546p-1},
+ {-0x1.e278ae1a1f51fp-23, 0x1.b3546p-1},
+ {-0x1.97552b7b5ea45p-23, 0x1.b88ccp-1},
+ {-0x1.97552b7b5ea45p-23, 0x1.b88ccp-1},
+ {-0x1.19b4f3c72c4f8p-24, 0x1.bdceap-1},
+ {0x1.f7402d26f1a12p-23, 0x1.c31a2p-1},
+ {0x1.f7402d26f1a12p-23, 0x1.c31a2p-1},
+ {-0x1.2056d5dd31d96p-23, 0x1.c86f8p-1},
+ {-0x1.2056d5dd31d96p-23, 0x1.c86f8p-1},
+ {-0x1.6e46335aae723p-24, 0x1.cdcecp-1},
+ {-0x1.beb244c59f331p-22, 0x1.d3382p-1},
+ {-0x1.beb244c59f331p-22, 0x1.d3382p-1},
+ {0x1.16c071e93fd97p-27, 0x1.d8abap-1},
+ {0x1.16c071e93fd97p-27, 0x1.d8abap-1},
+ {0x1.d8175819530c2p-22, 0x1.de298p-1},
+ {0x1.d8175819530c2p-22, 0x1.de298p-1},
+ {0x1.51bd552842c1cp-23, 0x1.e3b2p-1},
+ {0x1.51bd552842c1cp-23, 0x1.e3b2p-1},
+ {0x1.914e204f19d94p-22, 0x1.e9452p-1},
+ {0x1.914e204f19d94p-22, 0x1.e9452p-1},
+ {0x1.c55d997da24fdp-22, 0x1.eee32p-1},
+ {0x1.c55d997da24fdp-22, 0x1.eee32p-1},
+ {-0x1.685c2d2298a6ep-22, 0x1.f48c4p-1},
+ {-0x1.685c2d2298a6ep-22, 0x1.f48c4p-1},
+ {0x1.7a4887bd74039p-22, 0x1.fa406p-1},
+ {0.0, 1.0},
+};
+#else
+
#ifdef LIBC_TARGET_CPU_HAS_FMA
constexpr uint64_t ERR = 64;
#else
@@ -384,6 +517,7 @@ static constexpr DoubleDouble LOG2_R2_DD[] = {
{0x1.3d979ddf4746cp-61, 0x1.6cf6ddd2611d4p-7},
{-0x1.dc930484501f8p-63, 0x1.6fdf461d2e4f8p-7},
};
+#endif // LIBC_MATH_HAS_SKIP_ACCURATE_PASS
LIBC_INLINE bool is_odd_integer(float x) {
using FPBits = typename fputil::FPBits<float>;
@@ -407,6 +541,7 @@ LIBC_INLINE bool is_integer(float x) {
return (x_e + lsb >= UNIT_EXPONENT);
}
+#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
LIBC_INLINE bool larger_exponent(double a, double b) {
using DoubleBits = typename fputil::FPBits<double>;
return DoubleBits(a).get_biased_exponent() >=
@@ -506,6 +641,7 @@ double powf_double_double(int idx_x, double dx, double y6, double lo6_hi,
return cpp::bit_cast<double>(r_bits);
}
+#endif // LIBC_MATH_HAS_SKIP_ACCURATE_PASS
} // namespace
@@ -627,12 +763,14 @@ LLVM_LIBC_FUNCTION(float, powf, (float x, float y)) {
case 0x3f80'0000: // x = 1.0f
return 1.0f;
// TODO: Put these 2 entrypoint dependency under control flag.
+#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
case 0x4000'0000: // x = 2.0f
// pow(2, y) = exp2(y)
return generic::exp2f(y);
case 0x4120'0000: // x = 10.0f
// pow(10, y) = exp10(y)
return generic::exp10f(y);
+#endif // LIBC_MATH_HAS_SKIP_ACCURATE_PASS
}
const bool x_is_neg = x_u >= FloatBits::SIGN_MASK;
@@ -782,6 +920,16 @@ LLVM_LIBC_FUNCTION(float, powf, (float x, float y)) {
// lo6 = 2^6 * lo = 2^6 * (y - (hi + mid)) = y6 * log2(x) - hm.
double y6 = static_cast<double>(y * 0x1.0p6f); // Exact.
double hm = fputil::nearest_integer(s * y6);
+#ifdef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
+ // lo6 = 2^6 * lo.
+ double lo6_hi =
+ fputil::multiply_add(y6, e_x + LOG2_R_DD[idx_x].hi, -hm); // Exact
+ // Error bounds:
+ // | (y*log2(x) - hm * 2^-6 - lo) / y| < err(dx * p) + err(LOG2_R_DD.lo)
+ // < 2^-51 + 2^-75
+ double lo6 = fputil::multiply_add(
+ y6, fputil::multiply_add(dx, p, LOG2_R_DD[idx_x].lo), lo6_hi);
+#else
// lo6 = 2^6 * lo.
double lo6_hi =
fputil::multiply_add(y6, e_x + LOG2_R_TD[idx_x].hi, -hm); // Exact
@@ -790,6 +938,7 @@ LLVM_LIBC_FUNCTION(float, powf, (float x, float y)) {
// < 2^-51 + 2^-75
double lo6 = fputil::multiply_add(
y6, fputil::multiply_add(dx, p, LOG2_R_TD[idx_x].mid), lo6_hi);
+#endif
// |2^(hi + mid) - exp2_hi_mid| <= ulp(exp2_hi_mid) / 2
// Clamp the exponent part into smaller range that fits double precision.
@@ -830,6 +979,9 @@ LLVM_LIBC_FUNCTION(float, powf, (float x, float y)) {
double r = pp * exp2_hi_mid;
+#ifdef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
+ return static_cast<float>(r);
+#else
// Ziv accuracy test.
uint64_t r_u = cpp::bit_cast<uint64_t>(r);
float r_upper = static_cast<float>(cpp::bit_cast<double>(r_u + ERR));
@@ -861,6 +1013,7 @@ LLVM_LIBC_FUNCTION(float, powf, (float x, float y)) {
double r_dd = powf_double_double(idx_x, dx, y6, lo6_hi, exp2_hi_mid_dd);
return static_cast<float>(r_dd);
+#endif // LIBC_MATH_HAS_SKIP_ACCURATE_PASS
}
} // namespace LIBC_NAMESPACE_DECL
More information about the libc-commits
mailing list