[libc-commits] [libc] 51e9430 - [libc][math] Improve performance of double precision trig functions. (#111793)
via libc-commits
libc-commits at lists.llvm.org
Thu Oct 10 20:33:07 PDT 2024
Author: lntue
Date: 2024-10-10T23:33:02-04:00
New Revision: 51e9430a0c767243411d4b81c284700f89719277
URL: https://github.com/llvm/llvm-project/commit/51e9430a0c767243411d4b81c284700f89719277
DIFF: https://github.com/llvm/llvm-project/commit/51e9430a0c767243411d4b81c284700f89719277.diff
LOG: [libc][math] Improve performance of double precision trig functions. (#111793)
- Improve the accuracy of fast pass' range reduction.
- Provide tighter error estimations.
- Reduce the table size when `LIBC_MATH_SMALL_TABLES` flag is set.
Added:
Modified:
libc/src/__support/FPUtil/double_double.h
libc/src/__support/macros/optimization.h
libc/src/math/generic/cos.cpp
libc/src/math/generic/pow.cpp
libc/src/math/generic/range_reduction_double_common.h
libc/src/math/generic/range_reduction_double_fma.h
libc/src/math/generic/range_reduction_double_nofma.h
libc/src/math/generic/sin.cpp
libc/src/math/generic/sincos.cpp
libc/src/math/generic/sincos_eval.h
libc/src/math/generic/tan.cpp
libc/test/src/math/cos_test.cpp
libc/test/src/math/sin_test.cpp
libc/test/src/math/tan_test.cpp
Removed:
################################################################################
diff --git a/libc/src/__support/FPUtil/double_double.h b/libc/src/__support/FPUtil/double_double.h
index 25a4ee03387c67..db3c2c8a3d7a6e 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<double>;
// 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 <size_t N = 27> LIBC_INLINE constexpr DoubleDouble split(double a) {
+template <size_t N = DEFAULT_DOUBLE_SPLIT>
+LIBC_INLINE constexpr DoubleDouble split(double a) {
DoubleDouble r{0.0, 0.0};
// CN = 2^N.
constexpr double CN = static_cast<double>(1 << N);
@@ -73,6 +76,22 @@ template <size_t N = 27> LIBC_INLINE constexpr DoubleDouble split(double a) {
return r;
}
+// Helper for non-fma exact mult where the first number is already split.
+template <size_t SPLIT_B = DEFAULT_DOUBLE_SPLIT>
+LIBC_INLINE DoubleDouble exact_mult(const DoubleDouble &as, double a,
+ double b) {
+ DoubleDouble bs = split<SPLIT_B>(b);
+ DoubleDouble r{0.0, 0.0};
+
+ 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
@@ -80,7 +99,7 @@ template <size_t N = 27> LIBC_INLINE constexpr DoubleDouble split(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 <bool NO_FMA_ALL_ROUNDINGS = false>
+template <size_t SPLIT_B = 27>
LIBC_INLINE DoubleDouble exact_mult(double a, double b) {
DoubleDouble r{0.0, 0.0};
@@ -90,18 +109,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<SPLIT_B>(as, a, b);
#endif // LIBC_TARGET_CPU_HAS_FMA
return r;
@@ -113,10 +122,10 @@ LIBC_INLINE DoubleDouble quick_mult(double a, const DoubleDouble &b) {
return r;
}
-template <bool NO_FMA_ALL_ROUNDINGS = false>
+template <size_t SPLIT_B = 27>
LIBC_INLINE DoubleDouble quick_mult(const DoubleDouble &a,
const DoubleDouble &b) {
- DoubleDouble r = exact_mult<NO_FMA_ALL_ROUNDINGS>(a.hi, b.hi);
+ DoubleDouble r = exact_mult<SPLIT_B>(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;
@@ -157,8 +166,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</*NO_FMA=*/true>(b.hi, -r.hi);
- DoubleDouble b_lo_r_hi = fputil::exact_mult</*NO_FMA=*/true>(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 5ffd474d35c54d..41ecd2bd6d7191 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)
+#define LIBC_MATH_HAS_SKIP_ACCURATE_PASS
+#endif
+
+#if (LIBC_MATH & LIBC_MATH_SMALL_TABLES)
+#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 e61d800ce2dada..923ea96852d889 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<NO_FMA> 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,70 +74,51 @@ 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<uint64_t>((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<uint64_t>((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.
// 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<NO_FMA>(cos_y, cos_k);
- DoubleDouble msin_k_sin_y = fputil::quick_mult<NO_FMA>(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<false>(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<double>(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/pow.cpp b/libc/src/math/generic/pow.cpp
index 3a50e220154e51..181d3d40b3c9ad 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<true>(COEFFS[0], dx);
+ dx_c0 = fputil::exact_mult<28>(dx, COEFFS[0]); // Exact
#endif // LIBC_TARGET_CPU_HAS_FMA
double dx2 = dx * dx;
diff --git a/libc/src/math/generic/range_reduction_double_common.h b/libc/src/math/generic/range_reduction_double_common.h
index 290b642be4c69f..e23bbff144bee8 100644
--- a/libc/src/math/generic/range_reduction_double_common.h
+++ b/libc/src/math/generic/range_reduction_double_common.h
@@ -17,150 +17,272 @@
#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 unsigned SPLIT = DEFAULT_DOUBLE_SPLIT;
#else
-#include "range_reduction_double_nofma.h"
+// 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
diff erent splitting constant.
+static constexpr unsigned SPLIT = 28;
+#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<unsigned>(static_cast<int64_t>(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 <bool NO_FMA> struct LargeRangeReduction {
- // Calculate the high part of the range reduction exactly.
- LIBC_INLINE unsigned compute_high_part(double x) {
- using FPBits = typename fputil::FPBits<double>;
- 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<unsigned>((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<NO_FMA>(x_reduced, ONE_TWENTY_EIGHT_OVER_PI[idx][0]);
- // x * c_mid = pm.hi + pm.lo exactly.
- DoubleDouble pm =
- fputil::exact_mult<NO_FMA>(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<unsigned>(static_cast<int64_t>(kh) +
- static_cast<int64_t>(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<NO_FMA>(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 +300,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 +368,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 7448b5f63dfde2..cab031c28baa17 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<double>;
+ 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<unsigned>((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<SPLIT>(x_reduced, ONE_TWENTY_EIGHT_OVER_PI[idx][0]);
+ // x * c_mid = pm.hi + pm.lo exactly.
+ DoubleDouble pm =
+ fputil::exact_mult<SPLIT>(x_reduced, ONE_TWENTY_EIGHT_OVER_PI[idx][1]);
+ // x * c_lo = pl.hi + pl.lo exactly.
+ DoubleDouble pl =
+ fputil::exact_mult<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<SPLIT>(y, PI_OVER_128_DD);
+
+ return static_cast<unsigned>(static_cast<int64_t>(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<unsigned>(static_cast<int64_t>(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 445a45d3f9796a..56407329477989 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<double>;
+ 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<unsigned>((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<SPLIT>(x_split, x_reduced,
+ ONE_TWENTY_EIGHT_OVER_PI[idx][0]);
+ // x * c_mid = pm.hi + pm.lo exactly.
+ DoubleDouble pm = fputil::exact_mult<SPLIT>(x_split, x_reduced,
+ ONE_TWENTY_EIGHT_OVER_PI[idx][1]);
+ // x * c_lo = pl.hi + pl.lo exactly.
+ DoubleDouble pl = fputil::exact_mult<SPLIT>(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<SPLIT>(y, PI_OVER_128_DD);
+
+ return static_cast<unsigned>(static_cast<int64_t>(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<unsigned>(static_cast<int>(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 da3d1e94b5f645..2e1d3ffd5f37d8 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,39 @@ LLVM_LIBC_FUNCTION(double, sin, (double x)) {
DoubleDouble y;
unsigned k;
- generic::LargeRangeReduction<NO_FMA> 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, |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);
+ 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,69 +85,51 @@ 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<uint64_t>(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<uint64_t>((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.
// 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<NO_FMA>(cos_y, sin_k);
- DoubleDouble cos_k_sin_y = fputil::quick_mult<NO_FMA>(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<false>(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 +140,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 +148,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 +167,7 @@ LLVM_LIBC_FUNCTION(double, sin, (double x)) {
// https://github.com/llvm/llvm-project/issues/96452.
return static_cast<double>(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 1af0ee7b0eb2c8..166ce466031409 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<NO_FMA> 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,56 +95,46 @@ 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<uint64_t>(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<uint64_t>((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).
// 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<NO_FMA>(cos_y, sin_k);
- DoubleDouble cos_k_sin_y = fputil::quick_mult<NO_FMA>(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<NO_FMA>(cos_y, cos_k);
- DoubleDouble msin_k_sin_y = fputil::quick_mult<NO_FMA>(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<false>(sin_k_cos_y.hi, cos_k_sin_y.hi);
@@ -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 e491467c5663fd..6cd1da4721bf57 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<double>(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 45fd8bb9156be0..f9be25ed866e1d 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<double>(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,38 @@ LLVM_LIBC_FUNCTION(double, tan, (double x)) {
DoubleDouble y;
unsigned k;
- generic::LargeRangeReduction<NO_FMA> 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, |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);
+ 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 +172,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<uint64_t>(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<uint64_t>((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.
@@ -212,8 +207,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<NO_FMA>(tan_y, cos_k);
- DoubleDouble msin_k_tan_y = fputil::quick_mult<NO_FMA>(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<false>(cos_k_tan_y.hi, -msin_k.hi);
@@ -222,7 +217,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 +226,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<uint64_t>(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<int>((k & 63) == 0) * y_exp;
- int64_t tan_x_err = static_cast<int64_t>(FPBits(tan_x.hi).uintval()) -
- (static_cast<int64_t>(rel_err_exp) << 52);
- double tan_err = FPBits(static_cast<uint64_t>(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 +249,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 +257,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 +285,7 @@ LLVM_LIBC_FUNCTION(double, tan, (double x)) {
// https://github.com/llvm/llvm-project/issues/96452.
return static_cast<double>(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 484d47fd3e96c4..e2d47917e545e0 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 60f6ef5c844630..d4c6bd416a4099 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 1ca67afdaddf25..12dfc02bac111a 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,
More information about the libc-commits
mailing list