From fee360a05dcce594afb9f997c9f2904e00f089cc Mon Sep 17 00:00:00 2001 From: Tue Ly Date: Thu, 10 Oct 2024 05:08:33 +0000 Subject: [PATCH 1/7] [libc][math] Improve performance of double precision trig functions. --- libc/src/__support/FPUtil/double_double.h | 36 +- libc/src/__support/macros/optimization.h | 10 + libc/src/math/generic/cos.cpp | 113 +++--- .../generic/range_reduction_double_common.h | 345 ++++++++++++------ .../math/generic/range_reduction_double_fma.h | 254 +++---------- .../generic/range_reduction_double_nofma.h | 253 +++---------- libc/src/math/generic/sin.cpp | 124 +++---- libc/src/math/generic/sincos.cpp | 147 ++++---- libc/src/math/generic/sincos_eval.h | 27 +- libc/src/math/generic/tan.cpp | 142 ++++--- libc/test/src/math/cos_test.cpp | 3 +- libc/test/src/math/sin_test.cpp | 12 +- libc/test/src/math/tan_test.cpp | 21 +- 13 files changed, 649 insertions(+), 838 deletions(-) diff --git a/libc/src/__support/FPUtil/double_double.h b/libc/src/__support/FPUtil/double_double.h index 25a4ee03387c6..b51d3c4c63cf5 100644 --- a/libc/src/__support/FPUtil/double_double.h +++ b/libc/src/__support/FPUtil/double_double.h @@ -73,6 +73,26 @@ template LIBC_INLINE constexpr DoubleDouble split(double a) { return r; } +// Helper for non-fma exact mult where the first number is already splitted. +template +LIBC_INLINE DoubleDouble exact_mult(const DoubleDouble &as, double a, + double b) { + DoubleDouble bs, r; + + if constexpr (NO_FMA_ALL_ROUNDINGS) + bs = split<28>(b); + else + bs = split(b); + + r.hi = a * b; + double t1 = as.hi * bs.hi - r.hi; + double t2 = as.hi * bs.lo + t1; + double t3 = as.lo * bs.hi + t2; + r.lo = as.lo * bs.lo + t3; + + return r; +} + // Note: When FMA instruction is not available, the `exact_mult` function is // only correct for round-to-nearest mode. See: // Zimmermann, P., "Note on the Veltkamp/Dekker Algorithms with Directed @@ -90,18 +110,8 @@ LIBC_INLINE DoubleDouble exact_mult(double a, double b) { #else // Dekker's Product. DoubleDouble as = split(a); - DoubleDouble bs; - if constexpr (NO_FMA_ALL_ROUNDINGS) - bs = split<28>(b); - else - bs = split(b); - - r.hi = a * b; - double t1 = as.hi * bs.hi - r.hi; - double t2 = as.hi * bs.lo + t1; - double t3 = as.lo * bs.hi + t2; - r.lo = as.lo * bs.lo + t3; + r = exact_mult(as, a, b); #endif // LIBC_TARGET_CPU_HAS_FMA return r; @@ -157,8 +167,8 @@ LIBC_INLINE DoubleDouble div(const DoubleDouble &a, const DoubleDouble &b) { double e_hi = fputil::multiply_add(b.hi, -r.hi, a.hi); double e_lo = fputil::multiply_add(b.lo, -r.hi, a.lo); #else - DoubleDouble b_hi_r_hi = fputil::exact_mult(b.hi, -r.hi); - DoubleDouble b_lo_r_hi = fputil::exact_mult(b.lo, -r.hi); + DoubleDouble b_hi_r_hi = fputil::exact_mult(b.hi, -r.hi); + DoubleDouble b_lo_r_hi = fputil::exact_mult(b.lo, -r.hi); double e_hi = (a.hi + b_hi_r_hi.hi) + b_hi_r_hi.lo; double e_lo = (a.lo + b_lo_r_hi.hi) + b_lo_r_hi.lo; #endif // LIBC_TARGET_CPU_HAS_FMA diff --git a/libc/src/__support/macros/optimization.h b/libc/src/__support/macros/optimization.h index 5ffd474d35c54..ce99f1ce6e532 100644 --- a/libc/src/__support/macros/optimization.h +++ b/libc/src/__support/macros/optimization.h @@ -48,6 +48,16 @@ LIBC_INLINE constexpr bool expects_bool_condition(T value, T expected) { #ifndef LIBC_MATH #define LIBC_MATH 0 +#else + +#if ((LIBC_MATH & LIBC_MATH_SKIP_ACCURATE_PASS) != 0) +#define LIBC_MATH_HAS_SKIP_ACCURATE_PASS +#endif + +#if ((LIBC_MATH & LIBC_MATH_SMALL_TABLES) != 0) +#define LIBC_MATH_HAS_SMALL_TABLES +#endif + #endif // LIBC_MATH #endif // LLVM_LIBC_SRC___SUPPORT_MACROS_OPTIMIZATION_H diff --git a/libc/src/math/generic/cos.cpp b/libc/src/math/generic/cos.cpp index e61d800ce2dad..fa3c489aa15ec 100644 --- a/libc/src/math/generic/cos.cpp +++ b/libc/src/math/generic/cos.cpp @@ -17,17 +17,14 @@ #include "src/__support/macros/config.h" #include "src/__support/macros/optimization.h" // LIBC_UNLIKELY #include "src/__support/macros/properties/cpu_features.h" // LIBC_TARGET_CPU_HAS_FMA +#include "src/math/generic/range_reduction_double_common.h" #include "src/math/generic/sincos_eval.h" -// TODO: We might be able to improve the performance of large range reduction of -// non-FMA targets further by operating directly on 25-bit chunks of 128/pi and -// pre-split SIN_K_PI_OVER_128, but that might double the memory footprint of -// those lookup table. -#include "range_reduction_double_common.h" - -#if ((LIBC_MATH & LIBC_MATH_SKIP_ACCURATE_PASS) != 0) -#define LIBC_MATH_COS_SKIP_ACCURATE_PASS -#endif +#ifdef LIBC_TARGET_CPU_HAS_FMA +#include "range_reduction_double_fma.h" +#else +#include "range_reduction_double_nofma.h" +#endif // LIBC_TARGET_CPU_HAS_FMA namespace LIBC_NAMESPACE_DECL { @@ -42,22 +39,29 @@ LLVM_LIBC_FUNCTION(double, cos, (double x)) { DoubleDouble y; unsigned k; - generic::LargeRangeReduction range_reduction_large{}; + LargeRangeReduction range_reduction_large{}; - // |x| < 2^32 (with FMA) or |x| < 2^23 (w/o FMA) + // |x| < 2^16. if (LIBC_LIKELY(x_e < FPBits::EXP_BIAS + FAST_PASS_EXPONENT)) { - // |x| < 2^-27 - if (LIBC_UNLIKELY(x_e < FPBits::EXP_BIAS - 27)) { - // Signed zeros. - if (LIBC_UNLIKELY(x == 0.0)) - return 1.0; - - // For |x| < 2^-27, |cos(x) - 1| < |x|^2/2 < 2^-54 = ulp(1 - 2^-53)/2. - return fputil::round_result_slightly_down(1.0); + // |x| < 2^-7 + if (LIBC_UNLIKELY(x_e < FPBits::EXP_BIAS - 7)) { + // |x| < 2^-27 + if (LIBC_UNLIKELY(x_e < FPBits::EXP_BIAS - 27)) { + // Signed zeros. + if (LIBC_UNLIKELY(x == 0.0)) + return 1.0; + + // For |x| < 2^-27, |cos(x) - 1| < |x|^2/2 < 2^-54 = ulp(1 - 2^-53)/2. + return fputil::round_result_slightly_down(1.0); + } + // No range reduction needed. + k = 0; + y.lo = 0.0; + y.hi = x; + } else { + // Small range reduction. + k = range_reduction_small(x, y); } - - // // Small range reduction. - k = range_reduction_small(x, y); } else { // Inf or NaN if (LIBC_UNLIKELY(x_e > 2 * FPBits::EXP_BIAS)) { @@ -70,45 +74,34 @@ LLVM_LIBC_FUNCTION(double, cos, (double x)) { } // Large range reduction. - k = range_reduction_large.compute_high_part(x); - y = range_reduction_large.fast(); + k = range_reduction_large.fast(x, y); } DoubleDouble sin_y, cos_y; - generic::sincos_eval(y, sin_y, cos_y); + [[maybe_unused]] double err = generic::sincos_eval(y, sin_y, cos_y); // Look up sin(k * pi/128) and cos(k * pi/128) - // Memory saving versions: - - // Use 128-entry table instead: - // DoubleDouble sin_k = SIN_K_PI_OVER_128[k & 127]; - // uint64_t sin_s = static_cast((k + 128) & 128) << (63 - 7); - // sin_k.hi = FPBits(FPBits(sin_k.hi).uintval() ^ sin_s).get_val(); - // sin_k.lo = FPBits(FPBits(sin_k.hi).uintval() ^ sin_s).get_val(); - // DoubleDouble cos_k = SIN_K_PI_OVER_128[(k + 64) & 127]; - // uint64_t cos_s = static_cast((k + 64) & 128) << (63 - 7); - // cos_k.hi = FPBits(FPBits(cos_k.hi).uintval() ^ cos_s).get_val(); - // cos_k.lo = FPBits(FPBits(cos_k.hi).uintval() ^ cos_s).get_val(); - - // Use 64-entry table instead: - // auto get_idx_dd = [](unsigned kk) -> DoubleDouble { - // unsigned idx = (kk & 64) ? 64 - (kk & 63) : (kk & 63); - // DoubleDouble ans = SIN_K_PI_OVER_128[idx]; - // if (kk & 128) { - // ans.hi = -ans.hi; - // ans.lo = -ans.lo; - // } - // return ans; - // }; - // DoubleDouble sin_k = get_idx_dd(k + 128); - // DoubleDouble cos_k = get_idx_dd(k + 64); - +#ifdef LIBC_MATH_HAS_SMALL_TABLES + // Memory saving versions. Use 65-entry table. + auto get_idx_dd = [](unsigned kk) -> DoubleDouble { + unsigned idx = (kk & 64) ? 64 - (kk & 63) : (kk & 63); + DoubleDouble ans = SIN_K_PI_OVER_128[idx]; + if (kk & 128) { + ans.hi = -ans.hi; + ans.lo = -ans.lo; + } + return ans; + }; + DoubleDouble sin_k = get_idx_dd(k + 128); + DoubleDouble cos_k = get_idx_dd(k + 64); +#else // Fast look up version, but needs 256-entry table. // -sin(k * pi/128) = sin((k + 128) * pi/128) // cos(k * pi/128) = sin(k * pi/128 + pi/2) = sin((k + 64) * pi/128). DoubleDouble msin_k = SIN_K_PI_OVER_128[(k + 128) & 255]; DoubleDouble cos_k = SIN_K_PI_OVER_128[(k + 64) & 255]; +#endif // LIBC_MATH_HAS_SMALL_TABLES // After range reduction, k = round(x * 128 / pi) and y = x - k * (pi / 128). // So k is an integer and -pi / 256 <= y <= pi / 256. @@ -120,20 +113,12 @@ LLVM_LIBC_FUNCTION(double, cos, (double x)) { DoubleDouble rr = fputil::exact_add(cos_k_cos_y.hi, msin_k_sin_y.hi); rr.lo += msin_k_sin_y.lo + cos_k_cos_y.lo; -#ifdef LIBC_MATH_COS_SKIP_ACCURATE_PASS +#ifdef LIBC_MATH_HAS_SKIP_ACCURATE_PASS return rr.hi + rr.lo; #else - // Accurate test and pass for correctly rounded implementation. -#ifdef LIBC_TARGET_CPU_HAS_FMA - constexpr double ERR = 0x1.0p-70; -#else - // TODO: Improve non-FMA fast pass accuracy. - constexpr double ERR = 0x1.0p-66; -#endif // LIBC_TARGET_CPU_HAS_FMA - - double rlp = rr.lo + ERR; - double rlm = rr.lo - ERR; + double rlp = rr.lo + err; + double rlm = rr.lo - err; double r_upper = rr.hi + rlp; // (rr.lo + ERR); double r_lower = rr.hi + rlm; // (rr.lo - ERR); @@ -144,7 +129,7 @@ LLVM_LIBC_FUNCTION(double, cos, (double x)) { Float128 u_f128, sin_u, cos_u; if (LIBC_LIKELY(x_e < FPBits::EXP_BIAS + FAST_PASS_EXPONENT)) - u_f128 = generic::range_reduction_small_f128(x); + u_f128 = range_reduction_small_f128(x); else u_f128 = range_reduction_large.accurate(); @@ -152,7 +137,7 @@ LLVM_LIBC_FUNCTION(double, cos, (double x)) { auto get_sin_k = [](unsigned kk) -> Float128 { unsigned idx = (kk & 64) ? 64 - (kk & 63) : (kk & 63); - Float128 ans = generic::SIN_K_PI_OVER_128_F128[idx]; + Float128 ans = SIN_K_PI_OVER_128_F128[idx]; if (kk & 128) ans.sign = Sign::NEG; return ans; @@ -172,7 +157,7 @@ LLVM_LIBC_FUNCTION(double, cos, (double x)) { // https://github.com/llvm/llvm-project/issues/96452. return static_cast(r); -#endif // !LIBC_MATH_COS_SKIP_ACCURATE_PASS +#endif // !LIBC_MATH_HAS_SKIP_ACCURATE_PASS } } // namespace LIBC_NAMESPACE_DECL diff --git a/libc/src/math/generic/range_reduction_double_common.h b/libc/src/math/generic/range_reduction_double_common.h index 290b642be4c69..9469ecebc5580 100644 --- a/libc/src/math/generic/range_reduction_double_common.h +++ b/libc/src/math/generic/range_reduction_double_common.h @@ -17,150 +17,269 @@ #include "src/__support/common.h" #include "src/__support/integer_literals.h" #include "src/__support/macros/config.h" +#include "src/__support/macros/optimization.h" -#ifdef LIBC_TARGET_CPU_HAS_FMA -#include "range_reduction_double_fma.h" - -// With FMA, we limit the maxmimum exponent to be 2^16, so that the error bound -// from the fma::range_reduction_small is bounded by 2^-88 instead of 2^-72. -#define FAST_PASS_EXPONENT 16 -using LIBC_NAMESPACE::fma::ONE_TWENTY_EIGHT_OVER_PI; -using LIBC_NAMESPACE::fma::range_reduction_small; -using LIBC_NAMESPACE::fma::SIN_K_PI_OVER_128; +namespace LIBC_NAMESPACE_DECL { -LIBC_INLINE constexpr bool NO_FMA = false; +#ifdef LIBC_TARGET_CPU_HAS_FMA +static constexpr bool NO_FMA = false; #else -#include "range_reduction_double_nofma.h" +static constexpr bool NO_FMA = true; +#endif // LIBC_TARGET_CPU_HAS_FMA -using LIBC_NAMESPACE::nofma::FAST_PASS_EXPONENT; -using LIBC_NAMESPACE::nofma::ONE_TWENTY_EIGHT_OVER_PI; -using LIBC_NAMESPACE::nofma::range_reduction_small; -using LIBC_NAMESPACE::nofma::SIN_K_PI_OVER_128; +using LIBC_NAMESPACE::fputil::DoubleDouble; +using Float128 = LIBC_NAMESPACE::fputil::DyadicFloat<128>; -LIBC_INLINE constexpr bool NO_FMA = true; -#endif // LIBC_TARGET_CPU_HAS_FMA +#define FAST_PASS_EXPONENT 16 -namespace LIBC_NAMESPACE_DECL { +// For 2^-7 < |x| < 2^16, return k and u such that: +// k = round(x * 128/pi) +// x mod pi/128 = x - k * pi/128 ~ u.hi + u.lo +// Error bound: +// |(x - k * pi/128) - (u_hi + u_lo)| <= max(ulp(ulp(u_hi)), 2^-119) +// <= 2^-111. +LIBC_INLINE unsigned range_reduction_small(double x, DoubleDouble &u) { + // Values of -pi/128 used for inputs with absolute value <= 2^16. + // The first 3 parts are generated with (53 - 21 = 32)-bit precision, so that + // the product k * MPI_OVER_128[i] is exact. + // Generated by Sollya with: + // > display = hexadecimal!; + // > a = round(pi/128, 32, RN); + // > b = round(pi/128 - a, 32, RN); + // > c = round(pi/128 - a - b, D, RN); + // > print(-a, ",", -b, ",", -c); + constexpr double MPI_OVER_128[3] = {-0x1.921fb544p-6, -0x1.0b4611a6p-40, + -0x1.3198a2e037073p-75}; + constexpr double ONE_TWENTY_EIGHT_OVER_PI_D = 0x1.45f306dc9c883p5; + double prod_hi = x * ONE_TWENTY_EIGHT_OVER_PI_D; + double kd = fputil::nearest_integer(prod_hi); -namespace generic { + // Let y = x - k * (pi/128) + // Then |y| < pi / 256 + // With extra rounding errors, we can bound |y| < 1.6 * 2^-7. + double y_hi = fputil::multiply_add(kd, MPI_OVER_128[0], x); // Exact + // |u.hi| < 1.6*2^-7 + u.hi = fputil::multiply_add(kd, MPI_OVER_128[1], y_hi); + double u0 = y_hi - u.hi; // Exact; + // |u.lo| <= max(ulp(u.hi), |kd * MPI_OVER_128[2]|) + double u1 = fputil::multiply_add(kd, MPI_OVER_128[1], u0); // Exact + u.lo = fputil::multiply_add(kd, MPI_OVER_128[2], u1); + // Error bound: + // |x - k * pi/128| - (u.hi + u.lo) <= ulp(u.lo) + // <= ulp(max(ulp(u.hi), kd*MPI_OVER_128[2])) + // <= 2^(-7 - 104) = 2^-111. -using LIBC_NAMESPACE::fputil::DoubleDouble; -using Float128 = LIBC_NAMESPACE::fputil::DyadicFloat<128>; + return static_cast(static_cast(kd)); +} -LIBC_INLINE constexpr Float128 PI_OVER_128_F128 = { - Sign::POS, -133, 0xc90f'daa2'2168'c234'c4c6'628b'80dc'1cd1_u128}; +// Digits of 2^(16*i) / pi, generated by Sollya with: +// > procedure ulp(x, n) { return 2^(floor(log2(abs(x))) - n); }; +// > for i from 0 to 63 do { +// if i < 3 then { pi_inv = 0.25 + 2^(16*(i - 3)) / pi; } +// else { pi_inv = 2^(16*(i-3)) / pi; }; +// pn = nearestint(pi_inv); +// pi_frac = pi_inv - pn; +// a = round(pi_frac, 51, RN); +// b = round(pi_frac - a, 51, RN); +// c = round(pi_frac - a - b, 51, RN); +// d = round(pi_frac - a - b - c, D, RN); +// print("{", 2^7 * a, ",", 2^7 * b, ",", 2^7 * c, ",", 2^7 * d, "},"); +// }; +// +// Notice that for [0..2] the leading bit of 2^(16*(i - 3)) / pi is very small, +// so we add 0.25 so that the conditions for the algorithms are still satisfied, +// and one of those conditions guarantees that ulp(0.25 * x_reduced) >= 2, and +// will safely be discarded. -// Note: The look-up tables ONE_TWENTY_EIGHT_OVER_PI is selected to be either -// from fma:: or nofma:: namespace. +static constexpr double ONE_TWENTY_EIGHT_OVER_PI[64][4] = { + {0x1.0000000000014p5, 0x1.7cc1b727220a8p-49, 0x1.4fe13abe8fa9cp-101, + -0x1.911f924eb5336p-153}, + {0x1.0000000145f3p5, 0x1.b727220a94fep-49, 0x1.3abe8fa9a6eep-101, + 0x1.b6c52b3278872p-155}, + {0x1.000145f306dc8p5, 0x1.c882a53f84ebp-47, -0x1.70565911f925p-101, + 0x1.4acc9e21c821p-153}, + {0x1.45f306dc9c884p5, -0x1.5ac07b1505c14p-47, -0x1.96447e493ad4cp-99, + -0x1.b0ef1bef806bap-152}, + {-0x1.f246c6efab58p4, -0x1.ec5417056591p-49, -0x1.f924eb53361ep-101, + 0x1.c820ff28b1d5fp-153}, + {0x1.391054a7f09d4p4, 0x1.f47d4d377036cp-48, 0x1.8a5664f10e41p-100, + 0x1.fe5163abdebbcp-154}, + {0x1.529fc2757d1f4p2, 0x1.34ddc0db62958p-50, 0x1.93c439041fe5p-102, + 0x1.63abdebbc561bp-154}, + {-0x1.ec5417056591p-1, -0x1.f924eb53361ep-53, 0x1.c820ff28b1d6p-105, + -0x1.0a21d4f246dc9p-157}, + {-0x1.505c1596447e4p5, -0x1.275a99b0ef1cp-48, 0x1.07f9458eaf7bp-100, + -0x1.0ea79236e4717p-152}, + {-0x1.596447e493ad4p1, -0x1.9b0ef1bef806cp-52, 0x1.63abdebbc561cp-106, + -0x1.1b7238b7b645ap-159}, + {0x1.bb81b6c52b328p5, -0x1.de37df00d74e4p-49, 0x1.5ef5de2b0db94p-101, + -0x1.c8e2ded9169p-153}, + {0x1.b6c52b3278874p5, -0x1.f7c035d38a844p-47, 0x1.778ac36e48dc8p-99, + -0x1.6f6c8b47fe6dbp-152}, + {0x1.2b3278872084p5, -0x1.ae9c5421443a8p-50, -0x1.e48db91c5bdb4p-102, + 0x1.d2e006492eea1p-154}, + {-0x1.8778df7c035d4p5, 0x1.d5ef5de2b0db8p-49, 0x1.2371d2126e97p-101, + 0x1.924bba8274648p-160}, + {-0x1.bef806ba71508p4, -0x1.443a9e48db91cp-50, -0x1.6f6c8b47fe6dcp-104, + 0x1.77504e8c90e7fp-157}, + {-0x1.ae9c5421443a8p-2, -0x1.e48db91c5bdb4p-54, 0x1.d2e006492eeap-106, + 0x1.3a32439fc3bd6p-159}, + {-0x1.38a84288753c8p5, -0x1.1b7238b7b645cp-47, 0x1.c00c925dd413cp-99, + -0x1.cdbc603c429c7p-151}, + {-0x1.0a21d4f246dc8p3, -0x1.c5bdb22d1ff9cp-50, 0x1.25dd413a32438p-103, + 0x1.fc3bd63962535p-155}, + {-0x1.d4f246dc8e2ep3, 0x1.26e9700324978p-49, -0x1.5f62e6de301e4p-102, + 0x1.eb1cb129a73efp-154}, + {-0x1.236e4716f6c8cp4, 0x1.700324977505p-49, -0x1.736f180f10a7p-101, + -0x1.a76b2c608bbeep-153}, + {0x1.b8e909374b8p4, 0x1.924bba8274648p-48, 0x1.cfe1deb1cb128p-102, + 0x1.a73ee88235f53p-154}, + {0x1.09374b801924cp4, -0x1.15f62e6de302p-50, 0x1.deb1cb129a74p-102, + -0x1.177dca0ad144cp-154}, + {-0x1.68ffcdb688afcp3, 0x1.d1921cfe1debp-50, 0x1.cb129a73ee884p-102, + -0x1.ca0ad144bb7b1p-154}, + {0x1.924bba8274648p0, 0x1.cfe1deb1cb128p-54, 0x1.a73ee88235f54p-106, + -0x1.144bb7b16639p-158}, + {-0x1.a22bec5cdbc6p5, -0x1.e214e34ed658cp-50, -0x1.177dca0ad144cp-106, + 0x1.213a671c09ad1p-160}, + {0x1.3a32439fc3bd8p1, -0x1.c69dacb1822fp-51, 0x1.1afa975da2428p-105, + -0x1.6638fd94ba082p-158}, + {-0x1.b78c0788538d4p4, 0x1.29a73ee88236p-50, -0x1.5a28976f62cc8p-103, + 0x1.c09ad17df904ep-156}, + {0x1.fc3bd63962534p5, 0x1.cfba208d7d4bcp-48, -0x1.12edec598e3f8p-100, + 0x1.ad17df904e647p-152}, + {-0x1.4e34ed658c118p2, 0x1.046bea5d7689p-51, 0x1.3a671c09ad17cp-104, + 0x1.f904e64758e61p-156}, + {0x1.62534e7dd1048p5, -0x1.415a28976f62cp-47, -0x1.8e3f652e8207p-100, + 0x1.3991d63983534p-154}, + {-0x1.63045df7282b4p4, -0x1.44bb7b16638fcp-50, -0x1.94ba081bec67p-102, + 0x1.d639835339f4ap-154}, + {0x1.d1046bea5d768p5, 0x1.213a671c09adp-48, 0x1.7df904e64759p-100, + -0x1.9f2b3182d8defp-152}, + {0x1.afa975da24274p3, 0x1.9c7026b45f7e4p-50, 0x1.3991d63983534p-106, + -0x1.82d8dee81d108p-160}, + {-0x1.a28976f62cc7p5, -0x1.fb29741037d8cp-47, -0x1.b8a719f2b3184p-100, + 0x1.272117e2ef7e5p-152}, + {-0x1.76f62cc71fb28p5, -0x1.741037d8cdc54p-47, 0x1.cc1a99cfa4e44p-101, + -0x1.d03a21036be27p-153}, + {0x1.d338e04d68bfp5, -0x1.bec66e29c67ccp-50, 0x1.339f49c845f8cp-102, + -0x1.081b5f13801dap-156}, + {0x1.c09ad17df905p4, -0x1.9b8a719f2b318p-48, -0x1.6c6f740e8840cp-103, + -0x1.af89c00ed0004p-155}, + {0x1.68befc827323cp5, -0x1.38cf9598c16c8p-47, 0x1.08bf177bf2508p-99, + -0x1.3801da00087eap-152}, + {-0x1.037d8cdc538dp5, 0x1.a99cfa4e422fcp-49, 0x1.77bf250763ffp-103, + 0x1.2fffbc0b301fep-155}, + {-0x1.8cdc538cf9598p5, -0x1.82d8dee81d108p-48, -0x1.b5f13801dap-104, + -0x1.0fd33f8086877p-157}, + {-0x1.4e33e566305bp3, -0x1.bdd03a21036cp-49, 0x1.d8ffc4bffef04p-101, + -0x1.33f80868773a5p-153}, + {-0x1.f2b3182d8dee8p4, -0x1.d1081b5f138p-52, -0x1.da00087e99fcp-104, + -0x1.0d0ee74a5f593p-158}, + {-0x1.8c16c6f740e88p5, -0x1.036be27003b4p-49, -0x1.0fd33f8086878p-109, + 0x1.8b5a0a6d1f6d3p-162}, + {0x1.3908bf177bf24p5, 0x1.0763ff12fffbcp-47, 0x1.6603fbcbc462cp-104, + 0x1.6829b47db4dap-156}, + {0x1.7e2ef7e4a0ec8p4, -0x1.da00087e99fcp-56, -0x1.0d0ee74a5f594p-110, + 0x1.1f6d367ecf27dp-162}, + {-0x1.081b5f13801dcp4, 0x1.fff7816603fbcp-48, 0x1.788c5ad05369p-101, + -0x1.25930261b069fp-155}, + {-0x1.af89c00ed0004p5, -0x1.fa67f010d0ee8p-50, 0x1.6b414da3eda6cp-103, + 0x1.fb3c9f2c26dd4p-156}, + {-0x1.c00ed00043f4cp5, -0x1.fc04343b9d298p-48, 0x1.4da3eda6cfdap-103, + -0x1.b069ec9161738p-155}, + {0x1.2fffbc0b301fcp5, 0x1.e5e2316b414dcp-47, -0x1.c125930261b08p-99, + 0x1.6136e9e8c7ecdp-151}, + {-0x1.0fd33f8086878p3, 0x1.8b5a0a6d1f6d4p-50, -0x1.30261b069ec9p-103, + -0x1.61738132c3403p-155}, + {-0x1.9fc04343b9d28p4, -0x1.7d64b824b2604p-48, -0x1.86c1a7b24585cp-101, + -0x1.c09961a015d29p-154}, + {-0x1.0d0ee74a5f594p2, 0x1.1f6d367ecf27cp-50, 0x1.6136e9e8c7eccp-103, + 0x1.3cbfd45aea4f7p-155}, + {-0x1.dce94beb25c14p5, 0x1.a6cfd9e4f9614p-47, -0x1.22c2e70265868p-100, + -0x1.5d28ad8453814p-158}, + {-0x1.4beb25c12593p5, -0x1.30d834f648b0cp-50, 0x1.8fd9a797fa8b4p-104, + 0x1.d49eeb1faf97cp-156}, + {0x1.b47db4d9fb3c8p4, 0x1.f2c26dd3d18fcp-48, 0x1.9a797fa8b5d48p-100, + 0x1.eeb1faf97c5edp-152}, + {-0x1.25930261b06ap5, 0x1.36e9e8c7ecd3cp-47, 0x1.7fa8b5d49eebp-100, + 0x1.faf97c5ecf41dp-152}, + {0x1.fb3c9f2c26dd4p4, -0x1.738132c3402bcp-51, 0x1.aea4f758fd7ccp-103, + -0x1.d0985f18c10ebp-159}, + {-0x1.b069ec9161738p5, -0x1.32c3402ba515cp-51, 0x1.eeb1faf97c5ecp-104, + 0x1.e839cfbc52949p-157}, + {-0x1.ec9161738132cp5, -0x1.a015d28ad8454p-50, 0x1.faf97c5ecf41cp-104, + 0x1.cfbc529497536p-157}, + {-0x1.61738132c3404p5, 0x1.45aea4f758fd8p-47, -0x1.a0e84c2f8c608p-102, + -0x1.d6b5b45650128p-156}, + {0x1.fb34f2ff516bcp3, -0x1.6c229c0a0d074p-49, -0x1.30be31821d6b4p-104, + -0x1.b4565012813b8p-156}, + {0x1.3cbfd45aea4f8p5, -0x1.4e050683a130cp-48, 0x1.ce7de294a4ba8p-104, + 0x1.afed7ec47e357p-156}, + {-0x1.5d28ad8453814p2, -0x1.a0e84c2f8c608p-54, -0x1.d6b5b45650128p-108, + -0x1.3b81ca8bdea7fp-164}, + {-0x1.15b08a702834p5, -0x1.d0985f18c10ecp-47, 0x1.4a4ba9afed7ecp-100, + 0x1.1f8d5d0856033p-154}, +}; -// For large range |x| >= 2^32, we use the exponent of x to find 3 double-chunks -// of 128/pi c_hi, c_mid, c_lo such that: -// 1) ulp(round(x * c_hi, D, RN)) >= 256, +// For large range |x| >= 2^16, we perform the range reduction computations as: +// u = x - k * pi/128 = (pi/128) * (x * (128/pi) - k). +// We use the exponent of x to find 4 double-chunks of 128/pi: +// c_hi, c_mid, c_lo, c_lo_2 such that: +// 1) ulp(round(x * c_hi, D, RN)) >= 2^8 = 256, // 2) If x * c_hi = ph_hi + ph_lo and x * c_mid = pm_hi + pm_lo, then // min(ulp(ph_lo), ulp(pm_hi)) >= 2^-53. -// 3) ulp(round(x * c_lo, D, RN)) <= 2^-7x. -// This will allow us to do quick computations as: -// (x * 256/pi) ~ x * (c_hi + c_mid + c_lo) (mod 256) -// ~ ph_lo + pm_hi + pm_lo + (x * c_lo) +// This will allow us to drop the high part ph_hi and the addition: +// (ph_lo + pm_hi) mod 1 +// can be exactly representable in a double precision. +// This will allow us to do split the computations as: +// (x * 256/pi) ~ x * (c_hi + c_mid + c_lo + c_lo_2) (mod 256) +// ~ (ph_lo + pm_hi) + (pm_lo + x * c_lo) + x * c_lo_2. // Then, // round(x * 128/pi) = round(ph_lo + pm_hi) (mod 256) // And the high part of fractional part of (x * 128/pi) can simply be: // {x * 128/pi}_hi = {ph_lo + pm_hi}. // To prevent overflow when x is very large, we simply scale up -// (c_hi, c_mid, c_lo) by a fixed power of 2 (based on the index) and scale down -// x by the same amount. - -template struct LargeRangeReduction { - // Calculate the high part of the range reduction exactly. - LIBC_INLINE unsigned compute_high_part(double x) { - using FPBits = typename fputil::FPBits; - FPBits xbits(x); - - // TODO: The extra exponent gap of 62 below can be reduced a bit for non-FMA - // with a more careful analysis, which in turn will reduce the error bound - // for non-FMA - int x_e_m62 = xbits.get_biased_exponent() - (FPBits::EXP_BIAS + 62); - idx = static_cast((x_e_m62 >> 4) + 3); - // Scale x down by 2^(-(16 * (idx - 3)) - xbits.set_biased_exponent((x_e_m62 & 15) + FPBits::EXP_BIAS + 62); - // 2^62 <= |x_reduced| < 2^(62 + 16) = 2^78 - x_reduced = xbits.get_val(); - // x * c_hi = ph.hi + ph.lo exactly. - DoubleDouble ph = - fputil::exact_mult(x_reduced, ONE_TWENTY_EIGHT_OVER_PI[idx][0]); - // x * c_mid = pm.hi + pm.lo exactly. - DoubleDouble pm = - fputil::exact_mult(x_reduced, ONE_TWENTY_EIGHT_OVER_PI[idx][1]); - // Extract integral parts and fractional parts of (ph.lo + pm.hi). - double kh = fputil::nearest_integer(ph.lo); - double ph_lo_frac = ph.lo - kh; // Exact - double km = fputil::nearest_integer(pm.hi + ph_lo_frac); - double pm_hi_frac = pm.hi - km; // Exact - // x * 128/pi mod 1 ~ y_hi + y_lo - y_hi = ph_lo_frac + pm_hi_frac; // Exact - pm_lo = pm.lo; - return static_cast(static_cast(kh) + - static_cast(km)); - } +// (c_hi, c_mid, c_lo, c_lo_2) by a fixed power of 2 (based on the index) and +// scale down x by the same amount. - LIBC_INLINE DoubleDouble fast() const { - // y_lo = x * c_lo + pm.lo - double y_lo = fputil::multiply_add(x_reduced, - ONE_TWENTY_EIGHT_OVER_PI[idx][2], pm_lo); - DoubleDouble y = fputil::exact_add(y_hi, y_lo); - - // Digits of pi/128, generated by Sollya with: - // > a = round(pi/128, D, RN); - // > b = round(pi/128 - a, D, RN); - constexpr DoubleDouble PI_OVER_128_DD = {0x1.1a62633145c07p-60, - 0x1.921fb54442d18p-6}; - - // Error bound: with {a} denote the fractional part of a, i.e.: - // {a} = a - round(a) - // Then, - // | {x * 128/pi} - (y_hi + y_lo) | < 2 * ulp(x_reduced * - // * ONE_TWENTY_EIGHT_OVER_PI[idx][2]) - // For FMA: - // | {x * 128/pi} - (y_hi + y_lo) | <= 2 * 2^77 * 2^-103 * 2^-52 - // = 2^-77. - // | {x mod pi/128} - (u.hi + u.lo) | < 2 * 2^-6 * 2^-77. - // = 2^-82. - // For non-FMA: - // | {x * 128/pi} - (y_hi + y_lo) | <= 2 * 2^77 * 2^-99 * 2^-52 - // = 2^-73. - // | {x mod pi/128} - (u.hi + u.lo) | < 2 * 2^-6 * 2^-73. - // = 2^-78. - return fputil::quick_mult(y, PI_OVER_128_DD); - } +struct LargeRangeReduction { + + // To be implemented in range_reduction_double_fma.h and + // range_reduction_double_nofma.h. + unsigned fast(double x, DoubleDouble &u); +#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS LIBC_INLINE Float128 accurate() const { + constexpr Float128 PI_OVER_128_F128 = { + Sign::POS, -133, 0xc90f'daa2'2168'c234'c4c6'628b'80dc'1cd1_u128}; + // y_lo = x * c_lo + pm.lo Float128 y_lo_0(x_reduced * ONE_TWENTY_EIGHT_OVER_PI[idx][3]); - Float128 y_lo_1 = fputil::quick_mul( - Float128(x_reduced), Float128(ONE_TWENTY_EIGHT_OVER_PI[idx][2])); - Float128 y_lo_2(pm_lo); - Float128 y_hi_f128(y_hi); - - Float128 y = fputil::quick_add( - y_hi_f128, - fputil::quick_add(y_lo_2, fputil::quick_add(y_lo_1, y_lo_0))); + Float128 y_lo_1 = fputil::quick_add(Float128(y_lo), y_lo_0); + Float128 y_mid_f128 = fputil::quick_add(Float128(y_mid.lo), y_lo_1); + Float128 y_hi_f128 = fputil::quick_add(Float128(y_hi), Float128(y_mid.hi)); + Float128 y = fputil::quick_add(y_hi_f128, y_mid_f128); return fputil::quick_mul(y, PI_OVER_128_F128); } +#endif // !LIBC_MATH_HAS_SKIP_ACCURATE_PASS private: // Index of x in the look-up table ONE_TWENTY_EIGHT_OVER_PI. unsigned idx; // x scaled down by 2^(-16 *(idx - 3))). double x_reduced; - // High part of (x * 128/pi) mod 1. - double y_hi; - // Low part of x * ONE_TWENTY_EIGHT_OVER_PI[idx][1]. - double pm_lo; + // Parts of (x * 128/pi) mod 1. + double y_hi, y_lo; + DoubleDouble y_mid; }; -LIBC_INLINE Float128 range_reduction_small_f128(double x) { - double prod_hi = x * ONE_TWENTY_EIGHT_OVER_PI[3][0]; +static Float128 range_reduction_small_f128(double x) { + constexpr Float128 PI_OVER_128_F128 = { + Sign::POS, -133, 0xc90f'daa2'2168'c234'c4c6'628b'80dc'1cd1_u128}; + constexpr double ONE_TWENTY_EIGHT_OVER_PI_D = 0x1.45f306dc9c883p5; + double prod_hi = x * ONE_TWENTY_EIGHT_OVER_PI_D; double kd = fputil::nearest_integer(prod_hi); Float128 mk_f128(-kd); @@ -178,7 +297,8 @@ LIBC_INLINE Float128 range_reduction_small_f128(double x) { return fputil::quick_mul(y, PI_OVER_128_F128); } -LIBC_INLINE constexpr Float128 SIN_K_PI_OVER_128_F128[65] = { +#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS +static constexpr Float128 SIN_K_PI_OVER_128_F128[65] = { {Sign::POS, 0, 0}, {Sign::POS, -133, 0xc90a'afbd'1b33'efc9'c539'edcb'fda0'cf2c_u128}, {Sign::POS, -132, 0xc8fb'2f88'6ec0'9f37'6a17'954b'2b7c'5171_u128}, @@ -245,8 +365,7 @@ LIBC_INLINE constexpr Float128 SIN_K_PI_OVER_128_F128[65] = { {Sign::POS, -128, 0xffec'4304'2668'65d9'5657'5523'6696'1732_u128}, {Sign::POS, 0, 1}, }; - -} // namespace generic +#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 7448b5f63dfde..d220bef4cd349 100644 --- a/libc/src/math/generic/range_reduction_double_fma.h +++ b/libc/src/math/generic/range_reduction_double_fma.h @@ -15,174 +15,62 @@ #include "src/__support/FPUtil/nearest_integer.h" #include "src/__support/common.h" #include "src/__support/macros/config.h" +#include "src/__support/macros/optimization.h" +#include "src/math/generic/range_reduction_double_common.h" namespace LIBC_NAMESPACE_DECL { -namespace fma { - using LIBC_NAMESPACE::fputil::DoubleDouble; -LIBC_INLINE constexpr int FAST_PASS_EXPONENT = 32; +LIBC_INLINE unsigned LargeRangeReduction::fast(double x, DoubleDouble &u) { + using FPBits = typename fputil::FPBits; + FPBits xbits(x); -// Digits of 2^(16*i) / pi, generated by Sollya with: -// For [2..62]: -// > for i from 3 to 63 do { -// pi_inv = 2^(16*(i - 3)) / pi; -// pn = nearestint(pi_inv); -// pi_frac = pi_inv - pn; -// a = round(pi_frac, D, RN); -// b = round(pi_frac - a, D, RN); -// c = round(pi_frac - a - b, D, RN); -// d = round(pi_frac - a - b - c, D, RN); -// print("{", 2^7 * a, ",", 2^7 * b, ",", 2^7 * c, ",", 2^7 * d, "},"); -// }; -// For [0..1]: -// The leading bit of 2^(16*(i - 3)) / pi is very small, so we add 0.25 so that -// the conditions for the algorithms are still satisfied, and one of those -// conditions guarantees that ulp(0.25 * x_reduced) >= 2, and will safely be -// discarded. -// for i from 0 to 2 do { -// pi_frac = 0.25 + 2^(16*(i - 3)) / pi; -// a = round(pi_frac, D, RN); -// b = round(pi_frac - a, D, RN); -// c = round(pi_frac - a - b, D, RN); -// d = round(pi_frac - a - b - c, D, RN); -// print("{", 2^7 * a, ",", 2^7 * b, ",", 2^7 * c, ",", 2^7 * d, "},"); -// }; -// For The fast pass using double-double, we only need 3 parts (a, b, c), but -// for the accurate pass using Float128, instead of using another table of -// Float128s, we simply add the fourth path (a, b, c, d), which simplify the -// implementation a bit and saving some memory. -LIBC_INLINE constexpr double ONE_TWENTY_EIGHT_OVER_PI[64][4] = { - {0x1.0000000000014p5, 0x1.7cc1b727220a9p-49, 0x1.3f84eafa3ea6ap-103, - -0x1.11f924eb53362p-157}, - {0x1.0000000145f3p5, 0x1.b727220a94fe1p-49, 0x1.d5f47d4d37703p-104, - 0x1.b6295993c439p-158}, - {0x1.000145f306dcap5, -0x1.bbead603d8a83p-50, 0x1.f534ddc0db629p-106, - 0x1.664f10e4107f9p-160}, - {0x1.45f306dc9c883p5, -0x1.6b01ec5417056p-49, -0x1.6447e493ad4cep-103, - 0x1.e21c820ff28b2p-157}, - {-0x1.f246c6efab581p4, 0x1.3abe8fa9a6eep-53, 0x1.b6c52b3278872p-107, - 0x1.07f9458eaf7afp-164}, - {0x1.391054a7f09d6p4, -0x1.70565911f924fp-53, 0x1.2b3278872084p-107, - -0x1.ae9c5421443aap-162}, - {0x1.529fc2757d1f5p2, 0x1.a6ee06db14acdp-53, -0x1.8778df7c035d4p-107, - 0x1.d5ef5de2b0db9p-161}, - {-0x1.ec54170565912p-1, 0x1.b6c52b3278872p-59, 0x1.07f9458eaf7afp-116, - -0x1.d4f246dc8e2dfp-173}, - {-0x1.505c1596447e5p5, 0x1.b14acc9e21c82p-49, 0x1.fe5163abdebbcp-106, - 0x1.586dc91b8e909p-160}, - {-0x1.596447e493ad5p1, 0x1.93c439041fe51p-54, 0x1.8eaf7aef1586ep-108, - -0x1.b7238b7b645a4p-163}, - {0x1.bb81b6c52b328p5, -0x1.de37df00d74e3p-49, 0x1.7bd778ac36e49p-103, - -0x1.1c5bdb22d1ffap-158}, - {0x1.b6c52b3278872p5, 0x1.07f9458eaf7afp-52, -0x1.d4f246dc8e2dfp-109, - 0x1.374b801924bbbp-164}, - {0x1.2b3278872084p5, -0x1.ae9c5421443aap-50, 0x1.b7246e3a424ddp-106, - 0x1.700324977504fp-161}, - {-0x1.8778df7c035d4p5, 0x1.d5ef5de2b0db9p-49, 0x1.1b8e909374b8p-104, - 0x1.924bba8274648p-160}, - {-0x1.bef806ba71508p4, -0x1.443a9e48db91cp-50, -0x1.6f6c8b47fe6dbp-104, - -0x1.115f62e6de302p-158}, - {-0x1.ae9c5421443aap-2, 0x1.b7246e3a424ddp-58, 0x1.700324977504fp-113, - -0x1.cdbc603c429c7p-167}, - {-0x1.38a84288753c9p5, -0x1.b7238b7b645a4p-51, 0x1.924bba8274648p-112, - 0x1.cfe1deb1cb12ap-166}, - {-0x1.0a21d4f246dc9p3, 0x1.d2126e9700325p-53, -0x1.a22bec5cdbc6p-107, - -0x1.e214e34ed658cp-162}, - {-0x1.d4f246dc8e2dfp3, 0x1.374b801924bbbp-52, -0x1.f62e6de301e21p-106, - -0x1.38d3b5963045ep-160}, - {-0x1.236e4716f6c8bp4, -0x1.1ff9b6d115f63p-50, 0x1.921cfe1deb1cbp-106, - 0x1.29a73ee88235fp-162}, - {0x1.b8e909374b802p4, -0x1.b6d115f62e6dep-50, -0x1.80f10a71a76b3p-105, - 0x1.cfba208d7d4bbp-160}, - {0x1.09374b801924cp4, -0x1.15f62e6de301ep-50, -0x1.0a71a76b2c609p-105, - 0x1.1046bea5d7689p-159}, - {-0x1.68ffcdb688afbp3, -0x1.736f180f10a72p-53, 0x1.62534e7dd1047p-107, - -0x1.0568a25dbd8b3p-161}, - {0x1.924bba8274648p0, 0x1.cfe1deb1cb12ap-54, -0x1.63045df7282b4p-108, - -0x1.44bb7b16638fep-162}, - {-0x1.a22bec5cdbc6p5, -0x1.e214e34ed658cp-50, -0x1.177dca0ad144cp-106, - 0x1.213a671c09ad1p-160}, - {0x1.3a32439fc3bd6p1, 0x1.cb129a73ee882p-54, 0x1.afa975da24275p-109, - -0x1.8e3f652e8207p-164}, - {-0x1.b78c0788538d4p4, 0x1.29a73ee88235fp-50, 0x1.4baed1213a672p-104, - -0x1.fb29741037d8dp-159}, - {0x1.fc3bd63962535p5, -0x1.822efb9415a29p-51, 0x1.a24274ce38136p-105, - -0x1.741037d8cdc54p-159}, - {-0x1.4e34ed658c117p2, -0x1.f7282b4512edfp-52, 0x1.d338e04d68bfp-107, - -0x1.bec66e29c67cbp-162}, - {0x1.62534e7dd1047p5, -0x1.0568a25dbd8b3p-49, -0x1.c7eca5d040df6p-105, - -0x1.9b8a719f2b318p-160}, - {-0x1.63045df7282b4p4, -0x1.44bb7b16638fep-50, 0x1.ad17df904e647p-104, - 0x1.639835339f49dp-158}, - {0x1.d1046bea5d769p5, -0x1.bd8b31c7eca5dp-49, -0x1.037d8cdc538dp-107, - 0x1.a99cfa4e422fcp-161}, - {0x1.afa975da24275p3, -0x1.8e3f652e8207p-52, 0x1.3991d63983534p-106, - -0x1.82d8dee81d108p-160}, - {-0x1.a28976f62cc72p5, 0x1.35a2fbf209cc9p-53, -0x1.4e33e566305b2p-109, - 0x1.08bf177bf2507p-163}, - {-0x1.76f62cc71fb29p5, -0x1.d040df633714ep-49, -0x1.9f2b3182d8defp-104, - 0x1.f8bbdf9283b2p-158}, - {0x1.d338e04d68bfp5, -0x1.bec66e29c67cbp-50, 0x1.9cfa4e422fc5ep-105, - -0x1.036be27003b4p-161}, - {0x1.c09ad17df904ep4, 0x1.91d639835339fp-50, 0x1.272117e2ef7e5p-104, - -0x1.7c4e007680022p-158}, - {0x1.68befc827323bp5, -0x1.c67cacc60b638p-50, 0x1.17e2ef7e4a0ecp-104, - 0x1.ff897ffde0598p-158}, - {-0x1.037d8cdc538dp5, 0x1.a99cfa4e422fcp-49, 0x1.77bf250763ff1p-103, - 0x1.7ffde05980fefp-158}, - {-0x1.8cdc538cf9599p5, 0x1.f49c845f8bbep-50, -0x1.b5f13801da001p-104, - 0x1.e05980fef2f12p-158}, - {-0x1.4e33e566305b2p3, 0x1.08bf177bf2507p-51, 0x1.8ffc4bffef02dp-105, - -0x1.fc04343b9d298p-160}, - {-0x1.f2b3182d8dee8p4, -0x1.d1081b5f13802p-52, 0x1.2fffbc0b301fep-107, - -0x1.a1dce94beb25cp-163}, - {-0x1.8c16c6f740e88p5, -0x1.036be27003b4p-49, -0x1.0fd33f8086877p-109, - -0x1.d297d64b824b2p-164}, - {0x1.3908bf177bf25p5, 0x1.d8ffc4bffef03p-53, -0x1.9fc04343b9d29p-108, - -0x1.f592e092c9813p-162}, - {0x1.7e2ef7e4a0ec8p4, -0x1.da00087e99fcp-56, -0x1.0d0ee74a5f593p-110, - 0x1.f6d367ecf27cbp-166}, - {-0x1.081b5f13801dap4, -0x1.0fd33f8086877p-61, -0x1.d297d64b824b2p-116, - -0x1.8130d834f648bp-170}, - {-0x1.af89c00ed0004p5, -0x1.fa67f010d0ee7p-50, -0x1.297d64b824b26p-104, - -0x1.30d834f648b0cp-162}, - {-0x1.c00ed00043f4dp5, 0x1.fde5e2316b415p-55, -0x1.2e092c98130d8p-110, - -0x1.a7b24585ce04dp-165}, - {0x1.2fffbc0b301fep5, -0x1.a1dce94beb25cp-51, -0x1.25930261b069fp-107, - 0x1.b74f463f669e6p-162}, - {-0x1.0fd33f8086877p3, -0x1.d297d64b824b2p-52, -0x1.8130d834f648bp-106, - -0x1.738132c3402bap-163}, - {-0x1.9fc04343b9d29p4, -0x1.f592e092c9813p-50, -0x1.b069ec9161738p-107, - -0x1.32c3402ba515bp-163}, - {-0x1.0d0ee74a5f593p2, 0x1.f6d367ecf27cbp-54, 0x1.36e9e8c7ecd3dp-111, - -0x1.00ae9456c229cp-165}, - {-0x1.dce94beb25c12p5, -0x1.64c0986c1a7b2p-49, -0x1.161738132c34p-103, - -0x1.5d28ad8453814p-158}, - {-0x1.4beb25c12593p5, -0x1.30d834f648b0cp-50, 0x1.8fd9a797fa8b6p-104, - -0x1.5b08a7028341dp-159}, - {0x1.b47db4d9fb3cap4, -0x1.a7b24585ce04dp-53, 0x1.3cbfd45aea4f7p-107, - 0x1.63f5f2f8bd9e8p-161}, - {-0x1.25930261b069fp5, 0x1.b74f463f669e6p-50, -0x1.5d28ad8453814p-110, - -0x1.a0e84c2f8c608p-166}, - {0x1.fb3c9f2c26dd4p4, -0x1.738132c3402bap-51, -0x1.456c229c0a0dp-105, - -0x1.d0985f18c10ebp-159}, - {-0x1.b069ec9161738p5, -0x1.32c3402ba515bp-51, -0x1.14e050683a131p-108, - 0x1.0739f78a5292fp-162}, - {-0x1.ec9161738132cp5, -0x1.a015d28ad8454p-50, 0x1.faf97c5ecf41dp-104, - -0x1.821d6b5b4565p-160}, - {-0x1.61738132c3403p5, 0x1.16ba93dd63f5fp-49, 0x1.7c5ecf41ce7dep-104, - 0x1.4a525d4d7f6bfp-159}, - {0x1.fb34f2ff516bbp3, -0x1.b08a7028341d1p-51, 0x1.9e839cfbc5295p-105, - -0x1.a2b2809409dc1p-159}, - {0x1.3cbfd45aea4f7p5, 0x1.63f5f2f8bd9e8p-49, 0x1.ce7de294a4baap-104, - -0x1.404a04ee072a3p-158}, - {-0x1.5d28ad8453814p2, -0x1.a0e84c2f8c608p-54, -0x1.d6b5b45650128p-108, - -0x1.3b81ca8bdea7fp-164}, - {-0x1.15b08a7028342p5, 0x1.7b3d0739f78a5p-50, 0x1.497535fdafd89p-105, - -0x1.ca8bdea7f33eep-164}, -}; + int x_e_m62 = xbits.get_biased_exponent() - (FPBits::EXP_BIAS + 62); + idx = static_cast((x_e_m62 >> 4) + 3); + // Scale x down by 2^(-(16 * (idx - 3)) + xbits.set_biased_exponent((x_e_m62 & 15) + FPBits::EXP_BIAS + 62); + // 2^62 <= |x_reduced| < 2^(62 + 16) = 2^78 + x_reduced = xbits.get_val(); + // x * c_hi = ph.hi + ph.lo exactly. + DoubleDouble ph = + fputil::exact_mult(x_reduced, ONE_TWENTY_EIGHT_OVER_PI[idx][0]); + // x * c_mid = pm.hi + pm.lo exactly. + DoubleDouble pm = + fputil::exact_mult(x_reduced, ONE_TWENTY_EIGHT_OVER_PI[idx][1]); + // x * c_lo = pl.hi + pl.lo exactly. + DoubleDouble pl = + fputil::exact_mult(x_reduced, ONE_TWENTY_EIGHT_OVER_PI[idx][2]); + // Extract integral parts and fractional parts of (ph.lo + pm.hi). + double sum_hi = ph.lo + pm.hi; + double kd = fputil::nearest_integer(sum_hi); + + // x * 128/pi mod 1 ~ y_hi + y_mid + y_lo + y_hi = (ph.lo - kd) + pm.hi; // Exact + y_mid = fputil::exact_add(pm.lo, pl.hi); + y_lo = pl.lo; + + // y_l = x * c_lo_2 + pl.lo + double y_l = + fputil::multiply_add(x_reduced, ONE_TWENTY_EIGHT_OVER_PI[idx][3], y_lo); + DoubleDouble y = fputil::exact_add(y_hi, y_mid.hi); + y.lo += (y_mid.lo + y_l); + + // Digits of pi/128, generated by Sollya with: + // > a = round(pi/128, D, RN); + // > b = round(pi/128 - a, D, RN); + constexpr DoubleDouble PI_OVER_128_DD = {0x1.1a62633145c07p-60, + 0x1.921fb54442d18p-6}; + + // Error bound: with {a} denote the fractional part of a, i.e.: + // {a} = a - round(a) + // Then, + // | {x * 128/pi} - (y_hi + y_lo) | <= ulp(ulp(y_hi)) <= 2^-105 + // | {x mod pi/128} - (u.hi + u.lo) | < 2 * 2^-6 * 2^-105 = 2^-110 + u = fputil::quick_mult(y, PI_OVER_128_DD); + + return static_cast(static_cast(kd)); +} // Lookup table for sin(k * pi / 128) with k = 0, ..., 255. // Table is generated with Sollya as follow: @@ -258,6 +146,7 @@ LIBC_INLINE constexpr DoubleDouble SIN_K_PI_OVER_128[256] = { {-0x1.c57bc2e24aa15p-57, 0x1.ff621e3796d7ep-1}, {-0x1.1354d4556e4cbp-55, 0x1.ffd886084cd0dp-1}, {0, 1}, +#ifndef LIBC_MATH_HAS_SMALL_TABLES {-0x1.1354d4556e4cbp-55, 0x1.ffd886084cd0dp-1}, {-0x1.c57bc2e24aa15p-57, 0x1.ff621e3796d7ep-1}, {0x1.521ecd0c67e35p-57, 0x1.fe9cdad01883ap-1}, @@ -449,48 +338,9 @@ LIBC_INLINE constexpr DoubleDouble SIN_K_PI_OVER_128[256] = { {0x1.9a088a8bf6b2cp-59, -0x1.2d52092ce19f6p-4}, {0x1.912bd0d569a9p-61, -0x1.91f65f10dd814p-5}, {0x1.b1d63091a013p-64, -0x1.92155f7a3667ep-6}, +#endif // !LIBC_MATH_HAS_SMALL_TABLES }; -// For |x| < 2^-32, return k and u such that: -// k = round(x * 128/pi) -// x mod pi/128 = x - k * pi/128 ~ u.hi + u.lo -LIBC_INLINE unsigned range_reduction_small(double x, DoubleDouble &u) { - // Digits of pi/128, generated by Sollya with: - // > a = round(pi/128, D, RN); - // > b = round(pi/128 - a, D, RN); - constexpr DoubleDouble PI_OVER_128_DD = {0x1.1a62633145c07p-60, - 0x1.921fb54442d18p-6}; - - double prod_hi = x * ONE_TWENTY_EIGHT_OVER_PI[3][0]; - double kd = fputil::nearest_integer(prod_hi); - - // Let y = x - k * (pi/128) - // Then |y| < pi / 256 - // With extra rounding errors, we can bound |y| < 2^-6. - double y_hi = fputil::multiply_add(kd, -PI_OVER_128_DD.hi, x); // Exact - // u_hi + u_lo ~ (y_hi + kd*(-PI_OVER_128_DD[1])) - // and |u_lo| < 2* ulp(u_hi) - // The upper bound 2^-6 is over-estimated, we should still have: - // |u_hi + u_lo| < 2^-6. - u.hi = fputil::multiply_add(kd, -PI_OVER_128_DD.lo, y_hi); - u.lo = y_hi - u.hi; // Exact; - u.lo = fputil::multiply_add(kd, -PI_OVER_128_DD.lo, u.lo); - // Error bound: - // For |x| < 2^32: - // |x * high part of 128/pi| < 2^32 * 2^6 = 2^38 - // So |k| = |round(x * high part of 128/pi)| < 2^38 - // And hence, - // |(x mod pi/128) - (u.hi + u.lo)| <= ulp(2 * kd * PI_OVER_128_DD.lo) - // < 2 * 2^38 * 2^-59 * 2^-52 - // = 2^-72 - // Note: if we limit the input exponent to the same as in non-FMA version, - // i.e., |x| < 2^-23, then the output errors can be bounded by 2^-81, similar - // to the large range reduction bound. - return static_cast(static_cast(kd)); -} - -} // namespace fma - } // namespace LIBC_NAMESPACE_DECL #endif // LLVM_LIBC_SRC_MATH_GENERIC_RANGE_REDUCTION_DOUBLE_FMA_H diff --git a/libc/src/math/generic/range_reduction_double_nofma.h b/libc/src/math/generic/range_reduction_double_nofma.h index 445a45d3f9796..4ce35556b1ee8 100644 --- a/libc/src/math/generic/range_reduction_double_nofma.h +++ b/libc/src/math/generic/range_reduction_double_nofma.h @@ -15,174 +15,63 @@ #include "src/__support/FPUtil/nearest_integer.h" #include "src/__support/common.h" #include "src/__support/macros/config.h" +#include "src/__support/macros/optimization.h" +#include "src/math/generic/range_reduction_double_common.h" namespace LIBC_NAMESPACE_DECL { -namespace nofma { - using fputil::DoubleDouble; -LIBC_INLINE constexpr int FAST_PASS_EXPONENT = 23; +LIBC_INLINE unsigned LargeRangeReduction::fast(double x, DoubleDouble &u) { + using FPBits = typename fputil::FPBits; + FPBits xbits(x); -// Digits of 2^(16*i) / pi, generated by Sollya with: -// For [2..62]: -// > for i from 3 to 63 do { -// pi_inv = 2^(16*(i - 3)) / pi; -// pn = nearestint(pi_inv); -// pi_frac = pi_inv - pn; -// a = round(pi_frac, 51, RN); -// b = round(pi_frac - a, 51, RN); -// c = round(pi_frac - a - b, D, RN); -// d = round(pi_frac - a - b - c, D, RN); -// print("{", 2^7 * a, ",", 2^7 * b, ",", 2^7 * c, ",", 2^7 * d, "},"); -// }; -// For [0..1]: -// The leading bit of 2^(16*(i - 3)) / pi is very small, so we add 0.25 so that -// the conditions for the algorithms are still satisfied, and one of those -// conditions guarantees that ulp(0.25 * x_reduced) >= 2, and will safely be -// discarded. -// for i from 0 to 2 do { -// pi_frac = 0.25 + 2^(16*(i - 3)) / pi; -// a = round(pi_frac, 51, RN); -// b = round(pi_frac - a, 51, RN); -// c = round(pi_frac - a - b, D, RN); -// d = round(pi_frac - a - b - c, D, RN); -// print("{", 2^7 * a, ",", 2^7 * b, ",", 2^7 * c, ",", 2^7 * d, "},"); -// }; -// For The fast pass using double-double, we only need 3 parts (a, b, c), but -// for the accurate pass using Float128, instead of using another table of -// Float128s, we simply add the fourth path (a, b, c, d), which simplify the -// implementation a bit and saving some memory. -LIBC_INLINE constexpr double ONE_TWENTY_EIGHT_OVER_PI[64][4] = { - {0x1.0000000000014p5, 0x1.7cc1b727220a8p-49, 0x1.4fe13abe8fa9ap-101, - 0x1.bb81b6c52b328p-155}, - {0x1.0000000145f3p5, 0x1.b727220a94fep-49, 0x1.3abe8fa9a6eep-101, - 0x1.b6c52b3278872p-155}, - {0x1.000145f306dc8p5, 0x1.c882a53f84ebp-47, -0x1.70565911f924fp-101, - 0x1.2b3278872084p-155}, - {0x1.45f306dc9c884p5, -0x1.5ac07b1505c14p-47, -0x1.96447e493ad4dp-99, - 0x1.3c439041fe516p-154}, - {-0x1.f246c6efab58p4, -0x1.ec5417056591p-49, -0x1.f924eb53361dep-101, - -0x1.bef806ba71508p-156}, - {0x1.391054a7f09d4p4, 0x1.f47d4d377036cp-48, 0x1.8a5664f10e41p-100, - 0x1.fe5163abdebbcp-154}, - {0x1.529fc2757d1f4p2, 0x1.34ddc0db62958p-50, 0x1.93c439041fe51p-102, - 0x1.8eaf7aef1586ep-156}, - {-0x1.ec5417056591p-1, -0x1.f924eb53361ep-53, 0x1.c820ff28b1d5fp-105, - -0x1.443a9e48db91cp-162}, - {-0x1.505c1596447e4p5, -0x1.275a99b0ef1cp-48, 0x1.07f9458eaf7afp-100, - -0x1.d4f246dc8e2dfp-157}, - {-0x1.596447e493ad4p1, -0x1.9b0ef1bef806cp-52, 0x1.63abdebbc561bp-106, - 0x1.c91b8e909374cp-160}, - {0x1.bb81b6c52b328p5, -0x1.de37df00d74e4p-49, 0x1.5ef5de2b0db92p-101, - 0x1.b8e909374b802p-156}, - {0x1.b6c52b3278874p5, -0x1.f7c035d38a844p-47, 0x1.778ac36e48dc7p-99, - 0x1.2126e97003249p-153}, - {0x1.2b3278872084p5, -0x1.ae9c5421443a8p-50, -0x1.e48db91c5bdb2p-102, - -0x1.68ffcdb688afbp-157}, - {-0x1.8778df7c035d4p5, 0x1.d5ef5de2b0db8p-49, 0x1.2371d2126e97p-101, - 0x1.924bba8274648p-160}, - {-0x1.bef806ba71508p4, -0x1.443a9e48db91cp-50, -0x1.6f6c8b47fe6dbp-104, - -0x1.115f62e6de302p-158}, - {-0x1.ae9c5421443a8p-2, -0x1.e48db91c5bdb4p-54, 0x1.d2e006492eea1p-106, - -0x1.8b9b78c078854p-160}, - {-0x1.38a84288753c8p5, -0x1.1b7238b7b645cp-47, 0x1.c00c925dd413ap-99, - 0x1.921cfe1deb1cbp-154}, - {-0x1.0a21d4f246dc8p3, -0x1.c5bdb22d1ff9cp-50, 0x1.25dd413a3243ap-103, - -0x1.e214e34ed658cp-162}, - {-0x1.d4f246dc8e2ep3, 0x1.26e9700324978p-49, -0x1.5f62e6de301e2p-102, - -0x1.4e34ed658c117p-158}, - {-0x1.236e4716f6c8cp4, 0x1.700324977505p-49, -0x1.736f180f10a72p-101, - 0x1.62534e7dd1047p-155}, - {0x1.b8e909374b8p4, 0x1.924bba8274648p-48, 0x1.cfe1deb1cb12ap-102, - -0x1.63045df7282b4p-156}, - {0x1.09374b801924cp4, -0x1.15f62e6de302p-50, 0x1.deb1cb129a73fp-102, - -0x1.77dca0ad144bbp-158}, - {-0x1.68ffcdb688afcp3, 0x1.d1921cfe1debp-50, 0x1.cb129a73ee882p-102, - 0x1.afa975da24275p-157}, - {0x1.924bba8274648p0, 0x1.cfe1deb1cb128p-54, 0x1.a73ee88235f53p-106, - -0x1.44bb7b16638fep-162}, - {-0x1.a22bec5cdbc6p5, -0x1.e214e34ed658cp-50, -0x1.177dca0ad144cp-106, - 0x1.213a671c09ad1p-160}, - {0x1.3a32439fc3bd8p1, -0x1.c69dacb1822fp-51, 0x1.1afa975da2427p-105, - 0x1.338e04d68befdp-159}, - {-0x1.b78c0788538d4p4, 0x1.29a73ee88236p-50, -0x1.5a28976f62cc7p-103, - -0x1.fb29741037d8dp-159}, - {0x1.fc3bd63962534p5, 0x1.cfba208d7d4bcp-48, -0x1.12edec598e3f6p-100, - -0x1.4ba081bec66e3p-154}, - {-0x1.4e34ed658c118p2, 0x1.046bea5d7689p-51, 0x1.3a671c09ad17ep-104, - -0x1.bec66e29c67cbp-162}, - {0x1.62534e7dd1048p5, -0x1.415a28976f62cp-47, -0x1.8e3f652e8207p-100, - 0x1.3991d63983534p-154}, - {-0x1.63045df7282b4p4, -0x1.44bb7b16638fcp-50, -0x1.94ba081bec66ep-102, - -0x1.4e33e566305b2p-157}, - {0x1.d1046bea5d768p5, 0x1.213a671c09adp-48, 0x1.7df904e64758ep-100, - 0x1.835339f49c846p-154}, - {0x1.afa975da24274p3, 0x1.9c7026b45f7e4p-50, 0x1.3991d63983534p-106, - -0x1.82d8dee81d108p-160}, - {-0x1.a28976f62cc7p5, -0x1.fb29741037d8cp-47, -0x1.b8a719f2b3183p-100, - 0x1.3908bf177bf25p-155}, - {-0x1.76f62cc71fb28p5, -0x1.741037d8cdc54p-47, 0x1.cc1a99cfa4e42p-101, - 0x1.7e2ef7e4a0ec8p-156}, - {0x1.d338e04d68bfp5, -0x1.bec66e29c67ccp-50, 0x1.339f49c845f8cp-102, - -0x1.081b5f13801dap-156}, - {0x1.c09ad17df905p4, -0x1.9b8a719f2b318p-48, -0x1.6c6f740e8840ep-103, - 0x1.41d8ffc4bffefp-157}, - {0x1.68befc827323cp5, -0x1.38cf9598c16c8p-47, 0x1.08bf177bf2507p-99, - 0x1.8ffc4bffef02dp-153}, - {-0x1.037d8cdc538dp5, 0x1.a99cfa4e422fcp-49, 0x1.77bf250763ff1p-103, - 0x1.7ffde05980fefp-158}, - {-0x1.8cdc538cf9598p5, -0x1.82d8dee81d108p-48, -0x1.b5f13801da001p-104, - 0x1.e05980fef2f12p-158}, - {-0x1.4e33e566305bp3, -0x1.bdd03a21036cp-49, 0x1.d8ffc4bffef03p-101, - -0x1.9fc04343b9d29p-156}, - {-0x1.f2b3182d8dee8p4, -0x1.d1081b5f138p-52, -0x1.da00087e99fcp-104, - -0x1.0d0ee74a5f593p-158}, - {-0x1.8c16c6f740e88p5, -0x1.036be27003b4p-49, -0x1.0fd33f8086877p-109, - -0x1.d297d64b824b2p-164}, - {0x1.3908bf177bf24p5, 0x1.0763ff12fffbcp-47, 0x1.6603fbcbc462dp-104, - 0x1.a0a6d1f6d367fp-158}, - {0x1.7e2ef7e4a0ec8p4, -0x1.da00087e99fcp-56, -0x1.0d0ee74a5f593p-110, - 0x1.f6d367ecf27cbp-166}, - {-0x1.081b5f13801dcp4, 0x1.fff7816603fbcp-48, 0x1.788c5ad05369p-101, - -0x1.25930261b069fp-155}, - {-0x1.af89c00ed0004p5, -0x1.fa67f010d0ee8p-50, 0x1.6b414da3eda6dp-103, - -0x1.30d834f648b0cp-162}, - {-0x1.c00ed00043f4cp5, -0x1.fc04343b9d298p-48, 0x1.4da3eda6cfd9ep-103, - 0x1.3e584dba7a32p-157}, - {0x1.2fffbc0b301fcp5, 0x1.e5e2316b414dcp-47, -0x1.c125930261b07p-99, - 0x1.84dba7a31fb35p-153}, - {-0x1.0fd33f8086878p3, 0x1.8b5a0a6d1f6d4p-50, -0x1.30261b069ec91p-103, - -0x1.85ce04cb0d00bp-157}, - {-0x1.9fc04343b9d28p4, -0x1.7d64b824b2604p-48, -0x1.86c1a7b24585dp-101, - 0x1.fb34f2ff516bbp-157}, - {-0x1.0d0ee74a5f594p2, 0x1.1f6d367ecf27cp-50, 0x1.6136e9e8c7ecdp-103, - 0x1.e5fea2d7527bbp-158}, - {-0x1.dce94beb25c14p5, 0x1.a6cfd9e4f9614p-47, -0x1.22c2e70265868p-100, - -0x1.5d28ad8453814p-158}, - {-0x1.4beb25c12593p5, -0x1.30d834f648b0cp-50, 0x1.8fd9a797fa8b6p-104, - -0x1.5b08a7028341dp-159}, - {0x1.b47db4d9fb3c8p4, 0x1.f2c26dd3d18fcp-48, 0x1.9a797fa8b5d4ap-100, - -0x1.14e050683a131p-156}, - {-0x1.25930261b06ap5, 0x1.36e9e8c7ecd3cp-47, 0x1.7fa8b5d49eeb2p-100, - -0x1.41a0e84c2f8c6p-158}, - {0x1.fb3c9f2c26dd4p4, -0x1.738132c3402bcp-51, 0x1.aea4f758fd7ccp-103, - -0x1.d0985f18c10ebp-159}, - {-0x1.b069ec9161738p5, -0x1.32c3402ba515cp-51, 0x1.eeb1faf97c5edp-104, - -0x1.7c63043ad6b69p-161}, - {-0x1.ec9161738132cp5, -0x1.a015d28ad8454p-50, 0x1.faf97c5ecf41dp-104, - -0x1.821d6b5b4565p-160}, - {-0x1.61738132c3404p5, 0x1.45aea4f758fd8p-47, -0x1.a0e84c2f8c608p-102, - -0x1.d6b5b45650128p-156}, - {0x1.fb34f2ff516bcp3, -0x1.6c229c0a0d074p-49, -0x1.30be31821d6b6p-104, - 0x1.2ea6bfb5fb12p-158}, - {0x1.3cbfd45aea4f8p5, -0x1.4e050683a130cp-48, 0x1.ce7de294a4baap-104, - -0x1.404a04ee072a3p-158}, - {-0x1.5d28ad8453814p2, -0x1.a0e84c2f8c608p-54, -0x1.d6b5b45650128p-108, - -0x1.3b81ca8bdea7fp-164}, - {-0x1.15b08a702834p5, -0x1.d0985f18c10ecp-47, 0x1.4a4ba9afed7ecp-100, - 0x1.1f8d5d0856033p-154}, -}; + int x_e_m62 = xbits.get_biased_exponent() - (FPBits::EXP_BIAS + 62); + idx = static_cast((x_e_m62 >> 4) + 3); + // Scale x down by 2^(-(16 * (idx - 3)) + xbits.set_biased_exponent((x_e_m62 & 15) + FPBits::EXP_BIAS + 62); + // 2^62 <= |x_reduced| < 2^(62 + 16) = 2^78 + x_reduced = xbits.get_val(); + // x * c_hi = ph.hi + ph.lo exactly. + DoubleDouble x_split = fputil::split(x_reduced); + DoubleDouble ph = fputil::exact_mult( + x_split, x_reduced, ONE_TWENTY_EIGHT_OVER_PI[idx][0]); + // x * c_mid = pm.hi + pm.lo exactly. + DoubleDouble pm = fputil::exact_mult( + x_split, x_reduced, ONE_TWENTY_EIGHT_OVER_PI[idx][1]); + // x * c_lo = pl.hi + pl.lo exactly. + DoubleDouble pl = fputil::exact_mult( + x_split, x_reduced, ONE_TWENTY_EIGHT_OVER_PI[idx][2]); + // Extract integral parts and fractional parts of (ph.lo + pm.hi). + double sum_hi = ph.lo + pm.hi; + double kd = fputil::nearest_integer(sum_hi); + + // x * 128/pi mod 1 ~ y_hi + y_mid + y_lo + y_hi = (ph.lo - kd) + pm.hi; // Exact + y_mid = fputil::exact_add(pm.lo, pl.hi); + y_lo = pl.lo; + + // y_l = x * c_lo_2 + pl.lo + double y_l = + fputil::multiply_add(x_reduced, ONE_TWENTY_EIGHT_OVER_PI[idx][3], y_lo); + DoubleDouble y = fputil::exact_add(y_hi, y_mid.hi); + y.lo += (y_mid.lo + y_l); + + // Digits of pi/128, generated by Sollya with: + // > a = round(pi/128, D, RN); + // > b = round(pi/128 - a, D, RN); + constexpr DoubleDouble PI_OVER_128_DD = {0x1.1a62633145c07p-60, + 0x1.921fb54442d18p-6}; + + // Error bound: with {a} denote the fractional part of a, i.e.: + // {a} = a - round(a) + // Then, + // | {x * 128/pi} - (y_hi + y_lo) | <= ulp(ulp(y_hi)) <= 2^-105 + // | {x mod pi/128} - (u.hi + u.lo) | < 2 * 2^-6 * 2^-105 = 2^-110 + u = fputil::quick_mult(y, PI_OVER_128_DD); + + return static_cast(static_cast(kd)); +} // Lookup table for sin(k * pi / 128) with k = 0, ..., 255. // Table is generated with Sollya as follow: @@ -258,6 +147,7 @@ LIBC_INLINE constexpr DoubleDouble SIN_K_PI_OVER_128[256] = { {0x1.e3a843d1db55fp-53, 0x1.ff621e3796d7cp-1}, {0x1.765595d548d9ap-54, 0x1.ffd886084cd0cp-1}, {0, 1}, +#ifndef LIBC_MATH_HAS_SMALL_TABLES {0x1.765595d548d9ap-54, 0x1.ffd886084cd0cp-1}, {0x1.e3a843d1db55fp-53, 0x1.ff621e3796d7cp-1}, {-0x1.eade132f3981dp-53, 0x1.fe9cdad01883cp-1}, @@ -449,46 +339,9 @@ LIBC_INLINE constexpr DoubleDouble SIN_K_PI_OVER_128[256] = { {-0x1.ccbeeeae8129ap-56, -0x1.2d52092ce19f4p-4}, {0x1.912bd0d569a9p-61, -0x1.91f65f10dd814p-5}, {-0x1.f938a73db97fbp-58, -0x1.92155f7a3667cp-6}, +#endif // !LIBC_MATH_HAS_SMALL_TABLES }; -LIBC_INLINE unsigned range_reduction_small(double x, DoubleDouble &u) { - constexpr double ONE_TWENTY_EIGHT_OVER_PI = 0x1.45f306dc9c883p5; - - // Digits of -pi/128, generated by Sollya with: - // > a = round(-pi/128, 25, RN); - // > b = round(-pi/128 - a, 23, RN); - // > c = round(-pi/128 - a - b, 25, RN); - // > d = round(-pi/128 - a - b - c, D, RN); - // -pi/128 ~ a + b + c + d - // The precisions of the parts are chosen so that: - // 1) k * a, k * b, k * c are exact in double precision - // 2) k * b + (x - (k * a)) is exact in double precsion - constexpr double MPI_OVER_128[4] = {-0x1.921fb5p-6, -0x1.110b48p-32, - +0x1.ee59dap-56, -0x1.98a2e03707345p-83}; - - double prod_hi = x * ONE_TWENTY_EIGHT_OVER_PI; - double kd = fputil::nearest_integer(prod_hi); - - // With -pi/128 ~ a + b + c + d as in MPI_OVER_128 description: - // t = x + k * a - double t = fputil::multiply_add(kd, MPI_OVER_128[0], x); // Exact - // y_hi = t + k * b = (x + k * a) + k * b - double y_hi = fputil::multiply_add(kd, MPI_OVER_128[1], t); // Exact - // y_lo ~ k * c + k * d - double y_lo = fputil::multiply_add(kd, MPI_OVER_128[2], kd * MPI_OVER_128[3]); - // u.hi + u.lo ~ x + k * (a + b + c + d) - u = fputil::exact_add(y_hi, y_lo); - // Error bound: For |x| < 2^-23, - // |(x mod pi/128) - (u_hi + u_lo)| < ulp(y_lo) - // <= ulp(2 * x * c) - // <= ulp(2^24 * 2^-56) - // = 2^(24 - 56 - 52) - // = 2^-84 - return static_cast(static_cast(kd)); -} - -} // namespace nofma - } // namespace LIBC_NAMESPACE_DECL #endif // LLVM_LIBC_SRC_MATH_GENERIC_RANGE_REDUCTION_DOUBLE_NOFMA_H diff --git a/libc/src/math/generic/sin.cpp b/libc/src/math/generic/sin.cpp index da3d1e94b5f64..4cbcb35513efc 100644 --- a/libc/src/math/generic/sin.cpp +++ b/libc/src/math/generic/sin.cpp @@ -18,17 +18,14 @@ #include "src/__support/macros/config.h" #include "src/__support/macros/optimization.h" // LIBC_UNLIKELY #include "src/__support/macros/properties/cpu_features.h" // LIBC_TARGET_CPU_HAS_FMA +#include "src/math/generic/range_reduction_double_common.h" #include "src/math/generic/sincos_eval.h" -// TODO: We might be able to improve the performance of large range reduction of -// non-FMA targets further by operating directly on 25-bit chunks of 128/pi and -// pre-split SIN_K_PI_OVER_128, but that might double the memory footprint of -// those lookup table. -#include "range_reduction_double_common.h" - -#if ((LIBC_MATH & LIBC_MATH_SKIP_ACCURATE_PASS) != 0) -#define LIBC_MATH_SIN_SKIP_ACCURATE_PASS -#endif +#ifdef LIBC_TARGET_CPU_HAS_FMA +#include "range_reduction_double_fma.h" +#else +#include "range_reduction_double_nofma.h" +#endif // LIBC_TARGET_CPU_HAS_FMA namespace LIBC_NAMESPACE_DECL { @@ -43,33 +40,40 @@ LLVM_LIBC_FUNCTION(double, sin, (double x)) { DoubleDouble y; unsigned k; - generic::LargeRangeReduction range_reduction_large{}; + LargeRangeReduction range_reduction_large{}; - // |x| < 2^32 (with FMA) or |x| < 2^23 (w/o FMA) + // |x| < 2^16 if (LIBC_LIKELY(x_e < FPBits::EXP_BIAS + FAST_PASS_EXPONENT)) { - // |x| < 2^-26 - if (LIBC_UNLIKELY(x_e < FPBits::EXP_BIAS - 26)) { - // Signed zeros. - if (LIBC_UNLIKELY(x == 0.0)) - return x; + // |x| < 2^-7 + if (LIBC_UNLIKELY(x_e < FPBits::EXP_BIAS - 7)) { + // |x| < 2^-26 + if (LIBC_UNLIKELY(x_e < FPBits::EXP_BIAS - 26)) { + // Signed zeros. + if (LIBC_UNLIKELY(x == 0.0)) + return x; // For |x| < 2^-26, |sin(x) - x| < ulp(x)/2. #ifdef LIBC_TARGET_CPU_HAS_FMA - return fputil::multiply_add(x, -0x1.0p-54, x); + return fputil::multiply_add(x, -0x1.0p-54, x); #else - if (LIBC_UNLIKELY(x_e < 4)) { - int rounding_mode = fputil::quick_get_round(); - if (rounding_mode == FE_TOWARDZERO || - (xbits.sign() == Sign::POS && rounding_mode == FE_DOWNWARD) || - (xbits.sign() == Sign::NEG && rounding_mode == FE_UPWARD)) - return FPBits(xbits.uintval() - 1).get_val(); - } - return fputil::multiply_add(x, -0x1.0p-54, x); + if (LIBC_UNLIKELY(x_e < 4)) { + int rounding_mode = fputil::quick_get_round(); + if (rounding_mode == FE_TOWARDZERO || + (xbits.sign() == Sign::POS && rounding_mode == FE_DOWNWARD) || + (xbits.sign() == Sign::NEG && rounding_mode == FE_UPWARD)) + return FPBits(xbits.uintval() - 1).get_val(); + } + return fputil::multiply_add(x, -0x1.0p-54, x); #endif // LIBC_TARGET_CPU_HAS_FMA + } + // No range reduction needed. + k = 0; + y.lo = 0.0; + y.hi = x; + } else { + // Small range reduction. + k = range_reduction_small(x, y); } - - // // Small range reduction. - k = range_reduction_small(x, y); } else { // Inf or NaN if (LIBC_UNLIKELY(x_e > 2 * FPBits::EXP_BIAS)) { @@ -82,44 +86,33 @@ LLVM_LIBC_FUNCTION(double, sin, (double x)) { } // Large range reduction. - k = range_reduction_large.compute_high_part(x); - y = range_reduction_large.fast(); + k = range_reduction_large.fast(x, y); } DoubleDouble sin_y, cos_y; - generic::sincos_eval(y, sin_y, cos_y); + [[maybe_unused]] double err = generic::sincos_eval(y, sin_y, cos_y); // Look up sin(k * pi/128) and cos(k * pi/128) - // Memory saving versions: - - // Use 128-entry table instead: - // DoubleDouble sin_k = SIN_K_PI_OVER_128[k & 127]; - // uint64_t sin_s = static_cast(k & 128) << (63 - 7); - // sin_k.hi = FPBits(FPBits(sin_k.hi).uintval() ^ sin_s).get_val(); - // sin_k.lo = FPBits(FPBits(sin_k.hi).uintval() ^ sin_s).get_val(); - // DoubleDouble cos_k = SIN_K_PI_OVER_128[(k + 64) & 127]; - // uint64_t cos_s = static_cast((k + 64) & 128) << (63 - 7); - // cos_k.hi = FPBits(FPBits(cos_k.hi).uintval() ^ cos_s).get_val(); - // cos_k.lo = FPBits(FPBits(cos_k.hi).uintval() ^ cos_s).get_val(); - - // Use 64-entry table instead: - // auto get_idx_dd = [](unsigned kk) -> DoubleDouble { - // unsigned idx = (kk & 64) ? 64 - (kk & 63) : (kk & 63); - // DoubleDouble ans = SIN_K_PI_OVER_128[idx]; - // if (kk & 128) { - // ans.hi = -ans.hi; - // ans.lo = -ans.lo; - // } - // return ans; - // }; - // DoubleDouble sin_k = get_idx_dd(k); - // DoubleDouble cos_k = get_idx_dd(k + 64); - +#ifdef LIBC_MATH_HAS_SMALL_TABLES + // Memory saving versions. Use 65-entry table. + auto get_idx_dd = [](unsigned kk) -> DoubleDouble { + unsigned idx = (kk & 64) ? 64 - (kk & 63) : (kk & 63); + DoubleDouble ans = SIN_K_PI_OVER_128[idx]; + if (kk & 128) { + ans.hi = -ans.hi; + ans.lo = -ans.lo; + } + return ans; + }; + DoubleDouble sin_k = get_idx_dd(k); + DoubleDouble cos_k = get_idx_dd(k + 64); +#else // Fast look up version, but needs 256-entry table. // cos(k * pi/128) = sin(k * pi/128 + pi/2) = sin((k + 64) * pi/128). DoubleDouble sin_k = SIN_K_PI_OVER_128[k & 255]; DoubleDouble cos_k = SIN_K_PI_OVER_128[(k + 64) & 255]; +#endif // After range reduction, k = round(x * 128 / pi) and y = x - k * (pi / 128). // So k is an integer and -pi / 256 <= y <= pi / 256. @@ -131,20 +124,13 @@ LLVM_LIBC_FUNCTION(double, sin, (double x)) { DoubleDouble rr = fputil::exact_add(sin_k_cos_y.hi, cos_k_sin_y.hi); rr.lo += sin_k_cos_y.lo + cos_k_sin_y.lo; -#ifdef LIBC_MATH_SIN_SKIP_ACCURATE_PASS +#ifdef LIBC_MATH_HAS_SKIP_ACCURATE_PASS return rr.hi + rr.lo; #else // Accurate test and pass for correctly rounded implementation. -#ifdef LIBC_TARGET_CPU_HAS_FMA - constexpr double ERR = 0x1.0p-70; -#else - // TODO: Improve non-FMA fast pass accuracy. - constexpr double ERR = 0x1.0p-66; -#endif // LIBC_TARGET_CPU_HAS_FMA - - double rlp = rr.lo + ERR; - double rlm = rr.lo - ERR; + double rlp = rr.lo + err; + double rlm = rr.lo - err; double r_upper = rr.hi + rlp; // (rr.lo + ERR); double r_lower = rr.hi + rlm; // (rr.lo - ERR); @@ -155,7 +141,7 @@ LLVM_LIBC_FUNCTION(double, sin, (double x)) { Float128 u_f128, sin_u, cos_u; if (LIBC_LIKELY(x_e < FPBits::EXP_BIAS + FAST_PASS_EXPONENT)) - u_f128 = generic::range_reduction_small_f128(x); + u_f128 = range_reduction_small_f128(x); else u_f128 = range_reduction_large.accurate(); @@ -163,7 +149,7 @@ LLVM_LIBC_FUNCTION(double, sin, (double x)) { auto get_sin_k = [](unsigned kk) -> Float128 { unsigned idx = (kk & 64) ? 64 - (kk & 63) : (kk & 63); - Float128 ans = generic::SIN_K_PI_OVER_128_F128[idx]; + Float128 ans = SIN_K_PI_OVER_128_F128[idx]; if (kk & 128) ans.sign = Sign::NEG; return ans; @@ -182,7 +168,7 @@ LLVM_LIBC_FUNCTION(double, sin, (double x)) { // https://github.com/llvm/llvm-project/issues/96452. return static_cast(r); -#endif // !LIBC_MATH_SIN_SKIP_ACCURATE_PASS +#endif // !LIBC_MATH_HAS_SKIP_ACCURATE_PASS } } // namespace LIBC_NAMESPACE_DECL diff --git a/libc/src/math/generic/sincos.cpp b/libc/src/math/generic/sincos.cpp index 1af0ee7b0eb2c..314504b3350b5 100644 --- a/libc/src/math/generic/sincos.cpp +++ b/libc/src/math/generic/sincos.cpp @@ -19,17 +19,14 @@ #include "src/__support/macros/config.h" #include "src/__support/macros/optimization.h" // LIBC_UNLIKELY #include "src/__support/macros/properties/cpu_features.h" // LIBC_TARGET_CPU_HAS_FMA +#include "src/math/generic/range_reduction_double_common.h" #include "src/math/generic/sincos_eval.h" -// TODO: We might be able to improve the performance of large range reduction of -// non-FMA targets further by operating directly on 25-bit chunks of 128/pi and -// pre-split SIN_K_PI_OVER_128, but that might double the memory footprint of -// those lookup table. -#include "range_reduction_double_common.h" - -#if ((LIBC_MATH & LIBC_MATH_SKIP_ACCURATE_PASS) != 0) -#define LIBC_MATH_SINCOS_SKIP_ACCURATE_PASS -#endif +#ifdef LIBC_TARGET_CPU_HAS_FMA +#include "range_reduction_double_fma.h" +#else +#include "range_reduction_double_nofma.h" +#endif // LIBC_TARGET_CPU_HAS_FMA namespace LIBC_NAMESPACE_DECL { @@ -44,40 +41,47 @@ LLVM_LIBC_FUNCTION(void, sincos, (double x, double *sin_x, double *cos_x)) { DoubleDouble y; unsigned k; - generic::LargeRangeReduction range_reduction_large{}; + LargeRangeReduction range_reduction_large{}; - // |x| < 2^32 (with FMA) or |x| < 2^23 (w/o FMA) + // |x| < 2^16 if (LIBC_LIKELY(x_e < FPBits::EXP_BIAS + FAST_PASS_EXPONENT)) { - // |x| < 2^-27 - if (LIBC_UNLIKELY(x_e < FPBits::EXP_BIAS - 27)) { - // Signed zeros. - if (LIBC_UNLIKELY(x == 0.0)) { - *sin_x = x; - *cos_x = 1.0; - return; - } - - // For |x| < 2^-27, max(|sin(x) - x|, |cos(x) - 1|) < ulp(x)/2. + // |x| < 2^-7 + if (LIBC_UNLIKELY(x_e < FPBits::EXP_BIAS - 7)) { + // |x| < 2^-27 + if (LIBC_UNLIKELY(x_e < FPBits::EXP_BIAS - 27)) { + // Signed zeros. + if (LIBC_UNLIKELY(x == 0.0)) { + *sin_x = x; + *cos_x = 1.0; + return; + } + + // For |x| < 2^-27, max(|sin(x) - x|, |cos(x) - 1|) < ulp(x)/2. #ifdef LIBC_TARGET_CPU_HAS_FMA - *sin_x = fputil::multiply_add(x, -0x1.0p-54, x); - *cos_x = fputil::multiply_add(x, -x, 1.0); + *sin_x = fputil::multiply_add(x, -0x1.0p-54, x); + *cos_x = fputil::multiply_add(x, -x, 1.0); #else - *cos_x = fputil::round_result_slightly_down(1.0); - - if (LIBC_UNLIKELY(x_e < 4)) { - int rounding_mode = fputil::quick_get_round(); - if (rounding_mode == FE_TOWARDZERO || - (xbits.sign() == Sign::POS && rounding_mode == FE_DOWNWARD) || - (xbits.sign() == Sign::NEG && rounding_mode == FE_UPWARD)) - *sin_x = FPBits(xbits.uintval() - 1).get_val(); - } - *sin_x = fputil::multiply_add(x, -0x1.0p-54, x); + *cos_x = fputil::round_result_slightly_down(1.0); + + if (LIBC_UNLIKELY(x_e < 4)) { + int rounding_mode = fputil::quick_get_round(); + if (rounding_mode == FE_TOWARDZERO || + (xbits.sign() == Sign::POS && rounding_mode == FE_DOWNWARD) || + (xbits.sign() == Sign::NEG && rounding_mode == FE_UPWARD)) + *sin_x = FPBits(xbits.uintval() - 1).get_val(); + } + *sin_x = fputil::multiply_add(x, -0x1.0p-54, x); #endif // LIBC_TARGET_CPU_HAS_FMA - return; + return; + } + // No range reduction needed. + k = 0; + y.lo = 0.0; + y.hi = x; + } else { + // Small range reduction. + k = range_reduction_small(x, y); } - - // // Small range reduction. - k = range_reduction_small(x, y); } else { // Inf or NaN if (LIBC_UNLIKELY(x_e > 2 * FPBits::EXP_BIAS)) { @@ -91,44 +95,34 @@ LLVM_LIBC_FUNCTION(void, sincos, (double x, double *sin_x, double *cos_x)) { } // Large range reduction. - k = range_reduction_large.compute_high_part(x); - y = range_reduction_large.fast(); + k = range_reduction_large.fast(x, y); } DoubleDouble sin_y, cos_y; - generic::sincos_eval(y, sin_y, cos_y); + [[maybe_unused]] double err = generic::sincos_eval(y, sin_y, cos_y); // Look up sin(k * pi/128) and cos(k * pi/128) - // Memory saving versions: - - // Use 128-entry table instead: - // DoubleDouble sin_k = SIN_K_PI_OVER_128[k & 127]; - // uint64_t sin_s = static_cast(k & 128) << (63 - 7); - // sin_k.hi = FPBits(FPBits(sin_k.hi).uintval() ^ sin_s).get_val(); - // sin_k.lo = FPBits(FPBits(sin_k.hi).uintval() ^ sin_s).get_val(); - // DoubleDouble cos_k = SIN_K_PI_OVER_128[(k + 64) & 127]; - // uint64_t cos_s = static_cast((k + 64) & 128) << (63 - 7); - // cos_k.hi = FPBits(FPBits(cos_k.hi).uintval() ^ cos_s).get_val(); - // cos_k.lo = FPBits(FPBits(cos_k.hi).uintval() ^ cos_s).get_val(); - - // Use 64-entry table instead: - // auto get_idx_dd = [](unsigned kk) -> DoubleDouble { - // unsigned idx = (kk & 64) ? 64 - (kk & 63) : (kk & 63); - // DoubleDouble ans = SIN_K_PI_OVER_128[idx]; - // if (kk & 128) { - // ans.hi = -ans.hi; - // ans.lo = -ans.lo; - // } - // return ans; - // }; - // DoubleDouble sin_k = get_idx_dd(k); - // DoubleDouble cos_k = get_idx_dd(k + 64); - +#ifdef LIBC_MATH_HAS_SMALL_TABLES + // Memory saving versions. Use 65-entry table. + auto get_idx_dd = [](unsigned kk) -> DoubleDouble { + unsigned idx = (kk & 64) ? 64 - (kk & 63) : (kk & 63); + DoubleDouble ans = SIN_K_PI_OVER_128[idx]; + if (kk & 128) { + ans.hi = -ans.hi; + ans.lo = -ans.lo; + } + return ans; + }; + DoubleDouble sin_k = get_idx_dd(k); + DoubleDouble cos_k = get_idx_dd(k + 64); +#else // Fast look up version, but needs 256-entry table. // cos(k * pi/128) = sin(k * pi/128 + pi/2) = sin((k + 64) * pi/128). DoubleDouble sin_k = SIN_K_PI_OVER_128[k & 255]; DoubleDouble cos_k = SIN_K_PI_OVER_128[(k + 64) & 255]; +#endif // LIBC_MATH_HAS_SMALL_TABLES + DoubleDouble msin_k{-sin_k.lo, -sin_k.hi}; // After range reduction, k = round(x * 128 / pi) and y = x - k * (pi / 128). @@ -149,24 +143,17 @@ LLVM_LIBC_FUNCTION(void, sincos, (double x, double *sin_x, double *cos_x)) { sin_dd.lo += sin_k_cos_y.lo + cos_k_sin_y.lo; cos_dd.lo += msin_k_sin_y.lo + cos_k_cos_y.lo; -#ifdef LIBC_MATH_SINCOS_SKIP_ACCURATE_PASS +#ifdef LIBC_MATH_HAS_SKIP_ACCURATE_PASS *sin_x = sin_dd.hi + sin_dd.lo; *cos_x = cos_dd.hi + cos_dd.lo; return; #else // Accurate test and pass for correctly rounded implementation. -#ifdef LIBC_TARGET_CPU_HAS_FMA - constexpr double ERR = 0x1.0p-70; -#else - // TODO: Improve non-FMA fast pass accuracy. - constexpr double ERR = 0x1.0p-66; -#endif // LIBC_TARGET_CPU_HAS_FMA - - double sin_lp = sin_dd.lo + ERR; - double sin_lm = sin_dd.lo - ERR; - double cos_lp = cos_dd.lo + ERR; - double cos_lm = cos_dd.lo - ERR; + double sin_lp = sin_dd.lo + err; + double sin_lm = sin_dd.lo - err; + double cos_lp = cos_dd.lo + err; + double cos_lm = cos_dd.lo - err; double sin_upper = sin_dd.hi + sin_lp; double sin_lower = sin_dd.hi + sin_lm; @@ -182,7 +169,7 @@ LLVM_LIBC_FUNCTION(void, sincos, (double x, double *sin_x, double *cos_x)) { Float128 u_f128, sin_u, cos_u; if (LIBC_LIKELY(x_e < FPBits::EXP_BIAS + FAST_PASS_EXPONENT)) - u_f128 = generic::range_reduction_small_f128(x); + u_f128 = range_reduction_small_f128(x); else u_f128 = range_reduction_large.accurate(); @@ -190,7 +177,7 @@ LLVM_LIBC_FUNCTION(void, sincos, (double x, double *sin_x, double *cos_x)) { auto get_sin_k = [](unsigned kk) -> Float128 { unsigned idx = (kk & 64) ? 64 - (kk & 63) : (kk & 63); - Float128 ans = generic::SIN_K_PI_OVER_128_F128[idx]; + Float128 ans = SIN_K_PI_OVER_128_F128[idx]; if (kk & 128) ans.sign = Sign::NEG; return ans; @@ -222,7 +209,7 @@ LLVM_LIBC_FUNCTION(void, sincos, (double x, double *sin_x, double *cos_x)) { fputil::quick_add(fputil::quick_mul(cos_k_f128, cos_u), fputil::quick_mul(msin_k_f128, sin_u))); -#endif // !LIBC_MATH_SINCOS_SKIP_ACCURATE_PASS +#endif // !LIBC_MATH_HAS_SKIP_ACCURATE_PASS } } // namespace LIBC_NAMESPACE_DECL diff --git a/libc/src/math/generic/sincos_eval.h b/libc/src/math/generic/sincos_eval.h index e491467c5663f..6cd1da4721bf5 100644 --- a/libc/src/math/generic/sincos_eval.h +++ b/libc/src/math/generic/sincos_eval.h @@ -23,8 +23,8 @@ namespace generic { using fputil::DoubleDouble; using Float128 = fputil::DyadicFloat<128>; -LIBC_INLINE void sincos_eval(const DoubleDouble &u, DoubleDouble &sin_u, - DoubleDouble &cos_u) { +LIBC_INLINE double sincos_eval(const DoubleDouble &u, DoubleDouble &sin_u, + DoubleDouble &cos_u) { // Evaluate sin(y) = sin(x - k * (pi/128)) // We use the degree-7 Taylor approximation: // sin(y) ~ y - y^3/3! + y^5/5! - y^7/7! @@ -61,9 +61,19 @@ LIBC_INLINE void sincos_eval(const DoubleDouble &u, DoubleDouble &sin_u, // + u_hi u_lo (-1 + u_hi^2/6) // We compute 1 - u_hi^2 accurately: // v_hi + v_lo ~ 1 - u_hi^2/2 - double v_hi = fputil::multiply_add(u.hi, u.hi * (-0.5), 1.0); - double v_lo = 1.0 - v_hi; // Exact - v_lo = fputil::multiply_add(u.hi, u.hi * (-0.5), v_lo); + // with error <= 2^-105. + double u_hi_neg_half = (-0.5) * u.hi; + DoubleDouble v; + +#ifdef LIBC_TARGET_CPU_HAS_FMA + v.hi = fputil::multiply_add(u.hi, u_hi_neg_half, 1.0); + v.lo = 1.0 - v.hi; // Exact + v.lo = fputil::multiply_add(u.hi, u_hi_neg_half, v.lo); +#else + DoubleDouble u_hi_sq_neg_half = fputil::exact_mult(u.hi, u_hi_neg_half); + v = fputil::exact_add(1.0, u_hi_sq_neg_half.hi); + v.lo += u_hi_sq_neg_half.lo; +#endif // LIBC_TARGET_CPU_HAS_FMA // r1 ~ -1/720 + u_hi^2 / 40320 double r1 = fputil::multiply_add(u_hi_sq, 0x1.a01a01a01a01ap-16, @@ -75,12 +85,15 @@ LIBC_INLINE void sincos_eval(const DoubleDouble &u, DoubleDouble &sin_u, // r2 ~ 1/24 + u_hi^2 (-1/720 + u_hi^2 / 40320) double r2 = fputil::multiply_add(u_hi_sq, r1, 0x1.5555555555555p-5); // s2 ~ v_lo + u_hi * u_lo * (-1 + u_hi^2 / 6) - double s2 = fputil::multiply_add(u_hi_u_lo, s1, v_lo); + double s2 = fputil::multiply_add(u_hi_u_lo, s1, v.lo); double cos_lo = fputil::multiply_add(u_hi_4, r2, s2); // Overall, |cos(y) - (v_hi + cos_lo)| < 2*ulp(u_hi^4) < 2^-75. sin_u = fputil::exact_add(u.hi, sin_lo); - cos_u = fputil::exact_add(v_hi, cos_lo); + cos_u = fputil::exact_add(v.hi, cos_lo); + + return fputil::multiply_add(fputil::FPBits(u_hi_3).abs().get_val(), + 0x1.0p-51, 0x1.0p-105); } LIBC_INLINE void sincos_eval(const Float128 &u, Float128 &sin_u, diff --git a/libc/src/math/generic/tan.cpp b/libc/src/math/generic/tan.cpp index 45fd8bb9156be..8c19936068dda 100644 --- a/libc/src/math/generic/tan.cpp +++ b/libc/src/math/generic/tan.cpp @@ -20,16 +20,13 @@ #include "src/__support/macros/config.h" #include "src/__support/macros/optimization.h" // LIBC_UNLIKELY #include "src/__support/macros/properties/cpu_features.h" // LIBC_TARGET_CPU_HAS_FMA +#include "src/math/generic/range_reduction_double_common.h" -// TODO: We might be able to improve the performance of large range reduction of -// non-FMA targets further by operating directly on 25-bit chunks of 128/pi and -// pre-split SIN_K_PI_OVER_128, but that might double the memory footprint of -// those lookup table. -#include "range_reduction_double_common.h" - -#if ((LIBC_MATH & LIBC_MATH_SKIP_ACCURATE_PASS) != 0) -#define LIBC_MATH_TAN_SKIP_ACCURATE_PASS -#endif +#ifdef LIBC_TARGET_CPU_HAS_FMA +#include "range_reduction_double_fma.h" +#else +#include "range_reduction_double_nofma.h" +#endif // LIBC_TARGET_CPU_HAS_FMA namespace LIBC_NAMESPACE_DECL { @@ -38,7 +35,7 @@ using Float128 = typename fputil::DyadicFloat<128>; namespace { -LIBC_INLINE DoubleDouble tan_eval(const DoubleDouble &u) { +LIBC_INLINE double tan_eval(const DoubleDouble &u, DoubleDouble &result) { // Evaluate tan(y) = tan(x - k * (pi/128)) // We use the degree-9 Taylor approximation: // tan(y) ~ P(y) = y + y^3/3 + 2*y^5/15 + 17*y^7/315 + 62*y^9/2835 @@ -69,10 +66,12 @@ LIBC_INLINE DoubleDouble tan_eval(const DoubleDouble &u) { // Overall, |tan(y) - (u_hi + tan_lo)| < ulp(u_hi^3) <= 2^-71. // And the relative errors is: // |(tan(y) - (u_hi + tan_lo)) / tan(y) | <= 2*ulp(u_hi^2) < 2^-64 - - return fputil::exact_add(u.hi, tan_lo); + result = fputil::exact_add(u.hi, tan_lo); + return fputil::multiply_add(fputil::FPBits(u_hi_3).abs().get_val(), + 0x1.0p-51, 0x1.0p-102); } +#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS // Accurate evaluation of tan for small u. [[maybe_unused]] Float128 tan_eval(const Float128 &u) { Float128 u_sq = fputil::quick_mul(u, u); @@ -117,6 +116,7 @@ LIBC_INLINE DoubleDouble tan_eval(const DoubleDouble &u) { fputil::quick_mul(q1, fputil::quick_add(TWO, fputil::quick_mul(b, q1))); return fputil::quick_mul(a, q2); } +#endif // !LIBC_MATH_HAS_SKIP_ACCURATE_PASS } // anonymous namespace @@ -128,33 +128,39 @@ LLVM_LIBC_FUNCTION(double, tan, (double x)) { DoubleDouble y; unsigned k; - generic::LargeRangeReduction range_reduction_large{}; + LargeRangeReduction range_reduction_large{}; - // |x| < 2^32 (with FMA) or |x| < 2^23 (w/o FMA) + // |x| < 2^16 if (LIBC_LIKELY(x_e < FPBits::EXP_BIAS + FAST_PASS_EXPONENT)) { - // |x| < 2^-27 - if (LIBC_UNLIKELY(x_e < FPBits::EXP_BIAS - 27)) { - // Signed zeros. - if (LIBC_UNLIKELY(x == 0.0)) - return x; + // |x| < 2^-7 + if (LIBC_UNLIKELY(x_e < FPBits::EXP_BIAS - 7)) { + // |x| < 2^-27 + if (LIBC_UNLIKELY(x_e < FPBits::EXP_BIAS - 27)) { + // Signed zeros. + if (LIBC_UNLIKELY(x == 0.0)) + return x; // For |x| < 2^-27, |tan(x) - x| < ulp(x)/2. #ifdef LIBC_TARGET_CPU_HAS_FMA - return fputil::multiply_add(x, 0x1.0p-54, x); + return fputil::multiply_add(x, 0x1.0p-54, x); #else - if (LIBC_UNLIKELY(x_e < 4)) { - int rounding_mode = fputil::quick_get_round(); - if (rounding_mode == FE_TOWARDZERO || - (xbits.sign() == Sign::POS && rounding_mode == FE_DOWNWARD) || - (xbits.sign() == Sign::NEG && rounding_mode == FE_UPWARD)) - return FPBits(xbits.uintval() + 1).get_val(); - } - return fputil::multiply_add(x, 0x1.0p-54, x); + if (LIBC_UNLIKELY(x_e < 4)) { + int rounding_mode = fputil::quick_get_round(); + if ((xbits.sign() == Sign::POS && rounding_mode == FE_UPWARD) || + (xbits.sign() == Sign::NEG && rounding_mode == FE_DOWNWARD)) + return FPBits(xbits.uintval() + 1).get_val(); + } + return fputil::multiply_add(x, 0x1.0p-54, x); #endif // LIBC_TARGET_CPU_HAS_FMA + } + // No range reduction needed. + k = 0; + y.lo = 0.0; + y.hi = x; + } else { + // Small range reduction. + k = range_reduction_small(x, y); } - - // // Small range reduction. - k = range_reduction_small(x, y); } else { // Inf or NaN if (LIBC_UNLIKELY(x_e > 2 * FPBits::EXP_BIAS)) { @@ -167,42 +173,32 @@ LLVM_LIBC_FUNCTION(double, tan, (double x)) { } // Large range reduction. - k = range_reduction_large.compute_high_part(x); - y = range_reduction_large.fast(); + k = range_reduction_large.fast(x, y); } - DoubleDouble tan_y = tan_eval(y); + DoubleDouble tan_y; + [[maybe_unused]] double err = tan_eval(y, tan_y); // Look up sin(k * pi/128) and cos(k * pi/128) - // Memory saving versions: - - // Use 128-entry table instead: - // DoubleDouble sin_k = SIN_K_PI_OVER_128[k & 127]; - // uint64_t sin_s = static_cast(k & 128) << (63 - 7); - // sin_k.hi = FPBits(FPBits(sin_k.hi).uintval() ^ sin_s).get_val(); - // sin_k.lo = FPBits(FPBits(sin_k.hi).uintval() ^ sin_s).get_val(); - // DoubleDouble cos_k = SIN_K_PI_OVER_128[(k + 64) & 127]; - // uint64_t cos_s = static_cast((k + 64) & 128) << (63 - 7); - // cos_k.hi = FPBits(FPBits(cos_k.hi).uintval() ^ cos_s).get_val(); - // cos_k.lo = FPBits(FPBits(cos_k.hi).uintval() ^ cos_s).get_val(); - - // Use 64-entry table instead: - // auto get_idx_dd = [](unsigned kk) -> DoubleDouble { - // unsigned idx = (kk & 64) ? 64 - (kk & 63) : (kk & 63); - // DoubleDouble ans = SIN_K_PI_OVER_128[idx]; - // if (kk & 128) { - // ans.hi = -ans.hi; - // ans.lo = -ans.lo; - // } - // return ans; - // }; - // DoubleDouble msin_k = get_idx_dd(k + 128); - // DoubleDouble cos_k = get_idx_dd(k + 64); - +#ifdef LIBC_MATH_HAS_SMALL_TABLES + // Memory saving versions. Use 65-entry table: + auto get_idx_dd = [](unsigned kk) -> DoubleDouble { + unsigned idx = (kk & 64) ? 64 - (kk & 63) : (kk & 63); + DoubleDouble ans = SIN_K_PI_OVER_128[idx]; + if (kk & 128) { + ans.hi = -ans.hi; + ans.lo = -ans.lo; + } + return ans; + }; + DoubleDouble msin_k = get_idx_dd(k + 128); + DoubleDouble cos_k = get_idx_dd(k + 64); +#else // Fast look up version, but needs 256-entry table. // cos(k * pi/128) = sin(k * pi/128 + pi/2) = sin((k + 64) * pi/128). DoubleDouble msin_k = SIN_K_PI_OVER_128[(k + 128) & 255]; DoubleDouble cos_k = SIN_K_PI_OVER_128[(k + 64) & 255]; +#endif // LIBC_MATH_HAS_SMALL_TABLES // After range reduction, k = round(x * 128 / pi) and y = x - k * (pi / 128). // So k is an integer and -pi / 256 <= y <= pi / 256. @@ -222,7 +218,7 @@ LLVM_LIBC_FUNCTION(double, tan, (double x)) { num_dd.lo += cos_k_tan_y.lo - msin_k.lo; den_dd.lo += msin_k_tan_y.lo + cos_k.lo; -#ifdef LIBC_MATH_TAN_SKIP_ACCURATE_PASS +#ifdef LIBC_MATH_HAS_SKIP_ACCURATE_PASS double tan_x = (num_dd.hi + num_dd.lo) / (den_dd.hi + den_dd.lo); return tan_x; #else @@ -231,18 +227,16 @@ LLVM_LIBC_FUNCTION(double, tan, (double x)) { // Accurate double-double division DoubleDouble tan_x = fputil::div(num_dd, den_dd); - // Relative errors for k != 0 mod 64 is: - // absolute errors / min(sin(k*pi/128), cos(k*pi/128)) <= 2^-71 / 2^-7 - // = 2^-64. - // For k = 0 mod 64, the relative errors is bounded by: - // 2^-71 / 2^(exponent of x). - constexpr int ERR = 64; + // Simple error bound: |1 / den_dd| < 2^(1 + floor(-log2(den_dd)))). + uint64_t den_inv = (static_cast(FPBits::EXP_BIAS + 1) + << (FPBits::FRACTION_LEN + 1)) - + (FPBits(den_dd.hi).uintval() & FPBits::EXP_MASK); - int y_exp = 7 + FPBits(y.hi).get_exponent(); - int rel_err_exp = ERR + static_cast((k & 63) == 0) * y_exp; - int64_t tan_x_err = static_cast(FPBits(tan_x.hi).uintval()) - - (static_cast(rel_err_exp) << 52); - double tan_err = FPBits(static_cast(tan_x_err)).get_val(); + // For tan_x = (num_dd + err) / (den_dd + err), the error is bounded by: + // | tan_x - num_dd / den_dd | <= err * ( 1 + | tan_x * den_dd | ). + double tan_err = + err * fputil::multiply_add(FPBits(den_inv).get_val(), + FPBits(tan_x.hi).abs().get_val(), 1.0); double err_higher = tan_x.lo + tan_err; double err_lower = tan_x.lo - tan_err; @@ -256,7 +250,7 @@ LLVM_LIBC_FUNCTION(double, tan, (double x)) { Float128 u_f128; if (LIBC_LIKELY(x_e < FPBits::EXP_BIAS + FAST_PASS_EXPONENT)) - u_f128 = generic::range_reduction_small_f128(x); + u_f128 = range_reduction_small_f128(x); else u_f128 = range_reduction_large.accurate(); @@ -264,7 +258,7 @@ LLVM_LIBC_FUNCTION(double, tan, (double x)) { auto get_sin_k = [](unsigned kk) -> Float128 { unsigned idx = (kk & 64) ? 64 - (kk & 63) : (kk & 63); - Float128 ans = generic::SIN_K_PI_OVER_128_F128[idx]; + Float128 ans = SIN_K_PI_OVER_128_F128[idx]; if (kk & 128) ans.sign = Sign::NEG; return ans; @@ -292,7 +286,7 @@ LLVM_LIBC_FUNCTION(double, tan, (double x)) { // https://github.com/llvm/llvm-project/issues/96452. return static_cast(result); -#endif // !LIBC_MATH_TAN_SKIP_ACCURATE_PASS +#endif // !LIBC_MATH_HAS_SKIP_ACCURATE_PASS } } // namespace LIBC_NAMESPACE_DECL diff --git a/libc/test/src/math/cos_test.cpp b/libc/test/src/math/cos_test.cpp index 484d47fd3e96c..e2d47917e545e 100644 --- a/libc/test/src/math/cos_test.cpp +++ b/libc/test/src/math/cos_test.cpp @@ -50,8 +50,7 @@ TEST_F(LlvmLibcCosTest, TrickyInputs) { 0x1.2b5fe88a9d8d5p+903, 0x1.f6d7518808571p+1023, -0x1.a880417b7b119p+1023, 0x1.00a33764a0a83p-7, 0x1.fe81868fc47fep+1, 0x1.0da8cc189b47dp-10, - 0x1.da1838053b866p+5, - + 0x1.da1838053b866p+5, 0x1.ffffffffe854bp199, }; constexpr int N = sizeof(INPUTS) / sizeof(INPUTS[0]); diff --git a/libc/test/src/math/sin_test.cpp b/libc/test/src/math/sin_test.cpp index 60f6ef5c84463..d4c6bd416a409 100644 --- a/libc/test/src/math/sin_test.cpp +++ b/libc/test/src/math/sin_test.cpp @@ -20,11 +20,13 @@ using LIBC_NAMESPACE::testing::tlog; TEST_F(LlvmLibcSinTest, TrickyInputs) { constexpr double INPUTS[] = { - 0x1.940c877fb7dacp-7, 0x1.fffffffffdb6p24, 0x1.fd4da4ef37075p29, - 0x1.b951f1572eba5p+31, 0x1.55202aefde314p+31, 0x1.85fc0f04c0128p101, - 0x1.7776c2343ba4ep101, 0x1.678309fa50d58p110, 0x1.fffffffffef4ep199, - -0x1.ab514bfc61c76p+7, -0x1.f7898d5a756ddp+2, -0x1.f42fb19b5b9b2p-6, - 0x1.5f09cad750ab1p+3, -0x1.14823229799c2p+7, -0x1.0285070f9f1bcp-5, + 0x1.5f09cad750ab1p+3, 0x1.fff781921b61fp15, -0x1.f635b70b92407p-21, + -0x1.3ecf146c39c0cp-20, 0x1.6ac5b262ca1ffp849, 0x1.6c6cbc45dc8dep5, + 0x1.921fb5443p-7, 0x1.940c877fb7dacp-7, 0x1.fffffffffdb6p24, + 0x1.fd4da4ef37075p29, 0x1.b951f1572eba5p+31, 0x1.55202aefde314p+31, + 0x1.85fc0f04c0128p101, 0x1.7776c2343ba4ep101, 0x1.678309fa50d58p110, + 0x1.fffffffffef4ep199, -0x1.ab514bfc61c76p+7, -0x1.f7898d5a756ddp+2, + -0x1.f42fb19b5b9b2p-6, -0x1.14823229799c2p+7, -0x1.0285070f9f1bcp-5, 0x1.23f40dccdef72p+0, 0x1.43cf16358c9d7p+0, 0x1.addf3b9722265p+0, 0x1.48ff1782ca91dp+8, 0x1.a211877de55dbp+4, 0x1.dcbfda0c7559ep+8, 0x1.1ffb509f3db15p+5, 0x1.2345d1e090529p+5, 0x1.ae945054939c2p+10, diff --git a/libc/test/src/math/tan_test.cpp b/libc/test/src/math/tan_test.cpp index 1ca67afdaddf2..12dfc02bac111 100644 --- a/libc/test/src/math/tan_test.cpp +++ b/libc/test/src/math/tan_test.cpp @@ -20,17 +20,20 @@ using LIBC_NAMESPACE::testing::tlog; TEST_F(LlvmLibcTanTest, TrickyInputs) { constexpr double INPUTS[] = { - 0x1.d130383d17321p-27, 0x1.8000000000009p-23, 0x1.8000000000024p-22, - 0x1.800000000009p-21, 0x1.20000000000f3p-20, 0x1.800000000024p-20, - 0x1.e0000000001c2p-20, 0x1.00452f0e0134dp-13, 0x1.0da8cc189b47dp-10, - 0x1.00a33764a0a83p-7, 0x1.911a18779813fp-7, 0x1.940c877fb7dacp-7, - 0x1.f42fb19b5b9b2p-6, 0x1.0285070f9f1bcp-5, 0x1.89f0f5241255bp-2, + 0x0.0000000000001p-1022, 0x1.d130383d17321p-27, 0x1.8000000000009p-23, + 0x1.8000000000024p-22, 0x1.800000000009p-21, 0x1.20000000000f3p-20, + 0x1.800000000024p-20, 0x1.e0000000001c2p-20, 0x1.00452f0e0134dp-13, + 0x1.0da8cc189b47dp-10, 0x1.00a33764a0a83p-7, 0x1.911a18779813fp-7, + 0x1.940c877fb7dacp-7, 0x1.f42fb19b5b9b2p-6, 0x1.0285070f9f1bcp-5, + 0x1.90e833c6969c7p-4, 0x1.91d4b77c527eap-3, 0x1.89f0f5241255bp-2, 0x1.6ca9ef729af76p-1, 0x1.23f40dccdef72p+0, 0x1.43cf16358c9d7p+0, + 0x1.90f422b49115ep+0, 0x1.9220efee9fc7ep+0, 0x1.a224411cdebcep+0, 0x1.addf3b9722265p+0, 0x1.ae78d360afa15p+0, 0x1.fe81868fc47fep+1, - 0x1.e31b55306f22cp+2, 0x1.e639103a05997p+2, 0x1.f7898d5a756ddp+2, - 0x1.1685973506319p+3, 0x1.5f09cad750ab1p+3, 0x1.aaf85537ea4c7p+3, - 0x1.4f2b874135d27p+4, 0x1.13114266f9764p+4, 0x1.a211877de55dbp+4, - 0x1.a5eece87e8606p+4, 0x1.a65d441ea6dcep+4, 0x1.045457ae3994p+5, + 0x1.e31b55306f22cp+2, 0x1.e639103a05997p+2, 0x1.f69d074a3358fp+2, + 0x1.f7898d5a756ddp+2, 0x1.1685973506319p+3, 0x1.5f09cad750ab1p+3, + 0x1.aaf85537ea4c7p+3, 0x1.c50ddc4f513b4p+3, 0x1.13114266f9764p+4, + 0x1.4f2b874135d27p+4, 0x1.a211877de55dbp+4, 0x1.a5eece87e8606p+4, + 0x1.a65d441ea6dcep+4, 0x1.ab8c2f8ab5b7p+4, 0x1.045457ae3994p+5, 0x1.1ffb509f3db15p+5, 0x1.2345d1e090529p+5, 0x1.c96e28eb679f8p+5, 0x1.da1838053b866p+5, 0x1.be886d9c2324dp+6, 0x1.ab514bfc61c76p+7, 0x1.14823229799c2p+7, 0x1.48ff1782ca91dp+8, 0x1.dcbfda0c7559ep+8, From 2f10f9c7fa44ba6f6aeb0bce9b2eefd6387f8979 Mon Sep 17 00:00:00 2001 From: Tue Ly Date: Thu, 10 Oct 2024 05:41:00 +0000 Subject: [PATCH 2/7] Remove excessive `NO_FMA`. --- libc/src/math/generic/cos.cpp | 4 ++-- libc/src/math/generic/sin.cpp | 4 ++-- libc/src/math/generic/sincos.cpp | 8 ++++---- libc/src/math/generic/tan.cpp | 4 ++-- 4 files changed, 10 insertions(+), 10 deletions(-) diff --git a/libc/src/math/generic/cos.cpp b/libc/src/math/generic/cos.cpp index fa3c489aa15ec..923ea96852d88 100644 --- a/libc/src/math/generic/cos.cpp +++ b/libc/src/math/generic/cos.cpp @@ -107,8 +107,8 @@ LLVM_LIBC_FUNCTION(double, cos, (double x)) { // So k is an integer and -pi / 256 <= y <= pi / 256. // Then cos(x) = cos((k * pi/128 + y) // = cos(y) * cos(k*pi/128) - sin(y) * sin(k*pi/128) - DoubleDouble cos_k_cos_y = fputil::quick_mult(cos_y, cos_k); - DoubleDouble msin_k_sin_y = fputil::quick_mult(sin_y, msin_k); + DoubleDouble cos_k_cos_y = fputil::quick_mult(cos_y, cos_k); + DoubleDouble msin_k_sin_y = fputil::quick_mult(sin_y, msin_k); DoubleDouble rr = fputil::exact_add(cos_k_cos_y.hi, msin_k_sin_y.hi); rr.lo += msin_k_sin_y.lo + cos_k_cos_y.lo; diff --git a/libc/src/math/generic/sin.cpp b/libc/src/math/generic/sin.cpp index 4cbcb35513efc..b76b8ee56d34d 100644 --- a/libc/src/math/generic/sin.cpp +++ b/libc/src/math/generic/sin.cpp @@ -118,8 +118,8 @@ LLVM_LIBC_FUNCTION(double, sin, (double x)) { // So k is an integer and -pi / 256 <= y <= pi / 256. // Then sin(x) = sin((k * pi/128 + y) // = sin(y) * cos(k*pi/128) + cos(y) * sin(k*pi/128) - DoubleDouble sin_k_cos_y = fputil::quick_mult(cos_y, sin_k); - DoubleDouble cos_k_sin_y = fputil::quick_mult(sin_y, cos_k); + DoubleDouble sin_k_cos_y = fputil::quick_mult(cos_y, sin_k); + DoubleDouble cos_k_sin_y = fputil::quick_mult(sin_y, cos_k); DoubleDouble rr = fputil::exact_add(sin_k_cos_y.hi, cos_k_sin_y.hi); rr.lo += sin_k_cos_y.lo + cos_k_sin_y.lo; diff --git a/libc/src/math/generic/sincos.cpp b/libc/src/math/generic/sincos.cpp index 314504b3350b5..166ce46603140 100644 --- a/libc/src/math/generic/sincos.cpp +++ b/libc/src/math/generic/sincos.cpp @@ -129,12 +129,12 @@ LLVM_LIBC_FUNCTION(void, sincos, (double x, double *sin_x, double *cos_x)) { // So k is an integer and -pi / 256 <= y <= pi / 256. // Then sin(x) = sin((k * pi/128 + y) // = sin(y) * cos(k*pi/128) + cos(y) * sin(k*pi/128) - DoubleDouble sin_k_cos_y = fputil::quick_mult(cos_y, sin_k); - DoubleDouble cos_k_sin_y = fputil::quick_mult(sin_y, cos_k); + DoubleDouble sin_k_cos_y = fputil::quick_mult(cos_y, sin_k); + DoubleDouble cos_k_sin_y = fputil::quick_mult(sin_y, cos_k); // cos(x) = cos((k * pi/128 + y) // = cos(y) * cos(k*pi/128) - sin(y) * sin(k*pi/128) - DoubleDouble cos_k_cos_y = fputil::quick_mult(cos_y, cos_k); - DoubleDouble msin_k_sin_y = fputil::quick_mult(sin_y, msin_k); + DoubleDouble cos_k_cos_y = fputil::quick_mult(cos_y, cos_k); + DoubleDouble msin_k_sin_y = fputil::quick_mult(sin_y, msin_k); DoubleDouble sin_dd = fputil::exact_add(sin_k_cos_y.hi, cos_k_sin_y.hi); diff --git a/libc/src/math/generic/tan.cpp b/libc/src/math/generic/tan.cpp index 8c19936068dda..eeb859403620b 100644 --- a/libc/src/math/generic/tan.cpp +++ b/libc/src/math/generic/tan.cpp @@ -208,8 +208,8 @@ LLVM_LIBC_FUNCTION(double, tan, (double x)) { // / (cos(y) * cos(k*pi/128) - sin(y) * sin(k*pi/128)) // = (sin(k*pi/128) + tan(y) * cos(k*pi/128)) / // / (cos(k*pi/128) - tan(y) * sin(k*pi/128)) - DoubleDouble cos_k_tan_y = fputil::quick_mult(tan_y, cos_k); - DoubleDouble msin_k_tan_y = fputil::quick_mult(tan_y, msin_k); + DoubleDouble cos_k_tan_y = fputil::quick_mult(tan_y, cos_k); + DoubleDouble msin_k_tan_y = fputil::quick_mult(tan_y, msin_k); // num_dd = sin(k*pi/128) + tan(y) * cos(k*pi/128) DoubleDouble num_dd = fputil::exact_add(cos_k_tan_y.hi, -msin_k.hi); From 6c1cc4c3ad312d74fbf2cc05d35ba3a50fe6304d Mon Sep 17 00:00:00 2001 From: Tue Ly Date: Thu, 10 Oct 2024 15:31:21 +0000 Subject: [PATCH 3/7] Fix formatting. --- libc/src/math/generic/sin.cpp | 3 +-- libc/src/math/generic/tan.cpp | 3 +-- 2 files changed, 2 insertions(+), 4 deletions(-) diff --git a/libc/src/math/generic/sin.cpp b/libc/src/math/generic/sin.cpp index b76b8ee56d34d..2e1d3ffd5f37d 100644 --- a/libc/src/math/generic/sin.cpp +++ b/libc/src/math/generic/sin.cpp @@ -46,13 +46,12 @@ LLVM_LIBC_FUNCTION(double, sin, (double x)) { if (LIBC_LIKELY(x_e < FPBits::EXP_BIAS + FAST_PASS_EXPONENT)) { // |x| < 2^-7 if (LIBC_UNLIKELY(x_e < FPBits::EXP_BIAS - 7)) { - // |x| < 2^-26 + // |x| < 2^-26, |sin(x) - x| < ulp(x)/2. if (LIBC_UNLIKELY(x_e < FPBits::EXP_BIAS - 26)) { // Signed zeros. if (LIBC_UNLIKELY(x == 0.0)) return x; - // For |x| < 2^-26, |sin(x) - x| < ulp(x)/2. #ifdef LIBC_TARGET_CPU_HAS_FMA return fputil::multiply_add(x, -0x1.0p-54, x); #else diff --git a/libc/src/math/generic/tan.cpp b/libc/src/math/generic/tan.cpp index eeb859403620b..f9be25ed866e1 100644 --- a/libc/src/math/generic/tan.cpp +++ b/libc/src/math/generic/tan.cpp @@ -134,13 +134,12 @@ LLVM_LIBC_FUNCTION(double, tan, (double x)) { if (LIBC_LIKELY(x_e < FPBits::EXP_BIAS + FAST_PASS_EXPONENT)) { // |x| < 2^-7 if (LIBC_UNLIKELY(x_e < FPBits::EXP_BIAS - 7)) { - // |x| < 2^-27 + // |x| < 2^-27, |tan(x) - x| < ulp(x)/2. if (LIBC_UNLIKELY(x_e < FPBits::EXP_BIAS - 27)) { // Signed zeros. if (LIBC_UNLIKELY(x == 0.0)) return x; - // For |x| < 2^-27, |tan(x) - x| < ulp(x)/2. #ifdef LIBC_TARGET_CPU_HAS_FMA return fputil::multiply_add(x, 0x1.0p-54, x); #else From cbced2bad7263de8d92220a2072fca019ef9789d Mon Sep 17 00:00:00 2001 From: Tue Ly Date: Thu, 10 Oct 2024 19:12:53 +0000 Subject: [PATCH 4/7] Address comments. --- libc/src/__support/FPUtil/double_double.h | 21 +++++++++---------- libc/src/__support/macros/optimization.h | 4 ++-- .../generic/range_reduction_double_common.h | 7 +++++-- .../math/generic/range_reduction_double_fma.h | 8 +++---- .../generic/range_reduction_double_nofma.h | 14 ++++++------- 5 files changed, 28 insertions(+), 26 deletions(-) diff --git a/libc/src/__support/FPUtil/double_double.h b/libc/src/__support/FPUtil/double_double.h index b51d3c4c63cf5..cef8bd86013a8 100644 --- a/libc/src/__support/FPUtil/double_double.h +++ b/libc/src/__support/FPUtil/double_double.h @@ -18,6 +18,8 @@ namespace LIBC_NAMESPACE_DECL { namespace fputil { +#define DEFAULT_DOUBLE_SPLIT 27 + using DoubleDouble = LIBC_NAMESPACE::NumberPair; // The output of Dekker's FastTwoSum algorithm is correct, i.e.: @@ -61,7 +63,8 @@ LIBC_INLINE constexpr DoubleDouble add(const DoubleDouble &a, double b) { // Zimmermann, P., "Note on the Veltkamp/Dekker Algorithms with Directed // Roundings," https://inria.hal.science/hal-04480440. // Default splitting constant = 2^ceil(prec(double)/2) + 1 = 2^27 + 1. -template LIBC_INLINE constexpr DoubleDouble split(double a) { +template +LIBC_INLINE constexpr DoubleDouble split(double a) { DoubleDouble r{0.0, 0.0}; // CN = 2^N. constexpr double CN = static_cast(1 << N); @@ -73,16 +76,12 @@ template LIBC_INLINE constexpr DoubleDouble split(double a) { return r; } -// Helper for non-fma exact mult where the first number is already splitted. -template +// Helper for non-fma exact mult where the first number is already split. +template LIBC_INLINE DoubleDouble exact_mult(const DoubleDouble &as, double a, double b) { - DoubleDouble bs, r; - - if constexpr (NO_FMA_ALL_ROUNDINGS) - bs = split<28>(b); - else - bs = split(b); + DoubleDouble bs = split(b); + DoubleDouble r{0.0, 0.0}; r.hi = a * b; double t1 = as.hi * bs.hi - r.hi; @@ -100,7 +99,7 @@ LIBC_INLINE DoubleDouble exact_mult(const DoubleDouble &as, double a, // Using Theorem 1 in the paper above, without FMA instruction, if we restrict // the generated constants to precision <= 51, and splitting it by 2^28 + 1, // then a * b = r.hi + r.lo is exact for all rounding modes. -template +template LIBC_INLINE DoubleDouble exact_mult(double a, double b) { DoubleDouble r{0.0, 0.0}; @@ -111,7 +110,7 @@ LIBC_INLINE DoubleDouble exact_mult(double a, double b) { // Dekker's Product. DoubleDouble as = split(a); - r = exact_mult(as, a, b); + r = exact_mult(as, a, b); #endif // LIBC_TARGET_CPU_HAS_FMA return r; diff --git a/libc/src/__support/macros/optimization.h b/libc/src/__support/macros/optimization.h index ce99f1ce6e532..41ecd2bd6d719 100644 --- a/libc/src/__support/macros/optimization.h +++ b/libc/src/__support/macros/optimization.h @@ -50,11 +50,11 @@ LIBC_INLINE constexpr bool expects_bool_condition(T value, T expected) { #define LIBC_MATH 0 #else -#if ((LIBC_MATH & LIBC_MATH_SKIP_ACCURATE_PASS) != 0) +#if (LIBC_MATH & LIBC_MATH_SKIP_ACCURATE_PASS) #define LIBC_MATH_HAS_SKIP_ACCURATE_PASS #endif -#if ((LIBC_MATH & LIBC_MATH_SMALL_TABLES) != 0) +#if (LIBC_MATH & LIBC_MATH_SMALL_TABLES) #define LIBC_MATH_HAS_SMALL_TABLES #endif diff --git a/libc/src/math/generic/range_reduction_double_common.h b/libc/src/math/generic/range_reduction_double_common.h index 9469ecebc5580..1b9c5dfb2f071 100644 --- a/libc/src/math/generic/range_reduction_double_common.h +++ b/libc/src/math/generic/range_reduction_double_common.h @@ -22,9 +22,12 @@ namespace LIBC_NAMESPACE_DECL { #ifdef LIBC_TARGET_CPU_HAS_FMA -static constexpr bool NO_FMA = false; +static constexpr unsigned SPLIT = DEFAULT_DOUBLE_SPLIT; #else -static constexpr bool NO_FMA = true; +// When there is no-FMA instructions, in order to have exact product of 2 double +// precision with directional roundings, we need to lower the precision of the +// constants by at least 1 bit, and use a different splitting constant. +static constexpr unsigned SPLIT = 28; #endif // LIBC_TARGET_CPU_HAS_FMA using LIBC_NAMESPACE::fputil::DoubleDouble; diff --git a/libc/src/math/generic/range_reduction_double_fma.h b/libc/src/math/generic/range_reduction_double_fma.h index d220bef4cd349..cab031c28baa1 100644 --- a/libc/src/math/generic/range_reduction_double_fma.h +++ b/libc/src/math/generic/range_reduction_double_fma.h @@ -34,13 +34,13 @@ LIBC_INLINE unsigned LargeRangeReduction::fast(double x, DoubleDouble &u) { x_reduced = xbits.get_val(); // x * c_hi = ph.hi + ph.lo exactly. DoubleDouble ph = - fputil::exact_mult(x_reduced, ONE_TWENTY_EIGHT_OVER_PI[idx][0]); + fputil::exact_mult(x_reduced, ONE_TWENTY_EIGHT_OVER_PI[idx][0]); // x * c_mid = pm.hi + pm.lo exactly. DoubleDouble pm = - fputil::exact_mult(x_reduced, ONE_TWENTY_EIGHT_OVER_PI[idx][1]); + fputil::exact_mult(x_reduced, ONE_TWENTY_EIGHT_OVER_PI[idx][1]); // x * c_lo = pl.hi + pl.lo exactly. DoubleDouble pl = - fputil::exact_mult(x_reduced, ONE_TWENTY_EIGHT_OVER_PI[idx][2]); + fputil::exact_mult(x_reduced, ONE_TWENTY_EIGHT_OVER_PI[idx][2]); // Extract integral parts and fractional parts of (ph.lo + pm.hi). double sum_hi = ph.lo + pm.hi; double kd = fputil::nearest_integer(sum_hi); @@ -67,7 +67,7 @@ LIBC_INLINE unsigned LargeRangeReduction::fast(double x, DoubleDouble &u) { // Then, // | {x * 128/pi} - (y_hi + y_lo) | <= ulp(ulp(y_hi)) <= 2^-105 // | {x mod pi/128} - (u.hi + u.lo) | < 2 * 2^-6 * 2^-105 = 2^-110 - u = fputil::quick_mult(y, PI_OVER_128_DD); + u = fputil::quick_mult(y, PI_OVER_128_DD); return static_cast(static_cast(kd)); } diff --git a/libc/src/math/generic/range_reduction_double_nofma.h b/libc/src/math/generic/range_reduction_double_nofma.h index 4ce35556b1ee8..5640732947798 100644 --- a/libc/src/math/generic/range_reduction_double_nofma.h +++ b/libc/src/math/generic/range_reduction_double_nofma.h @@ -34,14 +34,14 @@ LIBC_INLINE unsigned LargeRangeReduction::fast(double x, DoubleDouble &u) { x_reduced = xbits.get_val(); // x * c_hi = ph.hi + ph.lo exactly. DoubleDouble x_split = fputil::split(x_reduced); - DoubleDouble ph = fputil::exact_mult( - x_split, x_reduced, ONE_TWENTY_EIGHT_OVER_PI[idx][0]); + DoubleDouble ph = fputil::exact_mult(x_split, x_reduced, + ONE_TWENTY_EIGHT_OVER_PI[idx][0]); // x * c_mid = pm.hi + pm.lo exactly. - DoubleDouble pm = fputil::exact_mult( - x_split, x_reduced, ONE_TWENTY_EIGHT_OVER_PI[idx][1]); + DoubleDouble pm = fputil::exact_mult(x_split, x_reduced, + ONE_TWENTY_EIGHT_OVER_PI[idx][1]); // x * c_lo = pl.hi + pl.lo exactly. - DoubleDouble pl = fputil::exact_mult( - x_split, x_reduced, ONE_TWENTY_EIGHT_OVER_PI[idx][2]); + DoubleDouble pl = fputil::exact_mult(x_split, x_reduced, + ONE_TWENTY_EIGHT_OVER_PI[idx][2]); // Extract integral parts and fractional parts of (ph.lo + pm.hi). double sum_hi = ph.lo + pm.hi; double kd = fputil::nearest_integer(sum_hi); @@ -68,7 +68,7 @@ LIBC_INLINE unsigned LargeRangeReduction::fast(double x, DoubleDouble &u) { // Then, // | {x * 128/pi} - (y_hi + y_lo) | <= ulp(ulp(y_hi)) <= 2^-105 // | {x mod pi/128} - (u.hi + u.lo) | < 2 * 2^-6 * 2^-105 = 2^-110 - u = fputil::quick_mult(y, PI_OVER_128_DD); + u = fputil::quick_mult(y, PI_OVER_128_DD); return static_cast(static_cast(kd)); } From ed5f0662c11793b77ffe5c0d2e1fae22e3350c29 Mon Sep 17 00:00:00 2001 From: Tue Ly Date: Thu, 10 Oct 2024 19:28:14 +0000 Subject: [PATCH 5/7] Remove extra ; in comments. --- libc/src/math/generic/range_reduction_double_common.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/libc/src/math/generic/range_reduction_double_common.h b/libc/src/math/generic/range_reduction_double_common.h index 1b9c5dfb2f071..e23bbff144bee 100644 --- a/libc/src/math/generic/range_reduction_double_common.h +++ b/libc/src/math/generic/range_reduction_double_common.h @@ -63,7 +63,7 @@ LIBC_INLINE unsigned range_reduction_small(double x, DoubleDouble &u) { double y_hi = fputil::multiply_add(kd, MPI_OVER_128[0], x); // Exact // |u.hi| < 1.6*2^-7 u.hi = fputil::multiply_add(kd, MPI_OVER_128[1], y_hi); - double u0 = y_hi - u.hi; // Exact; + double u0 = y_hi - u.hi; // Exact // |u.lo| <= max(ulp(u.hi), |kd * MPI_OVER_128[2]|) double u1 = fputil::multiply_add(kd, MPI_OVER_128[1], u0); // Exact u.lo = fputil::multiply_add(kd, MPI_OVER_128[2], u1); From 56d5c65cec7f71bef512cd3d5215743856add836 Mon Sep 17 00:00:00 2001 From: Tue Ly Date: Thu, 10 Oct 2024 19:37:04 +0000 Subject: [PATCH 6/7] Fix quick_mult. --- libc/src/__support/FPUtil/double_double.h | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/libc/src/__support/FPUtil/double_double.h b/libc/src/__support/FPUtil/double_double.h index cef8bd86013a8..db3c2c8a3d7a6 100644 --- a/libc/src/__support/FPUtil/double_double.h +++ b/libc/src/__support/FPUtil/double_double.h @@ -122,10 +122,10 @@ LIBC_INLINE DoubleDouble quick_mult(double a, const DoubleDouble &b) { return r; } -template +template LIBC_INLINE DoubleDouble quick_mult(const DoubleDouble &a, const DoubleDouble &b) { - DoubleDouble r = exact_mult(a.hi, b.hi); + DoubleDouble r = exact_mult(a.hi, b.hi); double t1 = multiply_add(a.hi, b.lo, r.lo); double t2 = multiply_add(a.lo, b.hi, t1); r.lo = t2; From 44505905763f5144e4bda1b1465e604aa1111d24 Mon Sep 17 00:00:00 2001 From: Tue Ly Date: Fri, 11 Oct 2024 02:38:16 +0000 Subject: [PATCH 7/7] Fix exact_mult call in pow.cpp. --- libc/src/math/generic/pow.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/libc/src/math/generic/pow.cpp b/libc/src/math/generic/pow.cpp index 3a50e220154e5..181d3d40b3c9a 100644 --- a/libc/src/math/generic/pow.cpp +++ b/libc/src/math/generic/pow.cpp @@ -398,7 +398,7 @@ LLVM_LIBC_FUNCTION(double, pow, (double x, double y)) { #else double c = FPBits(m_x.uintval() & 0x3fff'e000'0000'0000).get_val(); dx = fputil::multiply_add(RD[idx_x], m_x.get_val() - c, CD[idx_x]); // Exact - dx_c0 = fputil::exact_mult(COEFFS[0], dx); + dx_c0 = fputil::exact_mult<28>(dx, COEFFS[0]); // Exact #endif // LIBC_TARGET_CPU_HAS_FMA double dx2 = dx * dx;