aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorlntue <35648136+lntue@users.noreply.github.com>2024-06-24 17:57:08 -0400
committerGitHub <noreply@github.com>2024-06-24 17:57:08 -0400
commit16903ace180755b7558234ff2b2e8d89b00dcb88 (patch)
tree1158c0aef94cca9903d1238e0eef7abdf3e4cb91
parent5ae50698a0d6a3022af2e79d405a7eb6c8c790f0 (diff)
downloadllvm-16903ace180755b7558234ff2b2e8d89b00dcb88.zip
llvm-16903ace180755b7558234ff2b2e8d89b00dcb88.tar.gz
llvm-16903ace180755b7558234ff2b2e8d89b00dcb88.tar.bz2
[libc][math] Implement double precision sin correctly rounded to all rounding modes. (#95736)
- Algorithm: - Step 1 - Range reduction: for a double precision input `x`, return `k` and `u` such that - k is an integer - u = x - k * pi / 128, and |u| < pi/256 - Step 2 - Calculate `sin(u)` and `cos(u)` in double-double using Taylor polynomials with errors < 2^-70 with FMA or < 2^-66 w/o FMA. - Step 3 - Calculate `sin(x) = sin(k*pi/128) * cos(u) + cos(k*pi/128) * sin(u)` using look-up table for `sin(k*pi/128)` and `cos(k*pi/128)`. - Step 4 - Use Ziv's rounding test to decide if the result is correctly rounded. - Step 4' - If the Ziv's rounding test failed, redo step 1-3 using 128-bit precision. - Currently, without FMA instructions, the large range reduction only works correctly for the default rounding mode (FE_TONEAREST). - Provide `LIBC_MATH` flag so that users can set `LIBC_MATH = LIBC_MATH_SKIP_ACCURATE_PASS` to build the `sin` function without step 4 and 4'.
Diffstat
-rw-r--r--libc/config/darwin/arm/entrypoints.txt1
-rw-r--r--libc/config/linux/aarch64/entrypoints.txt1
-rw-r--r--libc/config/linux/arm/entrypoints.txt1
-rw-r--r--libc/config/linux/riscv/entrypoints.txt1
-rw-r--r--libc/docs/math/index.rst2
-rw-r--r--libc/src/__support/FPUtil/double_double.h52
-rw-r--r--libc/src/__support/FPUtil/dyadic_float.h10
-rw-r--r--libc/src/__support/macros/optimization.h14
-rw-r--r--libc/src/math/generic/CMakeLists.txt49
-rw-r--r--libc/src/math/generic/range_reduction_double_common.h162
-rw-r--r--libc/src/math/generic/range_reduction_double_fma.h495
-rw-r--r--libc/src/math/generic/range_reduction_double_nofma.h493
-rw-r--r--libc/src/math/generic/sin.cpp315
-rw-r--r--libc/src/math/generic/sincos_eval.h81
-rw-r--r--libc/src/math/x86_64/CMakeLists.txt10
-rw-r--r--libc/src/math/x86_64/sin.cpp19
-rw-r--r--libc/test/src/math/sin_test.cpp107
-rw-r--r--libc/test/src/math/smoke/CMakeLists.txt10
-rw-r--r--libc/test/src/math/smoke/sin_test.cpp26
19 files changed, 1790 insertions, 59 deletions
diff --git a/libc/config/darwin/arm/entrypoints.txt b/libc/config/darwin/arm/entrypoints.txt
index e130326..d486916 100644
--- a/libc/config/darwin/arm/entrypoints.txt
+++ b/libc/config/darwin/arm/entrypoints.txt
@@ -226,6 +226,7 @@ set(TARGET_LIBM_ENTRYPOINTS
libc.src.math.scalbnl
libc.src.math.sincosf
libc.src.math.sinhf
+ libc.src.math.sin
libc.src.math.sinf
libc.src.math.sqrt
libc.src.math.sqrtf
diff --git a/libc/config/linux/aarch64/entrypoints.txt b/libc/config/linux/aarch64/entrypoints.txt
index 63238ec..2d2e1b5 100644
--- a/libc/config/linux/aarch64/entrypoints.txt
+++ b/libc/config/linux/aarch64/entrypoints.txt
@@ -481,6 +481,7 @@ set(TARGET_LIBM_ENTRYPOINTS
libc.src.math.scalbnl
libc.src.math.sincosf
libc.src.math.sinhf
+ libc.src.math.sin
libc.src.math.sinf
libc.src.math.sqrt
libc.src.math.sqrtf
diff --git a/libc/config/linux/arm/entrypoints.txt b/libc/config/linux/arm/entrypoints.txt
index 00df445..53cdcc3 100644
--- a/libc/config/linux/arm/entrypoints.txt
+++ b/libc/config/linux/arm/entrypoints.txt
@@ -354,6 +354,7 @@ set(TARGET_LIBM_ENTRYPOINTS
libc.src.math.scalbnf
libc.src.math.scalbnl
libc.src.math.sincosf
+ libc.src.math.sin
libc.src.math.sinf
libc.src.math.sinhf
libc.src.math.sqrt
diff --git a/libc/config/linux/riscv/entrypoints.txt b/libc/config/linux/riscv/entrypoints.txt
index e12d6b3..18968f5 100644
--- a/libc/config/linux/riscv/entrypoints.txt
+++ b/libc/config/linux/riscv/entrypoints.txt
@@ -489,6 +489,7 @@ set(TARGET_LIBM_ENTRYPOINTS
libc.src.math.scalbnl
libc.src.math.sincosf
libc.src.math.sinhf
+ libc.src.math.sin
libc.src.math.sinf
libc.src.math.sqrt
libc.src.math.sqrtf
diff --git a/libc/docs/math/index.rst b/libc/docs/math/index.rst
index d6e642c..29723be 100644
--- a/libc/docs/math/index.rst
+++ b/libc/docs/math/index.rst
@@ -322,7 +322,7 @@ Higher Math Functions
+-----------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+----------------------------+
| rsqrt | | | | | | 7.12.7.9 | F.10.4.9 |
+-----------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+----------------------------+
-| sin | |check| | large | | | | 7.12.4.6 | F.10.1.6 |
+| sin | |check| | |check| | | | | 7.12.4.6 | F.10.1.6 |
+-----------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+----------------------------+
| sincos | |check| | large | | | | | |
+-----------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+----------------------------+
diff --git a/libc/src/__support/FPUtil/double_double.h b/libc/src/__support/FPUtil/double_double.h
index b9490b5..3d16a3c 100644
--- a/libc/src/__support/FPUtil/double_double.h
+++ b/libc/src/__support/FPUtil/double_double.h
@@ -21,12 +21,22 @@ using DoubleDouble = LIBC_NAMESPACE::NumberPair<double>;
// The output of Dekker's FastTwoSum algorithm is correct, i.e.:
// r.hi + r.lo = a + b exactly
// and |r.lo| < eps(r.lo)
-// if ssumption: |a| >= |b|, or a = 0.
+// Assumption: |a| >= |b|, or a = 0.
+template <bool FAST2SUM = true>
LIBC_INLINE constexpr DoubleDouble exact_add(double a, double b) {
DoubleDouble r{0.0, 0.0};
- r.hi = a + b;
- double t = r.hi - a;
- r.lo = b - t;
+ if constexpr (FAST2SUM) {
+ r.hi = a + b;
+ double t = r.hi - a;
+ r.lo = b - t;
+ } else {
+ r.hi = a + b;
+ double t1 = r.hi - a;
+ double t2 = r.hi - t1;
+ double t3 = b - t1;
+ double t4 = a - t2;
+ r.lo = t3 + t4;
+ }
return r;
}
@@ -40,15 +50,20 @@ LIBC_INLINE constexpr DoubleDouble add(const DoubleDouble &a,
// Assumption: |a.hi| >= |b|
LIBC_INLINE constexpr DoubleDouble add(const DoubleDouble &a, double b) {
- DoubleDouble r = exact_add(a.hi, b);
+ DoubleDouble r = exact_add<false>(a.hi, b);
return exact_add(r.hi, r.lo + a.lo);
}
-// Velkamp's Splitting for double precision.
-LIBC_INLINE constexpr DoubleDouble split(double a) {
+// Veltkamp's Splitting for double precision.
+// Note: This is proved to be correct for all rounding modes:
+// 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) {
DoubleDouble r{0.0, 0.0};
- // Splitting constant = 2^ceil(prec(double)/2) + 1 = 2^27 + 1.
- constexpr double C = 0x1.0p27 + 1.0;
+ // CN = 2^N.
+ constexpr double CN = static_cast<double>(1 << N);
+ constexpr double C = CN + 1.0;
double t1 = C * a;
double t2 = a - t1;
r.hi = t1 + t2;
@@ -56,6 +71,14 @@ LIBC_INLINE constexpr DoubleDouble split(double a) {
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
+// Roundings," https://inria.hal.science/hal-04480440.
+// 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>
LIBC_INLINE DoubleDouble exact_mult(double a, double b) {
DoubleDouble r{0.0, 0.0};
@@ -65,7 +88,13 @@ LIBC_INLINE DoubleDouble exact_mult(double a, double b) {
#else
// Dekker's Product.
DoubleDouble as = split(a);
- DoubleDouble bs = split(b);
+ 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;
@@ -82,9 +111,10 @@ LIBC_INLINE DoubleDouble quick_mult(double a, const DoubleDouble &b) {
return r;
}
+template <bool NO_FMA_ALL_ROUNDINGS = false>
LIBC_INLINE DoubleDouble quick_mult(const DoubleDouble &a,
const DoubleDouble &b) {
- DoubleDouble r = exact_mult(a.hi, b.hi);
+ DoubleDouble r = exact_mult<NO_FMA_ALL_ROUNDINGS>(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;
diff --git a/libc/src/__support/FPUtil/dyadic_float.h b/libc/src/__support/FPUtil/dyadic_float.h
index 63cb983..76786a2 100644
--- a/libc/src/__support/FPUtil/dyadic_float.h
+++ b/libc/src/__support/FPUtil/dyadic_float.h
@@ -278,11 +278,11 @@ LIBC_INLINE constexpr DyadicFloat<Bits> quick_add(DyadicFloat<Bits> a,
// don't need to normalize the inputs again in this function. If the inputs are
// not normalized, the results might lose precision significantly.
template <size_t Bits>
-LIBC_INLINE constexpr DyadicFloat<Bits> quick_mul(DyadicFloat<Bits> a,
- DyadicFloat<Bits> b) {
+LIBC_INLINE constexpr DyadicFloat<Bits> quick_mul(const DyadicFloat<Bits> &a,
+ const DyadicFloat<Bits> &b) {
DyadicFloat<Bits> result;
result.sign = (a.sign != b.sign) ? Sign::NEG : Sign::POS;
- result.exponent = a.exponent + b.exponent + int(Bits);
+ result.exponent = a.exponent + b.exponent + static_cast<int>(Bits);
if (!(a.mantissa.is_zero() || b.mantissa.is_zero())) {
result.mantissa = a.mantissa.quick_mul_hi(b.mantissa);
@@ -309,7 +309,7 @@ multiply_add(const DyadicFloat<Bits> &a, const DyadicFloat<Bits> &b,
// Simple exponentiation implementation for printf. Only handles positive
// exponents, since division isn't implemented.
template <size_t Bits>
-LIBC_INLINE constexpr DyadicFloat<Bits> pow_n(DyadicFloat<Bits> a,
+LIBC_INLINE constexpr DyadicFloat<Bits> pow_n(const DyadicFloat<Bits> &a,
uint32_t power) {
DyadicFloat<Bits> result = 1.0;
DyadicFloat<Bits> cur_power = a;
@@ -325,7 +325,7 @@ LIBC_INLINE constexpr DyadicFloat<Bits> pow_n(DyadicFloat<Bits> a,
}
template <size_t Bits>
-LIBC_INLINE constexpr DyadicFloat<Bits> mul_pow_2(DyadicFloat<Bits> a,
+LIBC_INLINE constexpr DyadicFloat<Bits> mul_pow_2(const DyadicFloat<Bits> &a,
int32_t pow_2) {
DyadicFloat<Bits> result = a;
result.exponent += pow_2;
diff --git a/libc/src/__support/macros/optimization.h b/libc/src/__support/macros/optimization.h
index 59886ca..05a47791 100644
--- a/libc/src/__support/macros/optimization.h
+++ b/libc/src/__support/macros/optimization.h
@@ -33,4 +33,18 @@ LIBC_INLINE constexpr bool expects_bool_condition(T value, T expected) {
#error "Unhandled compiler"
#endif
+// Defining optimization options for math functions.
+// TODO: Exporting this to public generated headers?
+#define LIBC_MATH_SKIP_ACCURATE_PASS 0x01
+#define LIBC_MATH_SMALL_TABLES 0x02
+#define LIBC_MATH_NO_ERRNO 0x04
+#define LIBC_MATH_NO_EXCEPT 0x08
+#define LIBC_MATH_FAST \
+ (LIBC_MATH_SKIP_ACCURATE_PASS | LIBC_MATH_SMALL_TABLES | \
+ LIBC_MATH_NO_ERRNO | LIBC_MATH_NO_EXCEPT)
+
+#ifndef LIBC_MATH
+#define LIBC_MATH 0
+#endif // LIBC_MATH
+
#endif // LLVM_LIBC_SRC___SUPPORT_MACROS_OPTIMIZATION_H
diff --git a/libc/src/math/generic/CMakeLists.txt b/libc/src/math/generic/CMakeLists.txt
index a0114aa..54a5b6a 100644
--- a/libc/src/math/generic/CMakeLists.txt
+++ b/libc/src/math/generic/CMakeLists.txt
@@ -136,6 +136,23 @@ add_header_library(
)
add_header_library(
+ range_reduction_double
+ HDRS
+ range_reduction_double_common.h
+ range_reduction_double_fma.h
+ range_reduction_double_nofma.h
+ DEPENDS
+ libc.src.__support.FPUtil.double_double
+ libc.src.__support.FPUtil.dyadic_float
+ libc.src.__support.FPUtil.fp_bits
+ libc.src.__support.FPUtil.fma
+ libc.src.__support.FPUtil.multiply_add
+ libc.src.__support.FPUtil.nearest_integer
+ libc.src.__support.common
+ libc.src.__support.integer_literals
+)
+
+add_header_library(
sincosf_utils
HDRS
sincosf_utils.h
@@ -146,6 +163,15 @@ add_header_library(
libc.src.__support.common
)
+add_header_library(
+ sincos_eval
+ HDRS
+ sincos_eval.h
+ DEPENDS
+ libc.src.__support.FPUtil.double_double
+ libc.src.__support.FPUtil.multiply_add
+)
+
add_entrypoint_object(
cosf
SRCS
@@ -168,6 +194,29 @@ add_entrypoint_object(
)
add_entrypoint_object(
+ sin
+ SRCS
+ sin.cpp
+ HDRS
+ ../sin.h
+ DEPENDS
+ libc.hdr.errno_macros
+ libc.src.errno.errno
+ libc.src.__support.FPUtil.double_double
+ libc.src.__support.FPUtil.dyadic_float
+ libc.src.__support.FPUtil.fenv_impl
+ libc.src.__support.FPUtil.fp_bits
+ libc.src.__support.FPUtil.fma
+ libc.src.__support.FPUtil.multiply_add
+ libc.src.__support.FPUtil.nearest_integer
+ libc.src.__support.FPUtil.polyeval
+ libc.src.__support.FPUtil.rounding_mode
+ libc.src.__support.macros.optimization
+ COMPILE_OPTIONS
+ -O3
+)
+
+add_entrypoint_object(
sinf
SRCS
sinf.cpp
diff --git a/libc/src/math/generic/range_reduction_double_common.h b/libc/src/math/generic/range_reduction_double_common.h
new file mode 100644
index 0000000..0e9edf8
--- /dev/null
+++ b/libc/src/math/generic/range_reduction_double_common.h
@@ -0,0 +1,162 @@
+//===-- Range reduction for double precision sin/cos/tan -*- C++ --------*-===//
+//
+// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
+// See https://llvm.org/LICENSE.txt for license information.
+// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
+//
+//===----------------------------------------------------------------------===//
+
+#ifndef LLVM_LIBC_SRC_MATH_GENERIC_RANGE_REDUCTION_DOUBLE_COMMON_H
+#define LLVM_LIBC_SRC_MATH_GENERIC_RANGE_REDUCTION_DOUBLE_COMMON_H
+
+#include "src/__support/FPUtil/FPBits.h"
+#include "src/__support/FPUtil/double_double.h"
+#include "src/__support/FPUtil/dyadic_float.h"
+#include "src/__support/FPUtil/multiply_add.h"
+#include "src/__support/FPUtil/nearest_integer.h"
+#include "src/__support/common.h"
+#include "src/__support/integer_literals.h"
+
+namespace LIBC_NAMESPACE {
+
+namespace generic {
+
+using LIBC_NAMESPACE::fputil::DoubleDouble;
+using Float128 = LIBC_NAMESPACE::fputil::DyadicFloat<128>;
+
+LIBC_INLINE constexpr Float128 PI_OVER_128_F128 = {
+ Sign::POS, -133, 0xc90f'daa2'2168'c234'c4c6'628b'80dc'1cd1_u128};
+
+// Note: The look-up tables ONE_TWENTY_EIGHT_OVER_PI is selected to be either
+// from fma:: or nofma:: namespace.
+
+// 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,
+// 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)
+// 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));
+ }
+
+ 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);
+ }
+
+ LIBC_INLINE Float128 accurate() const {
+ // 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)));
+
+ return fputil::quick_mul(y, PI_OVER_128_F128);
+ }
+
+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;
+};
+
+LIBC_INLINE Float128 range_reduction_small_f128(double x) {
+ double prod_hi = x * ONE_TWENTY_EIGHT_OVER_PI[3][0];
+ double kd = fputil::nearest_integer(prod_hi);
+
+ Float128 mk_f128(-kd);
+ Float128 x_f128(x);
+ Float128 p_hi =
+ fputil::quick_mul(x_f128, Float128(ONE_TWENTY_EIGHT_OVER_PI[3][0]));
+ Float128 p_mid =
+ fputil::quick_mul(x_f128, Float128(ONE_TWENTY_EIGHT_OVER_PI[3][1]));
+ Float128 p_lo =
+ fputil::quick_mul(x_f128, Float128(ONE_TWENTY_EIGHT_OVER_PI[3][2]));
+ Float128 s_hi = fputil::quick_add(p_hi, mk_f128);
+ Float128 s_lo = fputil::quick_add(p_mid, p_lo);
+ Float128 y = fputil::quick_add(s_hi, s_lo);
+
+ return fputil::quick_mul(y, PI_OVER_128_F128);
+}
+
+} // namespace generic
+
+} // namespace LIBC_NAMESPACE
+
+#endif // LLVM_LIBC_SRC_MATH_GENERIC_RANGE_REDUCTION_DOUBLE_COMMON_H
diff --git a/libc/src/math/generic/range_reduction_double_fma.h b/libc/src/math/generic/range_reduction_double_fma.h
new file mode 100644
index 0000000..c136de9
--- /dev/null
+++ b/libc/src/math/generic/range_reduction_double_fma.h
@@ -0,0 +1,495 @@
+//===-- Range reduction for double precision sin/cos/tan w/ FMA -*- C++ -*-===//
+//
+// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
+// See https://llvm.org/LICENSE.txt for license information.
+// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
+//
+//===----------------------------------------------------------------------===//
+
+#ifndef LLVM_LIBC_SRC_MATH_GENERIC_RANGE_REDUCTION_DOUBLE_FMA_H
+#define LLVM_LIBC_SRC_MATH_GENERIC_RANGE_REDUCTION_DOUBLE_FMA_H
+
+#include "src/__support/FPUtil/FPBits.h"
+#include "src/__support/FPUtil/double_double.h"
+#include "src/__support/FPUtil/multiply_add.h"
+#include "src/__support/FPUtil/nearest_integer.h"
+#include "src/__support/common.h"
+
+namespace LIBC_NAMESPACE {
+
+namespace fma {
+
+using LIBC_NAMESPACE::fputil::DoubleDouble;
+
+LIBC_INLINE constexpr int FAST_PASS_EXPONENT = 32;
+
+// 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},
+};
+
+// Lookup table for sin(k * pi / 128) with k = 0, ..., 255.
+// Table is generated with Sollya as follow:
+// > display = hexadecimal;
+// > for k from 0 to 255 do {
+// a = D(sin(k * pi/128)); };
+// b = D(sin(k * pi/128) - a);
+// print("{", b, ",", a, "},");
+// };
+LIBC_INLINE constexpr DoubleDouble SIN_K_PI_OVER_128[256] = {
+ {0, 0},
+ {-0x1.b1d63091a013p-64, 0x1.92155f7a3667ep-6},
+ {-0x1.912bd0d569a9p-61, 0x1.91f65f10dd814p-5},
+ {-0x1.9a088a8bf6b2cp-59, 0x1.2d52092ce19f6p-4},
+ {-0x1.e2718d26ed688p-60, 0x1.917a6bc29b42cp-4},
+ {0x1.a2704729ae56dp-59, 0x1.f564e56a9730ep-4},
+ {0x1.13000a89a11ep-58, 0x1.2c8106e8e613ap-3},
+ {0x1.531ff779ddac6p-57, 0x1.5e214448b3fc6p-3},
+ {-0x1.26d19b9ff8d82p-57, 0x1.8f8b83c69a60bp-3},
+ {-0x1.af1439e521935p-62, 0x1.c0b826a7e4f63p-3},
+ {-0x1.42deef11da2c4p-57, 0x1.f19f97b215f1bp-3},
+ {0x1.824c20ab7aa9ap-56, 0x1.111d262b1f677p-2},
+ {-0x1.5d28da2c4612dp-56, 0x1.294062ed59f06p-2},
+ {0x1.0c97c4afa2518p-56, 0x1.4135c94176601p-2},
+ {-0x1.efdc0d58cf62p-62, 0x1.58f9a75ab1fddp-2},
+ {-0x1.44b19e0864c5dp-56, 0x1.7088530fa459fp-2},
+ {-0x1.72cedd3d5a61p-57, 0x1.87de2a6aea963p-2},
+ {0x1.6da81290bdbabp-57, 0x1.9ef7943a8ed8ap-2},
+ {0x1.5b362cb974183p-57, 0x1.b5d1009e15ccp-2},
+ {0x1.6850e59c37f8fp-58, 0x1.cc66e9931c45ep-2},
+ {0x1.e0d891d3c6841p-58, 0x1.e2b5d3806f63bp-2},
+ {-0x1.2ec1fc1b776b8p-60, 0x1.f8ba4dbf89abap-2},
+ {-0x1.a5a014347406cp-55, 0x1.073879922ffeep-1},
+ {-0x1.ef23b69abe4f1p-55, 0x1.11eb3541b4b23p-1},
+ {0x1.b25dd267f66p-55, 0x1.1c73b39ae68c8p-1},
+ {-0x1.5da743ef3770cp-55, 0x1.26d054cdd12dfp-1},
+ {-0x1.efcc626f74a6fp-57, 0x1.30ff7fce17035p-1},
+ {0x1.e3e25e3954964p-56, 0x1.3affa292050b9p-1},
+ {0x1.8076a2cfdc6b3p-57, 0x1.44cf325091dd6p-1},
+ {0x1.3c293edceb327p-57, 0x1.4e6cabbe3e5e9p-1},
+ {-0x1.75720992bfbb2p-55, 0x1.57d69348cecap-1},
+ {-0x1.251b352ff2a37p-56, 0x1.610b7551d2cdfp-1},
+ {-0x1.bdd3413b26456p-55, 0x1.6a09e667f3bcdp-1},
+ {0x1.0d4ef0f1d915cp-55, 0x1.72d0837efff96p-1},
+ {-0x1.0f537acdf0ad7p-56, 0x1.7b5df226aafafp-1},
+ {-0x1.6f420f8ea3475p-56, 0x1.83b0e0bff976ep-1},
+ {-0x1.2c5e12ed1336dp-55, 0x1.8bc806b151741p-1},
+ {0x1.3d419a920df0bp-55, 0x1.93a22499263fbp-1},
+ {-0x1.30ee286712474p-55, 0x1.9b3e047f38741p-1},
+ {-0x1.128bb015df175p-56, 0x1.a29a7a0462782p-1},
+ {0x1.9f630e8b6dac8p-60, 0x1.a9b66290ea1a3p-1},
+ {-0x1.926da300ffccep-55, 0x1.b090a581502p-1},
+ {-0x1.bc69f324e6d61p-55, 0x1.b728345196e3ep-1},
+ {-0x1.825a732ac700ap-55, 0x1.bd7c0ac6f952ap-1},
+ {-0x1.6e0b1757c8d07p-56, 0x1.c38b2f180bdb1p-1},
+ {-0x1.2fb761e946603p-58, 0x1.c954b213411f5p-1},
+ {-0x1.e7b6bb5ab58aep-58, 0x1.ced7af43cc773p-1},
+ {-0x1.4ef5295d25af2p-55, 0x1.d4134d14dc93ap-1},
+ {0x1.457e610231ac2p-56, 0x1.d906bcf328d46p-1},
+ {0x1.83c37c6107db3p-55, 0x1.ddb13b6ccc23cp-1},
+ {-0x1.014c76c126527p-55, 0x1.e212104f686e5p-1},
+ {-0x1.16b56f2847754p-57, 0x1.e6288ec48e112p-1},
+ {0x1.760b1e2e3f81ep-55, 0x1.e9f4156c62ddap-1},
+ {0x1.e82c791f59cc2p-56, 0x1.ed740e7684963p-1},
+ {0x1.52c7adc6b4989p-56, 0x1.f0a7efb9230d7p-1},
+ {-0x1.d7bafb51f72e6p-56, 0x1.f38f3ac64e589p-1},
+ {0x1.562172a361fd3p-56, 0x1.f6297cff75cbp-1},
+ {0x1.ab256778ffcb6p-56, 0x1.f8764fa714ba9p-1},
+ {-0x1.7a0a8ca13571fp-55, 0x1.fa7557f08a517p-1},
+ {0x1.1ec8668ecaceep-55, 0x1.fc26470e19fd3p-1},
+ {-0x1.87df6378811c7p-55, 0x1.fd88da3d12526p-1},
+ {0x1.521ecd0c67e35p-57, 0x1.fe9cdad01883ap-1},
+ {-0x1.c57bc2e24aa15p-57, 0x1.ff621e3796d7ep-1},
+ {-0x1.1354d4556e4cbp-55, 0x1.ffd886084cd0dp-1},
+ {0, 1},
+ {-0x1.1354d4556e4cbp-55, 0x1.ffd886084cd0dp-1},
+ {-0x1.c57bc2e24aa15p-57, 0x1.ff621e3796d7ep-1},
+ {0x1.521ecd0c67e35p-57, 0x1.fe9cdad01883ap-1},
+ {-0x1.87df6378811c7p-55, 0x1.fd88da3d12526p-1},
+ {0x1.1ec8668ecaceep-55, 0x1.fc26470e19fd3p-1},
+ {-0x1.7a0a8ca13571fp-55, 0x1.fa7557f08a517p-1},
+ {0x1.ab256778ffcb6p-56, 0x1.f8764fa714ba9p-1},
+ {0x1.562172a361fd3p-56, 0x1.f6297cff75cbp-1},
+ {-0x1.d7bafb51f72e6p-56, 0x1.f38f3ac64e589p-1},
+ {0x1.52c7adc6b4989p-56, 0x1.f0a7efb9230d7p-1},
+ {0x1.e82c791f59cc2p-56, 0x1.ed740e7684963p-1},
+ {0x1.760b1e2e3f81ep-55, 0x1.e9f4156c62ddap-1},
+ {-0x1.16b56f2847754p-57, 0x1.e6288ec48e112p-1},
+ {-0x1.014c76c126527p-55, 0x1.e212104f686e5p-1},
+ {0x1.83c37c6107db3p-55, 0x1.ddb13b6ccc23cp-1},
+ {0x1.457e610231ac2p-56, 0x1.d906bcf328d46p-1},
+ {-0x1.4ef5295d25af2p-55, 0x1.d4134d14dc93ap-1},
+ {-0x1.e7b6bb5ab58aep-58, 0x1.ced7af43cc773p-1},
+ {-0x1.2fb761e946603p-58, 0x1.c954b213411f5p-1},
+ {-0x1.6e0b1757c8d07p-56, 0x1.c38b2f180bdb1p-1},
+ {-0x1.825a732ac700ap-55, 0x1.bd7c0ac6f952ap-1},
+ {-0x1.bc69f324e6d61p-55, 0x1.b728345196e3ep-1},
+ {-0x1.926da300ffccep-55, 0x1.b090a581502p-1},
+ {0x1.9f630e8b6dac8p-60, 0x1.a9b66290ea1a3p-1},
+ {-0x1.128bb015df175p-56, 0x1.a29a7a0462782p-1},
+ {-0x1.30ee286712474p-55, 0x1.9b3e047f38741p-1},
+ {0x1.3d419a920df0bp-55, 0x1.93a22499263fbp-1},
+ {-0x1.2c5e12ed1336dp-55, 0x1.8bc806b151741p-1},
+ {-0x1.6f420f8ea3475p-56, 0x1.83b0e0bff976ep-1},
+ {-0x1.0f537acdf0ad7p-56, 0x1.7b5df226aafafp-1},
+ {0x1.0d4ef0f1d915cp-55, 0x1.72d0837efff96p-1},
+ {-0x1.bdd3413b26456p-55, 0x1.6a09e667f3bcdp-1},
+ {-0x1.251b352ff2a37p-56, 0x1.610b7551d2cdfp-1},
+ {-0x1.75720992bfbb2p-55, 0x1.57d69348cecap-1},
+ {0x1.3c293edceb327p-57, 0x1.4e6cabbe3e5e9p-1},
+ {0x1.8076a2cfdc6b3p-57, 0x1.44cf325091dd6p-1},
+ {0x1.e3e25e3954964p-56, 0x1.3affa292050b9p-1},
+ {-0x1.efcc626f74a6fp-57, 0x1.30ff7fce17035p-1},
+ {-0x1.5da743ef3770cp-55, 0x1.26d054cdd12dfp-1},
+ {0x1.b25dd267f66p-55, 0x1.1c73b39ae68c8p-1},
+ {-0x1.ef23b69abe4f1p-55, 0x1.11eb3541b4b23p-1},
+ {-0x1.a5a014347406cp-55, 0x1.073879922ffeep-1},
+ {-0x1.2ec1fc1b776b8p-60, 0x1.f8ba4dbf89abap-2},
+ {0x1.e0d891d3c6841p-58, 0x1.e2b5d3806f63bp-2},
+ {0x1.6850e59c37f8fp-58, 0x1.cc66e9931c45ep-2},
+ {0x1.5b362cb974183p-57, 0x1.b5d1009e15ccp-2},
+ {0x1.6da81290bdbabp-57, 0x1.9ef7943a8ed8ap-2},
+ {-0x1.72cedd3d5a61p-57, 0x1.87de2a6aea963p-2},
+ {-0x1.44b19e0864c5dp-56, 0x1.7088530fa459fp-2},
+ {-0x1.efdc0d58cf62p-62, 0x1.58f9a75ab1fddp-2},
+ {0x1.0c97c4afa2518p-56, 0x1.4135c94176601p-2},
+ {-0x1.5d28da2c4612dp-56, 0x1.294062ed59f06p-2},
+ {0x1.824c20ab7aa9ap-56, 0x1.111d262b1f677p-2},
+ {-0x1.42deef11da2c4p-57, 0x1.f19f97b215f1bp-3},
+ {-0x1.af1439e521935p-62, 0x1.c0b826a7e4f63p-3},
+ {-0x1.26d19b9ff8d82p-57, 0x1.8f8b83c69a60bp-3},
+ {0x1.531ff779ddac6p-57, 0x1.5e214448b3fc6p-3},
+ {0x1.13000a89a11ep-58, 0x1.2c8106e8e613ap-3},
+ {0x1.a2704729ae56dp-59, 0x1.f564e56a9730ep-4},
+ {-0x1.e2718d26ed688p-60, 0x1.917a6bc29b42cp-4},
+ {-0x1.9a088a8bf6b2cp-59, 0x1.2d52092ce19f6p-4},
+ {-0x1.912bd0d569a9p-61, 0x1.91f65f10dd814p-5},
+ {-0x1.b1d63091a013p-64, 0x1.92155f7a3667ep-6},
+ {0, 0},
+ {0x1.b1d63091a013p-64, -0x1.92155f7a3667ep-6},
+ {0x1.912bd0d569a9p-61, -0x1.91f65f10dd814p-5},
+ {0x1.9a088a8bf6b2cp-59, -0x1.2d52092ce19f6p-4},
+ {0x1.e2718d26ed688p-60, -0x1.917a6bc29b42cp-4},
+ {-0x1.a2704729ae56dp-59, -0x1.f564e56a9730ep-4},
+ {-0x1.13000a89a11ep-58, -0x1.2c8106e8e613ap-3},
+ {-0x1.531ff779ddac6p-57, -0x1.5e214448b3fc6p-3},
+ {0x1.26d19b9ff8d82p-57, -0x1.8f8b83c69a60bp-3},
+ {0x1.af1439e521935p-62, -0x1.c0b826a7e4f63p-3},
+ {0x1.42deef11da2c4p-57, -0x1.f19f97b215f1bp-3},
+ {-0x1.824c20ab7aa9ap-56, -0x1.111d262b1f677p-2},
+ {0x1.5d28da2c4612dp-56, -0x1.294062ed59f06p-2},
+ {-0x1.0c97c4afa2518p-56, -0x1.4135c94176601p-2},
+ {0x1.efdc0d58cf62p-62, -0x1.58f9a75ab1fddp-2},
+ {0x1.44b19e0864c5dp-56, -0x1.7088530fa459fp-2},
+ {0x1.72cedd3d5a61p-57, -0x1.87de2a6aea963p-2},
+ {-0x1.6da81290bdbabp-57, -0x1.9ef7943a8ed8ap-2},
+ {-0x1.5b362cb974183p-57, -0x1.b5d1009e15ccp-2},
+ {-0x1.6850e59c37f8fp-58, -0x1.cc66e9931c45ep-2},
+ {-0x1.e0d891d3c6841p-58, -0x1.e2b5d3806f63bp-2},
+ {0x1.2ec1fc1b776b8p-60, -0x1.f8ba4dbf89abap-2},
+ {0x1.a5a014347406cp-55, -0x1.073879922ffeep-1},
+ {0x1.ef23b69abe4f1p-55, -0x1.11eb3541b4b23p-1},
+ {-0x1.b25dd267f66p-55, -0x1.1c73b39ae68c8p-1},
+ {0x1.5da743ef3770cp-55, -0x1.26d054cdd12dfp-1},
+ {0x1.efcc626f74a6fp-57, -0x1.30ff7fce17035p-1},
+ {-0x1.e3e25e3954964p-56, -0x1.3affa292050b9p-1},
+ {-0x1.8076a2cfdc6b3p-57, -0x1.44cf325091dd6p-1},
+ {-0x1.3c293edceb327p-57, -0x1.4e6cabbe3e5e9p-1},
+ {0x1.75720992bfbb2p-55, -0x1.57d69348cecap-1},
+ {0x1.251b352ff2a37p-56, -0x1.610b7551d2cdfp-1},
+ {0x1.bdd3413b26456p-55, -0x1.6a09e667f3bcdp-1},
+ {-0x1.0d4ef0f1d915cp-55, -0x1.72d0837efff96p-1},
+ {0x1.0f537acdf0ad7p-56, -0x1.7b5df226aafafp-1},
+ {0x1.6f420f8ea3475p-56, -0x1.83b0e0bff976ep-1},
+ {0x1.2c5e12ed1336dp-55, -0x1.8bc806b151741p-1},
+ {-0x1.3d419a920df0bp-55, -0x1.93a22499263fbp-1},
+ {0x1.30ee286712474p-55, -0x1.9b3e047f38741p-1},
+ {0x1.128bb015df175p-56, -0x1.a29a7a0462782p-1},
+ {-0x1.9f630e8b6dac8p-60, -0x1.a9b66290ea1a3p-1},
+ {0x1.926da300ffccep-55, -0x1.b090a581502p-1},
+ {0x1.bc69f324e6d61p-55, -0x1.b728345196e3ep-1},
+ {0x1.825a732ac700ap-55, -0x1.bd7c0ac6f952ap-1},
+ {0x1.6e0b1757c8d07p-56, -0x1.c38b2f180bdb1p-1},
+ {0x1.2fb761e946603p-58, -0x1.c954b213411f5p-1},
+ {0x1.e7b6bb5ab58aep-58, -0x1.ced7af43cc773p-1},
+ {0x1.4ef5295d25af2p-55, -0x1.d4134d14dc93ap-1},
+ {-0x1.457e610231ac2p-56, -0x1.d906bcf328d46p-1},
+ {-0x1.83c37c6107db3p-55, -0x1.ddb13b6ccc23cp-1},
+ {0x1.014c76c126527p-55, -0x1.e212104f686e5p-1},
+ {0x1.16b56f2847754p-57, -0x1.e6288ec48e112p-1},
+ {-0x1.760b1e2e3f81ep-55, -0x1.e9f4156c62ddap-1},
+ {-0x1.e82c791f59cc2p-56, -0x1.ed740e7684963p-1},
+ {-0x1.52c7adc6b4989p-56, -0x1.f0a7efb9230d7p-1},
+ {0x1.d7bafb51f72e6p-56, -0x1.f38f3ac64e589p-1},
+ {-0x1.562172a361fd3p-56, -0x1.f6297cff75cbp-1},
+ {-0x1.ab256778ffcb6p-56, -0x1.f8764fa714ba9p-1},
+ {0x1.7a0a8ca13571fp-55, -0x1.fa7557f08a517p-1},
+ {-0x1.1ec8668ecaceep-55, -0x1.fc26470e19fd3p-1},
+ {0x1.87df6378811c7p-55, -0x1.fd88da3d12526p-1},
+ {-0x1.521ecd0c67e35p-57, -0x1.fe9cdad01883ap-1},
+ {0x1.c57bc2e24aa15p-57, -0x1.ff621e3796d7ep-1},
+ {0x1.1354d4556e4cbp-55, -0x1.ffd886084cd0dp-1},
+ {0, -1},
+ {0x1.1354d4556e4cbp-55, -0x1.ffd886084cd0dp-1},
+ {0x1.c57bc2e24aa15p-57, -0x1.ff621e3796d7ep-1},
+ {-0x1.521ecd0c67e35p-57, -0x1.fe9cdad01883ap-1},
+ {0x1.87df6378811c7p-55, -0x1.fd88da3d12526p-1},
+ {-0x1.1ec8668ecaceep-55, -0x1.fc26470e19fd3p-1},
+ {0x1.7a0a8ca13571fp-55, -0x1.fa7557f08a517p-1},
+ {-0x1.ab256778ffcb6p-56, -0x1.f8764fa714ba9p-1},
+ {-0x1.562172a361fd3p-56, -0x1.f6297cff75cbp-1},
+ {0x1.d7bafb51f72e6p-56, -0x1.f38f3ac64e589p-1},
+ {-0x1.52c7adc6b4989p-56, -0x1.f0a7efb9230d7p-1},
+ {-0x1.e82c791f59cc2p-56, -0x1.ed740e7684963p-1},
+ {-0x1.760b1e2e3f81ep-55, -0x1.e9f4156c62ddap-1},
+ {0x1.16b56f2847754p-57, -0x1.e6288ec48e112p-1},
+ {0x1.014c76c126527p-55, -0x1.e212104f686e5p-1},
+ {-0x1.83c37c6107db3p-55, -0x1.ddb13b6ccc23cp-1},
+ {-0x1.457e610231ac2p-56, -0x1.d906bcf328d46p-1},
+ {0x1.4ef5295d25af2p-55, -0x1.d4134d14dc93ap-1},
+ {0x1.e7b6bb5ab58aep-58, -0x1.ced7af43cc773p-1},
+ {0x1.2fb761e946603p-58, -0x1.c954b213411f5p-1},
+ {0x1.6e0b1757c8d07p-56, -0x1.c38b2f180bdb1p-1},
+ {0x1.825a732ac700ap-55, -0x1.bd7c0ac6f952ap-1},
+ {0x1.bc69f324e6d61p-55, -0x1.b728345196e3ep-1},
+ {0x1.926da300ffccep-55, -0x1.b090a581502p-1},
+ {-0x1.9f630e8b6dac8p-60, -0x1.a9b66290ea1a3p-1},
+ {0x1.128bb015df175p-56, -0x1.a29a7a0462782p-1},
+ {0x1.30ee286712474p-55, -0x1.9b3e047f38741p-1},
+ {-0x1.3d419a920df0bp-55, -0x1.93a22499263fbp-1},
+ {0x1.2c5e12ed1336dp-55, -0x1.8bc806b151741p-1},
+ {0x1.6f420f8ea3475p-56, -0x1.83b0e0bff976ep-1},
+ {0x1.0f537acdf0ad7p-56, -0x1.7b5df226aafafp-1},
+ {-0x1.0d4ef0f1d915cp-55, -0x1.72d0837efff96p-1},
+ {0x1.bdd3413b26456p-55, -0x1.6a09e667f3bcdp-1},
+ {0x1.251b352ff2a37p-56, -0x1.610b7551d2cdfp-1},
+ {0x1.75720992bfbb2p-55, -0x1.57d69348cecap-1},
+ {-0x1.3c293edceb327p-57, -0x1.4e6cabbe3e5e9p-1},
+ {-0x1.8076a2cfdc6b3p-57, -0x1.44cf325091dd6p-1},
+ {-0x1.e3e25e3954964p-56, -0x1.3affa292050b9p-1},
+ {0x1.efcc626f74a6fp-57, -0x1.30ff7fce17035p-1},
+ {0x1.5da743ef3770cp-55, -0x1.26d054cdd12dfp-1},
+ {-0x1.b25dd267f66p-55, -0x1.1c73b39ae68c8p-1},
+ {0x1.ef23b69abe4f1p-55, -0x1.11eb3541b4b23p-1},
+ {0x1.a5a014347406cp-55, -0x1.073879922ffeep-1},
+ {0x1.2ec1fc1b776b8p-60, -0x1.f8ba4dbf89abap-2},
+ {-0x1.e0d891d3c6841p-58, -0x1.e2b5d3806f63bp-2},
+ {-0x1.6850e59c37f8fp-58, -0x1.cc66e9931c45ep-2},
+ {-0x1.5b362cb974183p-57, -0x1.b5d1009e15ccp-2},
+ {-0x1.6da81290bdbabp-57, -0x1.9ef7943a8ed8ap-2},
+ {0x1.72cedd3d5a61p-57, -0x1.87de2a6aea963p-2},
+ {0x1.44b19e0864c5dp-56, -0x1.7088530fa459fp-2},
+ {0x1.efdc0d58cf62p-62, -0x1.58f9a75ab1fddp-2},
+ {-0x1.0c97c4afa2518p-56, -0x1.4135c94176601p-2},
+ {0x1.5d28da2c4612dp-56, -0x1.294062ed59f06p-2},
+ {-0x1.824c20ab7aa9ap-56, -0x1.111d262b1f677p-2},
+ {0x1.42deef11da2c4p-57, -0x1.f19f97b215f1bp-3},
+ {0x1.af1439e521935p-62, -0x1.c0b826a7e4f63p-3},
+ {0x1.26d19b9ff8d82p-57, -0x1.8f8b83c69a60bp-3},
+ {-0x1.531ff779ddac6p-57, -0x1.5e214448b3fc6p-3},
+ {-0x1.13000a89a11ep-58, -0x1.2c8106e8e613ap-3},
+ {-0x1.a2704729ae56dp-59, -0x1.f564e56a9730ep-4},
+ {0x1.e2718d26ed688p-60, -0x1.917a6bc29b42cp-4},
+ {0x1.9a088a8bf6b2cp-59, -0x1.2d52092ce19f6p-4},
+ {0x1.912bd0d569a9p-61, -0x1.91f65f10dd814p-5},
+ {0x1.b1d63091a013p-64, -0x1.92155f7a3667ep-6},
+};
+
+// 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
+
+#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
new file mode 100644
index 0000000..b9d34d6
--- /dev/null
+++ b/libc/src/math/generic/range_reduction_double_nofma.h
@@ -0,0 +1,493 @@
+//===-- Range reduction for double precision sin/cos/tan --------*- C++ -*-===//
+//
+// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
+// See https://llvm.org/LICENSE.txt for license information.
+// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
+//
+//===----------------------------------------------------------------------===//
+
+#ifndef LLVM_LIBC_SRC_MATH_GENERIC_RANGE_REDUCTION_DOUBLE_NOFMA_H
+#define LLVM_LIBC_SRC_MATH_GENERIC_RANGE_REDUCTION_DOUBLE_NOFMA_H
+
+#include "src/__support/FPUtil/FPBits.h"
+#include "src/__support/FPUtil/double_double.h"
+#include "src/__support/FPUtil/multiply_add.h"
+#include "src/__support/FPUtil/nearest_integer.h"
+#include "src/__support/common.h"
+
+namespace LIBC_NAMESPACE {
+
+namespace nofma {
+
+using fputil::DoubleDouble;
+
+LIBC_INLINE constexpr int FAST_PASS_EXPONENT = 23;
+
+// 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},
+};
+
+// Lookup table for sin(k * pi / 128) with k = 0, ..., 255.
+// Table is generated with Sollya as follow:
+// > display = hexadecimal;
+// > for k from 0 to 255 do {
+// a = round(sin(k * pi/128), 51, RN);
+// b = round(sin(k * pi/128) - a, D, RN);
+// print("{", b, ",", a, "},");
+// };
+LIBC_INLINE constexpr DoubleDouble SIN_K_PI_OVER_128[256] = {
+ {0, 0},
+ {0x1.f938a73db97fbp-58, 0x1.92155f7a3667cp-6},
+ {-0x1.912bd0d569a9p-61, 0x1.91f65f10dd814p-5},
+ {0x1.ccbeeeae8129ap-56, 0x1.2d52092ce19f4p-4},
+ {-0x1.e2718d26ed688p-60, 0x1.917a6bc29b42cp-4},
+ {-0x1.cbb1f71aca352p-56, 0x1.f564e56a9731p-4},
+ {-0x1.dd9ffeaecbdc4p-55, 0x1.2c8106e8e613cp-3},
+ {-0x1.ab3802218894fp-55, 0x1.5e214448b3fc8p-3},
+ {-0x1.49b466e7fe36p-55, 0x1.8f8b83c69a60cp-3},
+ {-0x1.035e2873ca432p-55, 0x1.c0b826a7e4f64p-3},
+ {-0x1.50b7bbc4768b1p-55, 0x1.f19f97b215f1cp-3},
+ {-0x1.3ed9efaa42ab3p-55, 0x1.111d262b1f678p-2},
+ {0x1.a8b5c974ee7b5p-54, 0x1.294062ed59f04p-2},
+ {0x1.4325f12be8946p-54, 0x1.4135c941766p-2},
+ {0x1.fc2047e54e614p-55, 0x1.58f9a75ab1fdcp-2},
+ {-0x1.512c678219317p-54, 0x1.7088530fa45ap-2},
+ {-0x1.2e59dba7ab4c2p-54, 0x1.87de2a6aea964p-2},
+ {-0x1.d24afdade848bp-54, 0x1.9ef7943a8ed8cp-2},
+ {0x1.5b362cb974183p-57, 0x1.b5d1009e15ccp-2},
+ {-0x1.e97af1a63c807p-54, 0x1.cc66e9931c46p-2},
+ {-0x1.c3e4edc5872f8p-55, 0x1.e2b5d3806f63cp-2},
+ {0x1.fb44f80f92225p-54, 0x1.f8ba4dbf89ab8p-2},
+ {0x1.9697faf2e2fe5p-53, 0x1.073879922ffecp-1},
+ {-0x1.7bc8eda6af93cp-53, 0x1.11eb3541b4b24p-1},
+ {0x1.b25dd267f66p-55, 0x1.1c73b39ae68c8p-1},
+ {-0x1.5769d0fbcddc3p-53, 0x1.26d054cdd12ep-1},
+ {0x1.c20673b2116b2p-54, 0x1.30ff7fce17034p-1},
+ {0x1.3c7c4bc72a92cp-53, 0x1.3affa292050b8p-1},
+ {-0x1.e7f895d302395p-53, 0x1.44cf325091dd8p-1},
+ {0x1.13c293edceb32p-53, 0x1.4e6cabbe3e5e8p-1},
+ {-0x1.75720992bfbb2p-55, 0x1.57d69348cecap-1},
+ {-0x1.24a366a5fe547p-53, 0x1.610b7551d2cep-1},
+ {0x1.21165f626cdd5p-54, 0x1.6a09e667f3bccp-1},
+ {-0x1.bcac43c389ba9p-53, 0x1.72d0837efff98p-1},
+ {-0x1.21ea6f59be15bp-53, 0x1.7b5df226aafbp-1},
+ {0x1.d217be0e2b971p-53, 0x1.83b0e0bff976cp-1},
+ {0x1.69d0f6897664ap-54, 0x1.8bc806b15174p-1},
+ {-0x1.615f32b6f907ap-54, 0x1.93a22499263fcp-1},
+ {0x1.6788ebcc76dc6p-54, 0x1.9b3e047f3874p-1},
+ {0x1.ddae89fd441d1p-53, 0x1.a29a7a046278p-1},
+ {-0x1.f98273c5d2495p-54, 0x1.a9b66290ea1a4p-1},
+ {-0x1.926da300ffccep-55, 0x1.b090a581502p-1},
+ {0x1.90e58336c64a8p-53, 0x1.b728345196e3cp-1},
+ {0x1.9f6963354e3fep-53, 0x1.bd7c0ac6f9528p-1},
+ {0x1.a47d3a2a0dcbep-54, 0x1.c38b2f180bdbp-1},
+ {0x1.ed0489e16b9ap-54, 0x1.c954b213411f4p-1},
+ {-0x1.0f3db5dad5ac5p-53, 0x1.ced7af43cc774p-1},
+ {0x1.ac42b5a8b6943p-53, 0x1.d4134d14dc938p-1},
+ {-0x1.d75033dfb9ca8p-53, 0x1.d906bcf328d48p-1},
+ {0x1.83c37c6107db3p-55, 0x1.ddb13b6ccc23cp-1},
+ {0x1.7f59c49f6cd6dp-54, 0x1.e212104f686e4p-1},
+ {0x1.ee94a90d7b88bp-53, 0x1.e6288ec48e11p-1},
+ {-0x1.a27d3874701f9p-53, 0x1.e9f4156c62ddcp-1},
+ {-0x1.85f4e1b8298dp-54, 0x1.ed740e7684964p-1},
+ {-0x1.ab4e148e52d9ep-54, 0x1.f0a7efb9230d8p-1},
+ {0x1.8a11412b82346p-54, 0x1.f38f3ac64e588p-1},
+ {0x1.562172a361fd3p-56, 0x1.f6297cff75cbp-1},
+ {0x1.3564acef1ff97p-53, 0x1.f8764fa714ba8p-1},
+ {-0x1.5e82a3284d5c8p-53, 0x1.fa7557f08a518p-1},
+ {-0x1.709bccb89a989p-54, 0x1.fc26470e19fd4p-1},
+ {0x1.9e082721dfb8ep-53, 0x1.fd88da3d12524p-1},
+ {-0x1.eade132f3981dp-53, 0x1.fe9cdad01883cp-1},
+ {0x1.e3a843d1db55fp-53, 0x1.ff621e3796d7cp-1},
+ {0x1.765595d548d9ap-54, 0x1.ffd886084cd0cp-1},
+ {0, 1},
+ {0x1.765595d548d9ap-54, 0x1.ffd886084cd0cp-1},
+ {0x1.e3a843d1db55fp-53, 0x1.ff621e3796d7cp-1},
+ {-0x1.eade132f3981dp-53, 0x1.fe9cdad01883cp-1},
+ {0x1.9e082721dfb8ep-53, 0x1.fd88da3d12524p-1},
+ {-0x1.709bccb89a989p-54, 0x1.fc26470e19fd4p-1},
+ {-0x1.5e82a3284d5c8p-53, 0x1.fa7557f08a518p-1},
+ {0x1.3564acef1ff97p-53, 0x1.f8764fa714ba8p-1},
+ {0x1.562172a361fd3p-56, 0x1.f6297cff75cbp-1},
+ {0x1.8a11412b82346p-54, 0x1.f38f3ac64e588p-1},
+ {-0x1.ab4e148e52d9ep-54, 0x1.f0a7efb9230d8p-1},
+ {-0x1.85f4e1b8298dp-54, 0x1.ed740e7684964p-1},
+ {-0x1.a27d3874701f9p-53, 0x1.e9f4156c62ddcp-1},
+ {0x1.ee94a90d7b88bp-53, 0x1.e6288ec48e11p-1},
+ {0x1.7f59c49f6cd6dp-54, 0x1.e212104f686e4p-1},
+ {0x1.83c37c6107db3p-55, 0x1.ddb13b6ccc23cp-1},
+ {-0x1.d75033dfb9ca8p-53, 0x1.d906bcf328d48p-1},
+ {0x1.ac42b5a8b6943p-53, 0x1.d4134d14dc938p-1},
+ {-0x1.0f3db5dad5ac5p-53, 0x1.ced7af43cc774p-1},
+ {0x1.ed0489e16b9ap-54, 0x1.c954b213411f4p-1},
+ {0x1.a47d3a2a0dcbep-54, 0x1.c38b2f180bdbp-1},
+ {0x1.9f6963354e3fep-53, 0x1.bd7c0ac6f9528p-1},
+ {0x1.90e58336c64a8p-53, 0x1.b728345196e3cp-1},
+ {-0x1.926da300ffccep-55, 0x1.b090a581502p-1},
+ {-0x1.f98273c5d2495p-54, 0x1.a9b66290ea1a4p-1},
+ {0x1.ddae89fd441d1p-53, 0x1.a29a7a046278p-1},
+ {0x1.6788ebcc76dc6p-54, 0x1.9b3e047f3874p-1},
+ {-0x1.615f32b6f907ap-54, 0x1.93a22499263fcp-1},
+ {0x1.69d0f6897664ap-54, 0x1.8bc806b15174p-1},
+ {0x1.d217be0e2b971p-53, 0x1.83b0e0bff976cp-1},
+ {-0x1.21ea6f59be15bp-53, 0x1.7b5df226aafbp-1},
+ {-0x1.bcac43c389ba9p-53, 0x1.72d0837efff98p-1},
+ {0x1.21165f626cdd5p-54, 0x1.6a09e667f3bccp-1},
+ {-0x1.24a366a5fe547p-53, 0x1.610b7551d2cep-1},
+ {-0x1.75720992bfbb2p-55, 0x1.57d69348cecap-1},
+ {0x1.13c293edceb32p-53, 0x1.4e6cabbe3e5e8p-1},
+ {-0x1.e7f895d302395p-53, 0x1.44cf325091dd8p-1},
+ {0x1.3c7c4bc72a92cp-53, 0x1.3affa292050b8p-1},
+ {0x1.c20673b2116b2p-54, 0x1.30ff7fce17034p-1},
+ {-0x1.5769d0fbcddc3p-53, 0x1.26d054cdd12ep-1},
+ {0x1.b25dd267f66p-55, 0x1.1c73b39ae68c8p-1},
+ {-0x1.7bc8eda6af93cp-53, 0x1.11eb3541b4b24p-1},
+ {0x1.9697faf2e2fe5p-53, 0x1.073879922ffecp-1},
+ {0x1.fb44f80f92225p-54, 0x1.f8ba4dbf89ab8p-2},
+ {-0x1.c3e4edc5872f8p-55, 0x1.e2b5d3806f63cp-2},
+ {-0x1.e97af1a63c807p-54, 0x1.cc66e9931c46p-2},
+ {0x1.5b362cb974183p-57, 0x1.b5d1009e15ccp-2},
+ {-0x1.d24afdade848bp-54, 0x1.9ef7943a8ed8cp-2},
+ {-0x1.2e59dba7ab4c2p-54, 0x1.87de2a6aea964p-2},
+ {-0x1.512c678219317p-54, 0x1.7088530fa45ap-2},
+ {0x1.fc2047e54e614p-55, 0x1.58f9a75ab1fdcp-2},
+ {0x1.4325f12be8946p-54, 0x1.4135c941766p-2},
+ {0x1.a8b5c974ee7b5p-54, 0x1.294062ed59f04p-2},
+ {-0x1.3ed9efaa42ab3p-55, 0x1.111d262b1f678p-2},
+ {-0x1.50b7bbc4768b1p-55, 0x1.f19f97b215f1cp-3},
+ {-0x1.035e2873ca432p-55, 0x1.c0b826a7e4f64p-3},
+ {-0x1.49b466e7fe36p-55, 0x1.8f8b83c69a60cp-3},
+ {-0x1.ab3802218894fp-55, 0x1.5e214448b3fc8p-3},
+ {-0x1.dd9ffeaecbdc4p-55, 0x1.2c8106e8e613cp-3},
+ {-0x1.cbb1f71aca352p-56, 0x1.f564e56a9731p-4},
+ {-0x1.e2718d26ed688p-60, 0x1.917a6bc29b42cp-4},
+ {0x1.ccbeeeae8129ap-56, 0x1.2d52092ce19f4p-4},
+ {-0x1.912bd0d569a9p-61, 0x1.91f65f10dd814p-5},
+ {0x1.f938a73db97fbp-58, 0x1.92155f7a3667cp-6},
+ {0, 0},
+ {-0x1.f938a73db97fbp-58, -0x1.92155f7a3667cp-6},
+ {0x1.912bd0d569a9p-61, -0x1.91f65f10dd814p-5},
+ {-0x1.ccbeeeae8129ap-56, -0x1.2d52092ce19f4p-4},
+ {0x1.e2718d26ed688p-60, -0x1.917a6bc29b42cp-4},
+ {0x1.cbb1f71aca352p-56, -0x1.f564e56a9731p-4},
+ {0x1.dd9ffeaecbdc4p-55, -0x1.2c8106e8e613cp-3},
+ {0x1.ab3802218894fp-55, -0x1.5e214448b3fc8p-3},
+ {0x1.49b466e7fe36p-55, -0x1.8f8b83c69a60cp-3},
+ {0x1.035e2873ca432p-55, -0x1.c0b826a7e4f64p-3},
+ {0x1.50b7bbc4768b1p-55, -0x1.f19f97b215f1cp-3},
+ {0x1.3ed9efaa42ab3p-55, -0x1.111d262b1f678p-2},
+ {-0x1.a8b5c974ee7b5p-54, -0x1.294062ed59f04p-2},
+ {-0x1.4325f12be8946p-54, -0x1.4135c941766p-2},
+ {-0x1.fc2047e54e614p-55, -0x1.58f9a75ab1fdcp-2},
+ {0x1.512c678219317p-54, -0x1.7088530fa45ap-2},
+ {0x1.2e59dba7ab4c2p-54, -0x1.87de2a6aea964p-2},
+ {0x1.d24afdade848bp-54, -0x1.9ef7943a8ed8cp-2},
+ {-0x1.5b362cb974183p-57, -0x1.b5d1009e15ccp-2},
+ {0x1.e97af1a63c807p-54, -0x1.cc66e9931c46p-2},
+ {0x1.c3e4edc5872f8p-55, -0x1.e2b5d3806f63cp-2},
+ {-0x1.fb44f80f92225p-54, -0x1.f8ba4dbf89ab8p-2},
+ {-0x1.9697faf2e2fe5p-53, -0x1.073879922ffecp-1},
+ {0x1.7bc8eda6af93cp-53, -0x1.11eb3541b4b24p-1},
+ {-0x1.b25dd267f66p-55, -0x1.1c73b39ae68c8p-1},
+ {0x1.5769d0fbcddc3p-53, -0x1.26d054cdd12ep-1},
+ {-0x1.c20673b2116b2p-54, -0x1.30ff7fce17034p-1},
+ {-0x1.3c7c4bc72a92cp-53, -0x1.3affa292050b8p-1},
+ {0x1.e7f895d302395p-53, -0x1.44cf325091dd8p-1},
+ {-0x1.13c293edceb32p-53, -0x1.4e6cabbe3e5e8p-1},
+ {0x1.75720992bfbb2p-55, -0x1.57d69348cecap-1},
+ {0x1.24a366a5fe547p-53, -0x1.610b7551d2cep-1},
+ {-0x1.21165f626cdd5p-54, -0x1.6a09e667f3bccp-1},
+ {0x1.bcac43c389ba9p-53, -0x1.72d0837efff98p-1},
+ {0x1.21ea6f59be15bp-53, -0x1.7b5df226aafbp-1},
+ {-0x1.d217be0e2b971p-53, -0x1.83b0e0bff976cp-1},
+ {-0x1.69d0f6897664ap-54, -0x1.8bc806b15174p-1},
+ {0x1.615f32b6f907ap-54, -0x1.93a22499263fcp-1},
+ {-0x1.6788ebcc76dc6p-54, -0x1.9b3e047f3874p-1},
+ {-0x1.ddae89fd441d1p-53, -0x1.a29a7a046278p-1},
+ {0x1.f98273c5d2495p-54, -0x1.a9b66290ea1a4p-1},
+ {0x1.926da300ffccep-55, -0x1.b090a581502p-1},
+ {-0x1.90e58336c64a8p-53, -0x1.b728345196e3cp-1},
+ {-0x1.9f6963354e3fep-53, -0x1.bd7c0ac6f9528p-1},
+ {-0x1.a47d3a2a0dcbep-54, -0x1.c38b2f180bdbp-1},
+ {-0x1.ed0489e16b9ap-54, -0x1.c954b213411f4p-1},
+ {0x1.0f3db5dad5ac5p-53, -0x1.ced7af43cc774p-1},
+ {-0x1.ac42b5a8b6943p-53, -0x1.d4134d14dc938p-1},
+ {0x1.d75033dfb9ca8p-53, -0x1.d906bcf328d48p-1},
+ {-0x1.83c37c6107db3p-55, -0x1.ddb13b6ccc23cp-1},
+ {-0x1.7f59c49f6cd6dp-54, -0x1.e212104f686e4p-1},
+ {-0x1.ee94a90d7b88bp-53, -0x1.e6288ec48e11p-1},
+ {0x1.a27d3874701f9p-53, -0x1.e9f4156c62ddcp-1},
+ {0x1.85f4e1b8298dp-54, -0x1.ed740e7684964p-1},
+ {0x1.ab4e148e52d9ep-54, -0x1.f0a7efb9230d8p-1},
+ {-0x1.8a11412b82346p-54, -0x1.f38f3ac64e588p-1},
+ {-0x1.562172a361fd3p-56, -0x1.f6297cff75cbp-1},
+ {-0x1.3564acef1ff97p-53, -0x1.f8764fa714ba8p-1},
+ {0x1.5e82a3284d5c8p-53, -0x1.fa7557f08a518p-1},
+ {0x1.709bccb89a989p-54, -0x1.fc26470e19fd4p-1},
+ {-0x1.9e082721dfb8ep-53, -0x1.fd88da3d12524p-1},
+ {0x1.eade132f3981dp-53, -0x1.fe9cdad01883cp-1},
+ {-0x1.e3a843d1db55fp-53, -0x1.ff621e3796d7cp-1},
+ {-0x1.765595d548d9ap-54, -0x1.ffd886084cd0cp-1},
+ {0, -1},
+ {-0x1.765595d548d9ap-54, -0x1.ffd886084cd0cp-1},
+ {-0x1.e3a843d1db55fp-53, -0x1.ff621e3796d7cp-1},
+ {0x1.eade132f3981dp-53, -0x1.fe9cdad01883cp-1},
+ {-0x1.9e082721dfb8ep-53, -0x1.fd88da3d12524p-1},
+ {0x1.709bccb89a989p-54, -0x1.fc26470e19fd4p-1},
+ {0x1.5e82a3284d5c8p-53, -0x1.fa7557f08a518p-1},
+ {-0x1.3564acef1ff97p-53, -0x1.f8764fa714ba8p-1},
+ {-0x1.562172a361fd3p-56, -0x1.f6297cff75cbp-1},
+ {-0x1.8a11412b82346p-54, -0x1.f38f3ac64e588p-1},
+ {0x1.ab4e148e52d9ep-54, -0x1.f0a7efb9230d8p-1},
+ {0x1.85f4e1b8298dp-54, -0x1.ed740e7684964p-1},
+ {0x1.a27d3874701f9p-53, -0x1.e9f4156c62ddcp-1},
+ {-0x1.ee94a90d7b88bp-53, -0x1.e6288ec48e11p-1},
+ {-0x1.7f59c49f6cd6dp-54, -0x1.e212104f686e4p-1},
+ {-0x1.83c37c6107db3p-55, -0x1.ddb13b6ccc23cp-1},
+ {0x1.d75033dfb9ca8p-53, -0x1.d906bcf328d48p-1},
+ {-0x1.ac42b5a8b6943p-53, -0x1.d4134d14dc938p-1},
+ {0x1.0f3db5dad5ac5p-53, -0x1.ced7af43cc774p-1},
+ {-0x1.ed0489e16b9ap-54, -0x1.c954b213411f4p-1},
+ {-0x1.a47d3a2a0dcbep-54, -0x1.c38b2f180bdbp-1},
+ {-0x1.9f6963354e3fep-53, -0x1.bd7c0ac6f9528p-1},
+ {-0x1.90e58336c64a8p-53, -0x1.b728345196e3cp-1},
+ {0x1.926da300ffccep-55, -0x1.b090a581502p-1},
+ {0x1.f98273c5d2495p-54, -0x1.a9b66290ea1a4p-1},
+ {-0x1.ddae89fd441d1p-53, -0x1.a29a7a046278p-1},
+ {-0x1.6788ebcc76dc6p-54, -0x1.9b3e047f3874p-1},
+ {0x1.615f32b6f907ap-54, -0x1.93a22499263fcp-1},
+ {-0x1.69d0f6897664ap-54, -0x1.8bc806b15174p-1},
+ {-0x1.d217be0e2b971p-53, -0x1.83b0e0bff976cp-1},
+ {0x1.21ea6f59be15bp-53, -0x1.7b5df226aafbp-1},
+ {0x1.bcac43c389ba9p-53, -0x1.72d0837efff98p-1},
+ {-0x1.21165f626cdd5p-54, -0x1.6a09e667f3bccp-1},
+ {0x1.24a366a5fe547p-53, -0x1.610b7551d2cep-1},
+ {0x1.75720992bfbb2p-55, -0x1.57d69348cecap-1},
+ {-0x1.13c293edceb32p-53, -0x1.4e6cabbe3e5e8p-1},
+ {0x1.e7f895d302395p-53, -0x1.44cf325091dd8p-1},
+ {-0x1.3c7c4bc72a92cp-53, -0x1.3affa292050b8p-1},
+ {-0x1.c20673b2116b2p-54, -0x1.30ff7fce17034p-1},
+ {0x1.5769d0fbcddc3p-53, -0x1.26d054cdd12ep-1},
+ {-0x1.b25dd267f66p-55, -0x1.1c73b39ae68c8p-1},
+ {0x1.7bc8eda6af93cp-53, -0x1.11eb3541b4b24p-1},
+ {-0x1.9697faf2e2fe5p-53, -0x1.073879922ffecp-1},
+ {-0x1.fb44f80f92225p-54, -0x1.f8ba4dbf89ab8p-2},
+ {0x1.c3e4edc5872f8p-55, -0x1.e2b5d3806f63cp-2},
+ {0x1.e97af1a63c807p-54, -0x1.cc66e9931c46p-2},
+ {-0x1.5b362cb974183p-57, -0x1.b5d1009e15ccp-2},
+ {0x1.d24afdade848bp-54, -0x1.9ef7943a8ed8cp-2},
+ {0x1.2e59dba7ab4c2p-54, -0x1.87de2a6aea964p-2},
+ {0x1.512c678219317p-54, -0x1.7088530fa45ap-2},
+ {-0x1.fc2047e54e614p-55, -0x1.58f9a75ab1fdcp-2},
+ {-0x1.4325f12be8946p-54, -0x1.4135c941766p-2},
+ {-0x1.a8b5c974ee7b5p-54, -0x1.294062ed59f04p-2},
+ {0x1.3ed9efaa42ab3p-55, -0x1.111d262b1f678p-2},
+ {0x1.50b7bbc4768b1p-55, -0x1.f19f97b215f1cp-3},
+ {0x1.035e2873ca432p-55, -0x1.c0b826a7e4f64p-3},
+ {0x1.49b466e7fe36p-55, -0x1.8f8b83c69a60cp-3},
+ {0x1.ab3802218894fp-55, -0x1.5e214448b3fc8p-3},
+ {0x1.dd9ffeaecbdc4p-55, -0x1.2c8106e8e613cp-3},
+ {0x1.cbb1f71aca352p-56, -0x1.f564e56a9731p-4},
+ {0x1.e2718d26ed688p-60, -0x1.917a6bc29b42cp-4},
+ {-0x1.ccbeeeae8129ap-56, -0x1.2d52092ce19f4p-4},
+ {0x1.912bd0d569a9p-61, -0x1.91f65f10dd814p-5},
+ {-0x1.f938a73db97fbp-58, -0x1.92155f7a3667cp-6},
+};
+
+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
+
+#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
new file mode 100644
index 0000000..5f2d8e7
--- /dev/null
+++ b/libc/src/math/generic/sin.cpp
@@ -0,0 +1,315 @@
+//===-- Double-precision sin function -------------------------------------===//
+//
+// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
+// See https://llvm.org/LICENSE.txt for license information.
+// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
+//
+//===----------------------------------------------------------------------===//
+
+#include "src/math/sin.h"
+#include "hdr/errno_macros.h"
+#include "src/__support/FPUtil/FEnvImpl.h"
+#include "src/__support/FPUtil/FPBits.h"
+#include "src/__support/FPUtil/PolyEval.h"
+#include "src/__support/FPUtil/double_double.h"
+#include "src/__support/FPUtil/dyadic_float.h"
+#include "src/__support/FPUtil/multiply_add.h"
+#include "src/__support/FPUtil/nearest_integer.h"
+#include "src/__support/FPUtil/rounding_mode.h"
+#include "src/__support/common.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/sincos_eval.h"
+
+#ifdef LIBC_TARGET_CPU_HAS_FMA
+#include "range_reduction_double_fma.h"
+
+using LIBC_NAMESPACE::fma::FAST_PASS_EXPONENT;
+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;
+
+LIBC_INLINE constexpr bool NO_FMA = false;
+#else
+#include "range_reduction_double_nofma.h"
+
+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;
+
+LIBC_INLINE constexpr bool NO_FMA = true;
+#endif // LIBC_TARGET_CPU_HAS_FMA
+
+// 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
+
+namespace LIBC_NAMESPACE {
+
+using DoubleDouble = fputil::DoubleDouble;
+using Float128 = typename fputil::DyadicFloat<128>;
+
+namespace {
+
+#ifndef LIBC_MATH_SIN_SKIP_ACCURATE_PASS
+LIBC_INLINE 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},
+ {Sign::POS, -131, 0x96a9'0496'70cf'ae65'f775'7409'4d3c'35c4_u128},
+ {Sign::POS, -131, 0xc8bd'35e1'4da1'5f0e'c739'6c89'4bbf'7389_u128},
+ {Sign::POS, -131, 0xfab2'72b5'4b98'71a2'7047'29ae'56d7'8a37_u128},
+ {Sign::POS, -130, 0x9640'8374'7309'd113'000a'89a1'1e07'c1fe_u128},
+ {Sign::POS, -130, 0xaf10'a224'59fe'32a6'3fee'f3bb'58b1'f10d_u128},
+ {Sign::POS, -130, 0xc7c5'c1e3'4d30'55b2'5cc8'c00e'4fcc'd850_u128},
+ {Sign::POS, -130, 0xe05c'1353'f27b'17e5'0ebc'61ad'e6ca'83cd_u128},
+ {Sign::POS, -130, 0xf8cf'cbd9'0af8'd57a'4221'dc4b'a772'598d_u128},
+ {Sign::POS, -129, 0x888e'9315'8fb3'bb04'9841'56f5'5334'4306_u128},
+ {Sign::POS, -129, 0x94a0'3176'acf8'2d45'ae4b'a773'da6b'f754_u128},
+ {Sign::POS, -129, 0xa09a'e4a0'bb30'0a19'2f89'5f44'a303'cc0b_u128},
+ {Sign::POS, -129, 0xac7c'd3ad'58fe'e7f0'811f'9539'84ef'f83e_u128},
+ {Sign::POS, -129, 0xb844'2987'd22c'f576'9cc3'ef36'746d'e3b8_u128},
+ {Sign::POS, -129, 0xc3ef'1535'754b'168d'3122'c2a5'9efd'dc37_u128},
+ {Sign::POS, -129, 0xcf7b'ca1d'476c'516d'a812'90bd'baad'62e4_u128},
+ {Sign::POS, -129, 0xdae8'804f'0ae6'015b'362c'b974'182e'3030_u128},
+ {Sign::POS, -129, 0xe633'74c9'8e22'f0b4'2872'ce1b'fc7a'd1cd_u128},
+ {Sign::POS, -129, 0xf15a'e9c0'37b1'd8f0'6c48'e9e3'420b'0f1e_u128},
+ {Sign::POS, -129, 0xfc5d'26df'c4d5'cfda'27c0'7c91'1290'b8d1_u128},
+ {Sign::POS, -128, 0x839c'3cc9'17ff'6cb4'bfd7'9717'f288'0abf_u128},
+ {Sign::POS, -128, 0x88f5'9aa0'da59'1421'b892'ca83'61d8'c84c_u128},
+ {Sign::POS, -128, 0x8e39'd9cd'7346'4364'bba4'cfec'bff5'4867_u128},
+ {Sign::POS, -128, 0x9368'2a66'e896'f544'b178'2191'1e71'c16e_u128},
+ {Sign::POS, -128, 0x987f'bfe7'0b81'a708'19ce'c845'ac87'a5c6_u128},
+ {Sign::POS, -128, 0x9d7f'd149'0285'c9e3'e25e'3954'9638'ae68_u128},
+ {Sign::POS, -128, 0xa267'9928'48ee'b0c0'3b51'67ee'359a'234e_u128},
+ {Sign::POS, -128, 0xa736'55df'1f2f'489e'149f'6e75'9934'68a3_u128},
+ {Sign::POS, -128, 0xabeb'49a4'6764'fd15'1bec'da80'89c1'a94c_u128},
+ {Sign::POS, -128, 0xb085'baa8'e966'f6da'e4ca'd00d'5c94'bcd2_u128},
+ {Sign::POS, -128, 0xb504'f333'f9de'6484'597d'89b3'754a'be9f_u128},
+ {Sign::POS, -128, 0xb968'41bf'7ffc'b21a'9de1'e3b2'2b8b'f4db_u128},
+ {Sign::POS, -128, 0xbdae'f913'557d'76f0'ac85'320f'528d'6d5d_u128},
+ {Sign::POS, -128, 0xc1d8'705f'fcbb'6e90'bdf0'715c'b8b2'0bd7_u128},
+ {Sign::POS, -128, 0xc5e4'0358'a8ba'05a7'43da'25d9'9267'326b_u128},
+ {Sign::POS, -128, 0xc9d1'124c'931f'da7a'8335'241b'e169'3225_u128},
+ {Sign::POS, -128, 0xcd9f'023f'9c3a'059e'23af'31db'7179'a4aa_u128},
+ {Sign::POS, -128, 0xd14d'3d02'313c'0eed'744f'ea20'e8ab'ef92_u128},
+ {Sign::POS, -128, 0xd4db'3148'750d'1819'f630'e8b6'dac8'3e69_u128},
+ {Sign::POS, -128, 0xd848'52c0'a80f'fcdb'24b9'fe00'6635'74a4_u128},
+ {Sign::POS, -128, 0xdb94'1a28'cb71'ec87'2c19'b632'53da'43fc_u128},
+ {Sign::POS, -128, 0xdebe'0563'7ca9'4cfb'4b19'aa71'fec3'ae6d_u128},
+ {Sign::POS, -128, 0xe1c5'978c'05ed'8691'f4e8'a837'2f8c'5810_u128},
+ {Sign::POS, -128, 0xe4aa'5909'a08f'a7b4'1227'85ae'67f5'515d_u128},
+ {Sign::POS, -128, 0xe76b'd7a1'e63b'9786'1251'2952'9d48'a92f_u128},
+ {Sign::POS, -128, 0xea09'a68a'6e49'cd62'15ad'45b4'a1b5'e823_u128},
+ {Sign::POS, -128, 0xec83'5e79'946a'3145'7e61'0231'ac1d'6181_u128},
+ {Sign::POS, -128, 0xeed8'9db6'6611'e307'86f8'c20f'b664'b01b_u128},
+ {Sign::POS, -128, 0xf109'0827'b437'25fd'6712'7db3'5b28'7316_u128},
+ {Sign::POS, -128, 0xf314'4762'4708'8f74'a548'6bdc'455d'56a2_u128},
+ {Sign::POS, -128, 0xf4fa'0ab6'316e'd2ec'163c'5c7f'03b7'18c5_u128},
+ {Sign::POS, -128, 0xf6ba'073b'424b'19e8'2c79'1f59'cc1f'fc23_u128},
+ {Sign::POS, -128, 0xf853'f7dc'9186'b952'c7ad'c6b4'9888'91bb_u128},
+ {Sign::POS, -128, 0xf9c7'9d63'272c'4628'4504'ae08'd19b'2980_u128},
+ {Sign::POS, -128, 0xfb14'be7f'bae5'8156'2172'a361'fd2a'722f_u128},
+ {Sign::POS, -128, 0xfc3b'27d3'8a5d'49ab'2567'78ff'cb5c'1769_u128},
+ {Sign::POS, -128, 0xfd3a'abf8'4528'b50b'eae6'bd95'1c1d'abbe_u128},
+ {Sign::POS, -128, 0xfe13'2387'0cfe'9a3d'90cd'1d95'9db6'74ef_u128},
+ {Sign::POS, -128, 0xfec4'6d1e'8929'2cf0'4139'0efd'c726'e9ef_u128},
+ {Sign::POS, -128, 0xff4e'6d68'0c41'd0a9'0f66'8633'f1ab'858a_u128},
+ {Sign::POS, -128, 0xffb1'0f1b'cb6b'ef1d'421e'8eda'af59'453e_u128},
+ {Sign::POS, -128, 0xffec'4304'2668'65d9'5657'5523'6696'1732_u128},
+ {Sign::POS, 0, 1},
+};
+
+#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
+
+#endif // !LIBC_MATH_SIN_SKIP_ACCURATE_PASS
+
+} // anonymous namespace
+
+LLVM_LIBC_FUNCTION(double, sin, (double x)) {
+ using FPBits = typename fputil::FPBits<double>;
+ FPBits xbits(x);
+
+ uint16_t x_e = xbits.get_biased_exponent();
+
+ DoubleDouble y;
+ unsigned k;
+ generic::LargeRangeReduction<NO_FMA> range_reduction_large;
+
+ // |x| < 2^32 (with FMA) or |x| < 2^23 (w/o FMA)
+ 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;
+
+ // 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
+ 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
+ }
+
+ // // Small range reduction.
+ k = range_reduction_small(x, y);
+ } else {
+ // Inf or NaN
+ if (LIBC_UNLIKELY(x_e > 2 * FPBits::EXP_BIAS)) {
+ // sin(+-Inf) = NaN
+ if (xbits.get_mantissa() == 0) {
+ fputil::set_errno_if_required(EDOM);
+ fputil::raise_except_if_required(FE_INVALID);
+ }
+ return x + FPBits::quiet_nan().get_val();
+ }
+
+ // Large range reduction.
+ k = range_reduction_large.compute_high_part(x);
+ y = range_reduction_large.fast();
+ }
+
+ DoubleDouble sin_y, cos_y;
+
+ 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);
+
+ // 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];
+
+ // 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);
+
+ FPBits sk_cy(sin_k_cos_y.hi);
+ FPBits ck_sy(cos_k_sin_y.hi);
+ 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
+ return rr.hi + rr.lo;
+#else
+ // Accurate test and pass for correctly rounded implementation.
+ 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);
+
+ // Ziv's rounding test.
+ if (LIBC_LIKELY(r_upper == r_lower))
+ return r_upper;
+
+ Float128 u_f128;
+ if (LIBC_LIKELY(x_e < FPBits::EXP_BIAS + FAST_PASS_EXPONENT))
+ u_f128 = generic::range_reduction_small_f128(x);
+ else
+ u_f128 = range_reduction_large.accurate();
+
+ Float128 u_sq = fputil::quick_mul(u_f128, u_f128);
+
+ // sin(u) ~ x - x^3/3! + x^5/5! - x^7/7! + x^9/9! - x^11/11! + x^13/13!
+ constexpr Float128 SIN_COEFFS[] = {
+ {Sign::POS, -127, 0x80000000'00000000'00000000'00000000_u128}, // 1
+ {Sign::NEG, -130, 0xaaaaaaaa'aaaaaaaa'aaaaaaaa'aaaaaaab_u128}, // -1/3!
+ {Sign::POS, -134, 0x88888888'88888888'88888888'88888889_u128}, // 1/5!
+ {Sign::NEG, -140, 0xd00d00d0'0d00d00d'00d00d00'd00d00d0_u128}, // -1/7!
+ {Sign::POS, -146, 0xb8ef1d2a'b6399c7d'560e4472'800b8ef2_u128}, // 1/9!
+ {Sign::NEG, -153, 0xd7322b3f'aa271c7f'3a3f25c1'bee38f10_u128}, // -1/11!
+ {Sign::POS, -160, 0xb092309d'43684be5'1c198e91'd7b4269e_u128}, // 1/13!
+ };
+
+ // cos(u) ~ 1 - x^2/2 + x^4/4! - x^6/6! + x^8/8! - x^10/10! + x^12/12!
+ constexpr Float128 COS_COEFFS[] = {
+ {Sign::POS, -127, 0x80000000'00000000'00000000'00000000_u128}, // 1.0
+ {Sign::NEG, -128, 0x80000000'00000000'00000000'00000000_u128}, // 1/2
+ {Sign::POS, -132, 0xaaaaaaaa'aaaaaaaa'aaaaaaaa'aaaaaaab_u128}, // 1/4!
+ {Sign::NEG, -137, 0xb60b60b6'0b60b60b'60b60b60'b60b60b6_u128}, // 1/6!
+ {Sign::POS, -143, 0xd00d00d0'0d00d00d'00d00d00'd00d00d0_u128}, // 1/8!
+ {Sign::NEG, -149, 0x93f27dbb'c4fae397'780b69f5'333c725b_u128}, // 1/10!
+ {Sign::POS, -156, 0x8f76c77f'c6c4bdaa'26d4c3d6'7f425f60_u128}, // 1/12!
+ };
+
+ Float128 sin_u = fputil::quick_mul(
+ u_f128, fputil::polyeval(u_sq, SIN_COEFFS[0], SIN_COEFFS[1],
+ SIN_COEFFS[2], SIN_COEFFS[3], SIN_COEFFS[4],
+ SIN_COEFFS[5], SIN_COEFFS[6]));
+ Float128 cos_u = fputil::polyeval(u_sq, COS_COEFFS[0], COS_COEFFS[1],
+ COS_COEFFS[2], COS_COEFFS[3], COS_COEFFS[4],
+ COS_COEFFS[5], COS_COEFFS[6]);
+
+ auto get_sin_k = [](unsigned kk) -> Float128 {
+ unsigned idx = (kk & 64) ? 64 - (kk & 63) : (kk & 63);
+ Float128 ans = SIN_K_PI_OVER_128_F128[idx];
+ if (kk & 128)
+ ans.sign = Sign::NEG;
+ return ans;
+ };
+
+ // cos(k * pi/128) = sin(k * pi/128 + pi/2) = sin((k + 64) * pi/128).
+ Float128 sin_k_f128 = get_sin_k(k);
+ Float128 cos_k_f128 = get_sin_k(k + 64);
+
+ // sin(x) = sin((k * pi/128 + u)
+ // = sin(u) * cos(k*pi/128) + cos(u) * sin(k*pi/128)
+ Float128 r = fputil::quick_add(fputil::quick_mul(sin_k_f128, cos_u),
+ fputil::quick_mul(cos_k_f128, sin_u));
+
+ // TODO: Add assertion if Ziv's accuracy tests fail in debug mode.
+ // https://github.com/llvm/llvm-project/issues/96452.
+
+ return static_cast<double>(r);
+#endif // !LIBC_MATH_SIN_SKIP_ACCURATE_PASS
+}
+
+} // namespace LIBC_NAMESPACE
diff --git a/libc/src/math/generic/sincos_eval.h b/libc/src/math/generic/sincos_eval.h
new file mode 100644
index 0000000..d5db18f
--- /dev/null
+++ b/libc/src/math/generic/sincos_eval.h
@@ -0,0 +1,81 @@
+//===-- Compute sin + cos for small angles ----------------------*- C++ -*-===//
+//
+// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
+// See https://llvm.org/LICENSE.txt for license information.
+// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
+//
+//===----------------------------------------------------------------------===//
+
+#ifndef LLVM_LIBC_SRC_MATH_GENERIC_SINCOS_EVAL_H
+#define LLVM_LIBC_SRC_MATH_GENERIC_SINCOS_EVAL_H
+
+#include "src/__support/FPUtil/double_double.h"
+#include "src/__support/FPUtil/multiply_add.h"
+
+namespace LIBC_NAMESPACE {
+
+using fputil::DoubleDouble;
+
+LIBC_INLINE void 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!
+ // Then the error is bounded by:
+ // |sin(y) - (y - y^3/3! + y^5/5! - y^7/7!)| < |y|^9/9! < 2^-54/9! < 2^-72.
+ // For y ~ u_hi + u_lo, fully expanding the polynomial and drop any terms
+ // < ulp(u_hi^3) gives us:
+ // y - y^3/3! + y^5/5! - y^7/7! = ...
+ // ~ u_hi + u_hi^3 * (-1/6 + u_hi^2 * (1/120 - u_hi^2 * 1/5040)) +
+ // + u_lo (1 + u_hi^2 * (-1/2 + u_hi^2 / 24))
+ double u_hi_sq = u.hi * u.hi; // Error < ulp(u_hi^2) < 2^(-6 - 52) = 2^-58.
+ // p1 ~ 1/120 + u_hi^2 / 5040.
+ double p1 = fputil::multiply_add(u_hi_sq, -0x1.a01a01a01a01ap-13,
+ 0x1.1111111111111p-7);
+ // q1 ~ -1/2 + u_hi^2 / 24.
+ double q1 = fputil::multiply_add(u_hi_sq, 0x1.5555555555555p-5, -0x1.0p-1);
+ double u_hi_3 = u_hi_sq * u.hi;
+ // p2 ~ -1/6 + u_hi^2 (1/120 - u_hi^2 * 1/5040)
+ double p2 = fputil::multiply_add(u_hi_sq, p1, -0x1.5555555555555p-3);
+ // q2 ~ 1 + u_hi^2 (-1/2 + u_hi^2 / 24)
+ double q2 = fputil::multiply_add(u_hi_sq, q1, 1.0);
+ double sin_lo = fputil::multiply_add(u_hi_3, p2, u.lo * q2);
+ // Overall, |sin(y) - (u_hi + sin_lo)| < 2*ulp(u_hi^3) < 2^-69.
+
+ // Evaluate cos(y) = cos(x - k * (pi/128))
+ // We use the degree-8 Taylor approximation:
+ // cos(y) ~ 1 - y^2/2 + y^4/4! - y^6/6! + y^8/8!
+ // Then the error is bounded by:
+ // |cos(y) - (...)| < |y|^10/10! < 2^-81
+ // For y ~ u_hi + u_lo, fully expanding the polynomial and drop any terms
+ // < ulp(u_hi^3) gives us:
+ // 1 - y^2/2 + y^4/4! - y^6/6! + y^8/8! = ...
+ // ~ 1 - u_hi^2/2 + u_hi^4(1/24 + u_hi^2 (-1/720 + u_hi^2/40320)) +
+ // + 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);
+
+ // r1 ~ -1/720 + u_hi^2 / 40320
+ double r1 = fputil::multiply_add(u_hi_sq, 0x1.a01a01a01a01ap-16,
+ -0x1.6c16c16c16c17p-10);
+ // s1 ~ -1 + u_hi^2 / 6
+ double s1 = fputil::multiply_add(u_hi_sq, 0x1.5555555555555p-3, -1.0);
+ double u_hi_4 = u_hi_sq * u_hi_sq;
+ double u_hi_u_lo = u.hi * u.lo;
+ // 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 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);
+}
+
+} // namespace LIBC_NAMESPACE
+
+#endif // LLVM_LIBC_SRC_MATH_GENERIC_SINCOSF_EVAL_H
diff --git a/libc/src/math/x86_64/CMakeLists.txt b/libc/src/math/x86_64/CMakeLists.txt
index cd129e3..882181b 100644
--- a/libc/src/math/x86_64/CMakeLists.txt
+++ b/libc/src/math/x86_64/CMakeLists.txt
@@ -9,16 +9,6 @@ add_entrypoint_object(
)
add_entrypoint_object(
- sin
- SRCS
- sin.cpp
- HDRS
- ../sin.h
- COMPILE_OPTIONS
- -O2
-)
-
-add_entrypoint_object(
tan
SRCS
tan.cpp
diff --git a/libc/src/math/x86_64/sin.cpp b/libc/src/math/x86_64/sin.cpp
deleted file mode 100644
index 2c7b8aa..0000000
--- a/libc/src/math/x86_64/sin.cpp
+++ /dev/null
@@ -1,19 +0,0 @@
-//===-- Implementation of the sin function for x86_64 ---------------------===//
-//
-// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
-// See https://llvm.org/LICENSE.txt for license information.
-// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
-//
-//===----------------------------------------------------------------------===//
-
-#include "src/math/sin.h"
-#include "src/__support/common.h"
-
-namespace LIBC_NAMESPACE {
-
-LLVM_LIBC_FUNCTION(double, sin, (double x)) {
- __asm__ __volatile__("fsin" : "+t"(x));
- return x;
-}
-
-} // namespace LIBC_NAMESPACE
diff --git a/libc/test/src/math/sin_test.cpp b/libc/test/src/math/sin_test.cpp
index 0171b79..89534ae 100644
--- a/libc/test/src/math/sin_test.cpp
+++ b/libc/test/src/math/sin_test.cpp
@@ -12,22 +12,103 @@
#include "test/UnitTest/Test.h"
#include "utils/MPFRWrapper/MPFRUtils.h"
-#include "hdr/math_macros.h"
-
using LlvmLibcSinTest = LIBC_NAMESPACE::testing::FPTest<double>;
namespace mpfr = LIBC_NAMESPACE::testing::mpfr;
-TEST_F(LlvmLibcSinTest, Range) {
- static constexpr double _2pi = 6.283185307179586;
- constexpr StorageType COUNT = 100'000;
- constexpr StorageType STEP = STORAGE_MAX / COUNT;
- for (StorageType i = 0, v = 0; i <= COUNT; ++i, v += STEP) {
- double x = FPBits(v).get_val();
- // TODO: Expand the range of testing after range reduction is implemented.
- if (isnan(x) || isinf(x) || x > _2pi || x < -_2pi)
- continue;
-
- ASSERT_MPFR_MATCH(mpfr::Operation::Sin, x, LIBC_NAMESPACE::sin(x), 1.0);
+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.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,
+ 0x1.2e566149bf5fdp+9, 0x1.be886d9c2324dp+6, -0x1.119471e9216cdp+10,
+ -0x1.aaf85537ea4c7p+3, 0x1.cb996c60f437ep+9, 0x1.c96e28eb679f8p+5,
+ -0x1.a5eece87e8606p+4, 0x1.e31b55306f22cp+2, 0x1.ae78d360afa15p+0,
+ 0x1.1685973506319p+3, 0x1.4f2b874135d27p+4, 0x1.ae945054939c2p+10,
+ 0x1.3eec5912ea7cdp+331, 0x1.dcbfda0c7559ep+8, 0x1.a65d441ea6dcep+4,
+ 0x1.e639103a05997p+2, 0x1.13114266f9764p+4, -0x1.3eec5912ea7cdp+331,
+ 0x1.08087e9aad90bp+887, 0x1.2b5fe88a9d8d5p+903, -0x1.a880417b7b119p+1023,
+ -0x1.6deb37da81129p+205, 0x1.08087e9aad90bp+887, 0x1.f6d7518808571p+1023,
+ -0x1.8bb5847d49973p+845, 0x1.f08b14e1c4d0fp+890, 0x1.6ac5b262ca1ffp+849,
+ 0x1.e0000000001c2p-20,
+ };
+ constexpr int N = sizeof(INPUTS) / sizeof(INPUTS[0]);
+
+ for (int i = 0; i < N; ++i) {
+ double x = INPUTS[i];
+ ASSERT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Sin, x,
+ LIBC_NAMESPACE::sin(x), 0.5);
}
}
+
+TEST_F(LlvmLibcSinTest, InDoubleRange) {
+ constexpr uint64_t COUNT = 1'234'51;
+ uint64_t START = LIBC_NAMESPACE::fputil::FPBits<double>(0x1.0p-50).uintval();
+ uint64_t STOP = LIBC_NAMESPACE::fputil::FPBits<double>(0x1.0p200).uintval();
+ uint64_t STEP = (STOP - START) / COUNT;
+
+ auto test = [&](mpfr::RoundingMode rounding_mode) {
+ mpfr::ForceRoundingMode __r(rounding_mode);
+ if (!__r.success)
+ return;
+
+ uint64_t fails = 0;
+ uint64_t count = 0;
+ uint64_t cc = 0;
+ double mx, mr = 0.0;
+ double tol = 0.5;
+
+ for (uint64_t i = 0, v = START; i <= COUNT; ++i, v += STEP) {
+ double x = FPBits(v).get_val();
+ if (isnan(x) || isinf(x))
+ continue;
+ LIBC_NAMESPACE::libc_errno = 0;
+ double result = LIBC_NAMESPACE::sin(x);
+ ++cc;
+ if (isnan(result) || isinf(result))
+ continue;
+
+ ++count;
+
+ if (!TEST_MPFR_MATCH_ROUNDING_SILENTLY(mpfr::Operation::Sin, x, result,
+ 0.5, rounding_mode)) {
+ ++fails;
+ while (!TEST_MPFR_MATCH_ROUNDING_SILENTLY(mpfr::Operation::Sin, x,
+ result, tol, rounding_mode)) {
+ mx = x;
+ mr = result;
+
+ if (tol > 1000.0)
+ break;
+
+ tol *= 2.0;
+ }
+ }
+ }
+ if (fails) {
+ tlog << " Sin failed: " << fails << "/" << count << "/" << cc
+ << " tests.\n";
+ tlog << " Max ULPs is at most: " << static_cast<uint64_t>(tol) << ".\n";
+ EXPECT_MPFR_MATCH(mpfr::Operation::Sin, mx, mr, 0.5, rounding_mode);
+ }
+ };
+
+ tlog << " Test Rounding To Nearest...\n";
+ test(mpfr::RoundingMode::Nearest);
+
+ tlog << " Test Rounding Downward...\n";
+ test(mpfr::RoundingMode::Downward);
+
+ tlog << " Test Rounding Upward...\n";
+ test(mpfr::RoundingMode::Upward);
+
+ tlog << " Test Rounding Toward Zero...\n";
+ test(mpfr::RoundingMode::TowardZero);
+}
diff --git a/libc/test/src/math/smoke/CMakeLists.txt b/libc/test/src/math/smoke/CMakeLists.txt
index 8cae52f..ffa8ba1 100644
--- a/libc/test/src/math/smoke/CMakeLists.txt
+++ b/libc/test/src/math/smoke/CMakeLists.txt
@@ -3653,3 +3653,13 @@ add_fp_unittest(
DEPENDS
libc.src.math.f16sqrtf
)
+
+add_fp_unittest(
+ sin_test
+ SUITE
+ libc-math-smoke-tests
+ SRCS
+ sin_test.cpp
+ DEPENDS
+ libc.src.math.sin
+)
diff --git a/libc/test/src/math/smoke/sin_test.cpp b/libc/test/src/math/smoke/sin_test.cpp
new file mode 100644
index 0000000..16ced68
--- /dev/null
+++ b/libc/test/src/math/smoke/sin_test.cpp
@@ -0,0 +1,26 @@
+//===-- Unittests for sin -------------------------------------------------===//
+//
+// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
+// See https://llvm.org/LICENSE.txt for license information.
+// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
+//
+//===----------------------------------------------------------------------===//
+
+#include "src/math/sin.h"
+#include "test/UnitTest/FPMatcher.h"
+#include "test/UnitTest/Test.h"
+
+using LlvmLibcSinTest = LIBC_NAMESPACE::testing::FPTest<double>;
+
+using LIBC_NAMESPACE::testing::tlog;
+
+TEST_F(LlvmLibcSinTest, SpecialNumbers) {
+ EXPECT_FP_EQ_ALL_ROUNDING(aNaN, LIBC_NAMESPACE::sin(aNaN));
+ EXPECT_FP_EQ_ALL_ROUNDING(aNaN, LIBC_NAMESPACE::sin(inf));
+ EXPECT_FP_EQ_ALL_ROUNDING(aNaN, LIBC_NAMESPACE::sin(neg_inf));
+ EXPECT_FP_EQ_ALL_ROUNDING(zero, LIBC_NAMESPACE::sin(zero));
+ EXPECT_FP_EQ_ALL_ROUNDING(neg_zero, LIBC_NAMESPACE::sin(neg_zero));
+ EXPECT_FP_EQ(0x1.0p-50, LIBC_NAMESPACE::sin(0x1.0p-50));
+ EXPECT_FP_EQ(min_normal, LIBC_NAMESPACE::sin(min_normal));
+ EXPECT_FP_EQ(min_denormal, LIBC_NAMESPACE::sin(min_denormal));
+}
generated by cgit v1.1 at 2024-07-08 10:16:07 +0000