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 LOG_R_DD[128] = { +alignas(16) const NumberPair 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 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(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(hi) << FPBits::FRACTION_LEN; + double r = + cpp::bit_cast(exp_hi + cpp::bit_cast(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(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(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(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(hi) << FPBits::FRACTION_LEN; + double r = + cpp::bit_cast(exp_hi + cpp::bit_cast(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(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(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(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(hi) << FPBits::FRACTION_LEN; + double r = + cpp::bit_cast(exp_hi + cpp::bit_cast(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(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 ziv_test_denorm(int hi, double mid, double lo, - double err) { +template +LIBC_INLINE static cpp::optional +ziv_test_denorm(int hi, double mid, double lo, double err) { using FPBits = typename fputil::FPBits; // Scaling factor = 1/(min normal number) = 2^1022 @@ -360,22 +361,27 @@ LIBC_INLINE cpp::optional 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(exp_hi + cpp::bit_cast(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(cpp::bit_cast(r) - scale_down); + } else { + double err_scaled = + cpp::bit_cast(exp_hi + cpp::bit_cast(err)); - if (LIBC_LIKELY(upper == lower)) { - return cpp::bit_cast(cpp::bit_cast(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(cpp::bit_cast(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(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(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(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(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(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; @@ -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; 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(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(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(r); +#else // Ziv accuracy test. uint64_t r_u = cpp::bit_cast(r); float r_upper = static_cast(cpp::bit_cast(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(r_dd); +#endif // LIBC_MATH_HAS_SKIP_ACCURATE_PASS } } // namespace LIBC_NAMESPACE_DECL diff --git a/libc/src/math/generic/range_reduction_double_fma.h b/libc/src/math/generic/range_reduction_double_fma.h index 8e0bc3a42462c..160fb2461fe21 100644 --- a/libc/src/math/generic/range_reduction_double_fma.h +++ b/libc/src/math/generic/range_reduction_double_fma.h @@ -80,7 +80,7 @@ LIBC_INLINE unsigned LargeRangeReduction::fast(double x, DoubleDouble &u) { // b = D(sin(k * pi/128) - a); // print("{", b, ",", a, "},"); // }; -LIBC_INLINE constexpr DoubleDouble SIN_K_PI_OVER_128[256] = { +LIBC_INLINE constexpr DoubleDouble SIN_K_PI_OVER_128[] = { {0, 0}, {-0x1.b1d63091a013p-64, 0x1.92155f7a3667ep-6}, {-0x1.912bd0d569a9p-61, 0x1.91f65f10dd814p-5}, diff --git a/libc/src/math/generic/range_reduction_double_nofma.h b/libc/src/math/generic/range_reduction_double_nofma.h index 606c3f8185d61..9d13d246ce91f 100644 --- a/libc/src/math/generic/range_reduction_double_nofma.h +++ b/libc/src/math/generic/range_reduction_double_nofma.h @@ -81,7 +81,7 @@ LIBC_INLINE unsigned LargeRangeReduction::fast(double x, DoubleDouble &u) { // b = round(sin(k * pi/128) - a, D, RN); // print("{", b, ",", a, "},"); // }; -LIBC_INLINE constexpr DoubleDouble SIN_K_PI_OVER_128[256] = { +LIBC_INLINE constexpr DoubleDouble SIN_K_PI_OVER_128[] = { {0, 0}, {0x1.f938a73db97fbp-58, 0x1.92155f7a3667cp-6}, {-0x1.912bd0d569a9p-61, 0x1.91f65f10dd814p-5},