From 16903ace180755b7558234ff2b2e8d89b00dcb88 Mon Sep 17 00:00:00 2001 From: lntue <35648136+lntue@users.noreply.github.com> Date: Mon, 24 Jun 2024 17:57:08 -0400 Subject: [PATCH] [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'. --- libc/config/darwin/arm/entrypoints.txt | 1 + libc/config/linux/aarch64/entrypoints.txt | 1 + libc/config/linux/arm/entrypoints.txt | 1 + libc/config/linux/riscv/entrypoints.txt | 1 + libc/docs/math/index.rst | 2 +- libc/src/__support/FPUtil/double_double.h | 52 +- libc/src/__support/FPUtil/dyadic_float.h | 10 +- libc/src/__support/macros/optimization.h | 14 + libc/src/math/generic/CMakeLists.txt | 49 ++ .../generic/range_reduction_double_common.h | 162 ++++++ .../math/generic/range_reduction_double_fma.h | 495 ++++++++++++++++++ .../generic/range_reduction_double_nofma.h | 493 +++++++++++++++++ libc/src/math/generic/sin.cpp | 315 +++++++++++ libc/src/math/generic/sincos_eval.h | 81 +++ libc/src/math/x86_64/CMakeLists.txt | 10 - libc/src/math/x86_64/sin.cpp | 19 - libc/test/src/math/sin_test.cpp | 107 +++- libc/test/src/math/smoke/CMakeLists.txt | 10 + libc/test/src/math/smoke/sin_test.cpp | 26 + 19 files changed, 1790 insertions(+), 59 deletions(-) create mode 100644 libc/src/math/generic/range_reduction_double_common.h create mode 100644 libc/src/math/generic/range_reduction_double_fma.h create mode 100644 libc/src/math/generic/range_reduction_double_nofma.h create mode 100644 libc/src/math/generic/sin.cpp create mode 100644 libc/src/math/generic/sincos_eval.h delete mode 100644 libc/src/math/x86_64/sin.cpp create mode 100644 libc/test/src/math/smoke/sin_test.cpp diff --git a/libc/config/darwin/arm/entrypoints.txt b/libc/config/darwin/arm/entrypoints.txt index e1303265b9ac41..d486916a542d6d 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 63238ec0094fbf..2d2e1b58c0b271 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 00df44511a4c1c..53cdcc3caab1c9 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 e12d6b3957e518..18968f5b07b596 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 d6e642c4907f04..29723be7ce7526 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 b9490b52f6b41e..3d16a3cce3a99c 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; // 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 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(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 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(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 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 LIBC_INLINE DoubleDouble quick_mult(const DoubleDouble &a, const DoubleDouble &b) { - DoubleDouble r = exact_mult(a.hi, b.hi); + DoubleDouble r = exact_mult(a.hi, b.hi); double t1 = multiply_add(a.hi, b.lo, r.lo); double t2 = multiply_add(a.lo, b.hi, t1); r.lo = t2; diff --git a/libc/src/__support/FPUtil/dyadic_float.h b/libc/src/__support/FPUtil/dyadic_float.h index 63cb98351201f6..76786a24f3e922 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 quick_add(DyadicFloat 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 -LIBC_INLINE constexpr DyadicFloat quick_mul(DyadicFloat a, - DyadicFloat b) { +LIBC_INLINE constexpr DyadicFloat quick_mul(const DyadicFloat &a, + const DyadicFloat &b) { DyadicFloat 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(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 &a, const DyadicFloat &b, // Simple exponentiation implementation for printf. Only handles positive // exponents, since division isn't implemented. template -LIBC_INLINE constexpr DyadicFloat pow_n(DyadicFloat a, +LIBC_INLINE constexpr DyadicFloat pow_n(const DyadicFloat &a, uint32_t power) { DyadicFloat result = 1.0; DyadicFloat cur_power = a; @@ -325,7 +325,7 @@ LIBC_INLINE constexpr DyadicFloat pow_n(DyadicFloat a, } template -LIBC_INLINE constexpr DyadicFloat mul_pow_2(DyadicFloat a, +LIBC_INLINE constexpr DyadicFloat mul_pow_2(const DyadicFloat &a, int32_t pow_2) { DyadicFloat result = a; result.exponent += pow_2; diff --git a/libc/src/__support/macros/optimization.h b/libc/src/__support/macros/optimization.h index 59886ca44be12a..05a47791deed8a 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 a0114aafcad39e..54a5b6ac0a70c4 100644 --- a/libc/src/math/generic/CMakeLists.txt +++ b/libc/src/math/generic/CMakeLists.txt @@ -135,6 +135,23 @@ add_header_library( libc.src.__support.common ) +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 @@ -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 @@ -167,6 +193,29 @@ add_entrypoint_object( -O3 ) +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 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 00000000000000..0e9edf811891d6 --- /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 struct LargeRangeReduction { + // Calculate the high part of the range reduction exactly. + LIBC_INLINE unsigned compute_high_part(double x) { + using FPBits = typename fputil::FPBits; + FPBits xbits(x); + + // TODO: The extra exponent gap of 62 below can be reduced a bit for non-FMA + // with a more careful analysis, which in turn will reduce the error bound + // for non-FMA + int x_e_m62 = xbits.get_biased_exponent() - (FPBits::EXP_BIAS + 62); + idx = static_cast((x_e_m62 >> 4) + 3); + // Scale x down by 2^(-(16 * (idx - 3)) + xbits.set_biased_exponent((x_e_m62 & 15) + FPBits::EXP_BIAS + 62); + // 2^62 <= |x_reduced| < 2^(62 + 16) = 2^78 + x_reduced = xbits.get_val(); + // x * c_hi = ph.hi + ph.lo exactly. + DoubleDouble ph = + fputil::exact_mult(x_reduced, ONE_TWENTY_EIGHT_OVER_PI[idx][0]); + // x * c_mid = pm.hi + pm.lo exactly. + DoubleDouble pm = + fputil::exact_mult(x_reduced, ONE_TWENTY_EIGHT_OVER_PI[idx][1]); + // Extract integral parts and fractional parts of (ph.lo + pm.hi). + double kh = fputil::nearest_integer(ph.lo); + double ph_lo_frac = ph.lo - kh; // Exact + double km = fputil::nearest_integer(pm.hi + ph_lo_frac); + double pm_hi_frac = pm.hi - km; // Exact + // x * 128/pi mod 1 ~ y_hi + y_lo + y_hi = ph_lo_frac + pm_hi_frac; // Exact + pm_lo = pm.lo; + return static_cast(static_cast(kh) + + static_cast(km)); + } + + LIBC_INLINE DoubleDouble fast() const { + // y_lo = x * c_lo + pm.lo + double y_lo = fputil::multiply_add(x_reduced, + ONE_TWENTY_EIGHT_OVER_PI[idx][2], pm_lo); + DoubleDouble y = fputil::exact_add(y_hi, y_lo); + + // Digits of pi/128, generated by Sollya with: + // > a = round(pi/128, D, RN); + // > b = round(pi/128 - a, D, RN); + constexpr DoubleDouble PI_OVER_128_DD = {0x1.1a62633145c07p-60, + 0x1.921fb54442d18p-6}; + + // Error bound: with {a} denote the fractional part of a, i.e.: + // {a} = a - round(a) + // Then, + // | {x * 128/pi} - (y_hi + y_lo) | < 2 * ulp(x_reduced * + // * ONE_TWENTY_EIGHT_OVER_PI[idx][2]) + // For FMA: + // | {x * 128/pi} - (y_hi + y_lo) | <= 2 * 2^77 * 2^-103 * 2^-52 + // = 2^-77. + // | {x mod pi/128} - (u.hi + u.lo) | < 2 * 2^-6 * 2^-77. + // = 2^-82. + // For non-FMA: + // | {x * 128/pi} - (y_hi + y_lo) | <= 2 * 2^77 * 2^-99 * 2^-52 + // = 2^-73. + // | {x mod pi/128} - (u.hi + u.lo) | < 2 * 2^-6 * 2^-73. + // = 2^-78. + return fputil::quick_mult(y, PI_OVER_128_DD); + } + + 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 00000000000000..c136de957d2f05 --- /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(static_cast(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 00000000000000..b9d34d6d9f1a8f --- /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(static_cast(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 00000000000000..5f2d8e7bd2d648 --- /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; + FPBits xbits(x); + + uint16_t x_e = xbits.get_biased_exponent(); + + DoubleDouble y; + unsigned k; + generic::LargeRangeReduction 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(k & 128) << (63 - 7); + // sin_k.hi = FPBits(FPBits(sin_k.hi).uintval() ^ sin_s).get_val(); + // sin_k.lo = FPBits(FPBits(sin_k.hi).uintval() ^ sin_s).get_val(); + // DoubleDouble cos_k = SIN_K_PI_OVER_128[(k + 64) & 127]; + // uint64_t cos_s = static_cast((k + 64) & 128) << (63 - 7); + // cos_k.hi = FPBits(FPBits(cos_k.hi).uintval() ^ cos_s).get_val(); + // cos_k.lo = FPBits(FPBits(cos_k.hi).uintval() ^ cos_s).get_val(); + + // Use 64-entry table instead: + // auto get_idx_dd = [](unsigned kk) -> DoubleDouble { + // unsigned idx = (kk & 64) ? 64 - (kk & 63) : (kk & 63); + // DoubleDouble ans = SIN_K_PI_OVER_128[idx]; + // if (kk & 128) { + // ans.hi = -ans.hi; + // ans.lo = -ans.lo; + // } + // return ans; + // }; + // DoubleDouble sin_k = get_idx_dd(k); + // DoubleDouble cos_k = get_idx_dd(k + 64); + + // 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(cos_y, sin_k); + DoubleDouble cos_k_sin_y = fputil::quick_mult(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(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(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 00000000000000..d5db18f04a8f17 --- /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 cd129e3eefb75c..882181b33b9f88 100644 --- a/libc/src/math/x86_64/CMakeLists.txt +++ b/libc/src/math/x86_64/CMakeLists.txt @@ -8,16 +8,6 @@ add_entrypoint_object( -O2 ) -add_entrypoint_object( - sin - SRCS - sin.cpp - HDRS - ../sin.h - COMPILE_OPTIONS - -O2 -) - add_entrypoint_object( tan SRCS diff --git a/libc/src/math/x86_64/sin.cpp b/libc/src/math/x86_64/sin.cpp deleted file mode 100644 index 2c7b8aa6e8c836..00000000000000 --- 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 0171b79810d4e9..89534ae82b56e3 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; 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(0x1.0p-50).uintval(); + uint64_t STOP = LIBC_NAMESPACE::fputil::FPBits(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(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 8cae52fb0cebe5..ffa8ba138e309c 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 00000000000000..16ced68709ca75 --- /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; + +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)); +}