[libc-commits] [libc] [libc][math] Implement double precision sin correctly rounded to all rounding modes. (PR #95736)
via libc-commits
libc-commits at lists.llvm.org
Wed Jun 19 06:45:30 PDT 2024
https://github.com/lntue updated https://github.com/llvm/llvm-project/pull/95736
>From ffaf0ba65fdc4503de1c305433cf68091de2da13 Mon Sep 17 00:00:00 2001
From: Tue Ly <lntue.h at gmail.com>
Date: Mon, 17 Jun 2024 04:51:20 +0000
Subject: [PATCH 1/3] [libc][math] Implement double precision sin correctly
rounded to all rounding modes.
---
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 | 10 +-
libc/src/__support/FPUtil/dyadic_float.h | 10 +-
libc/src/__support/macros/optimization.h | 14 +
libc/src/math/generic/CMakeLists.txt | 48 ++
.../src/math/generic/range_reduction_double.h | 67 +++
.../math/generic/range_reduction_double_fma.h | 334 +++++++++++
libc/src/math/generic/sin.cpp | 567 ++++++++++++++++++
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 | 106 +++-
libc/test/src/math/smoke/CMakeLists.txt | 10 +
libc/test/src/math/smoke/sin_test.cpp | 26 +
18 files changed, 1259 insertions(+), 49 deletions(-)
create mode 100644 libc/src/math/generic/range_reduction_double.h
create mode 100644 libc/src/math/generic/range_reduction_double_fma.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 e1303265b9ac4..d486916a542d6 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 dfed6acbdf257..07d6583c0a816 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 d4f932416bd9f..5716b777c8c96 100644
--- a/libc/config/linux/arm/entrypoints.txt
+++ b/libc/config/linux/arm/entrypoints.txt
@@ -345,6 +345,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 e12d6b3957e51..18968f5b07b59 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 293edd1c15100..48867d7012b51 100644
--- a/libc/docs/math/index.rst
+++ b/libc/docs/math/index.rst
@@ -314,7 +314,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 b9490b52f6b41..a94902b7ec7ae 100644
--- a/libc/src/__support/FPUtil/double_double.h
+++ b/libc/src/__support/FPUtil/double_double.h
@@ -44,7 +44,12 @@ LIBC_INLINE constexpr DoubleDouble add(const DoubleDouble &a, double b) {
return exact_add(r.hi, r.lo + a.lo);
}
-// Velkamp's Splitting for double precision.
+// Veltkamp's Splitting for double precision.
+// Note: This is the original version of Veltkamp's Splitting, which is only
+// correct for round-to-nearest mode. See:
+// Graillat, S., Lafevre, V., and Muller, J.-M., "Alternative Split Functions
+// and Dekker's Product," ARITH'2020.
+// http://arith2020.arithsymposium.org/resources/paper_31.pdf
LIBC_INLINE constexpr DoubleDouble split(double a) {
DoubleDouble r{0.0, 0.0};
// Splitting constant = 2^ceil(prec(double)/2) + 1 = 2^27 + 1.
@@ -56,6 +61,9 @@ LIBC_INLINE constexpr DoubleDouble split(double a) {
return r;
}
+// Note: When FMA instruction is not available, the `exact_mult` function relies
+// on Veltkamp's Splitting algorithm, and is only correct for round-to-nearest
+// mode.
LIBC_INLINE DoubleDouble exact_mult(double a, double b) {
DoubleDouble r{0.0, 0.0};
diff --git a/libc/src/__support/FPUtil/dyadic_float.h b/libc/src/__support/FPUtil/dyadic_float.h
index 12a69228d36c7..9e25cc6486011 100644
--- a/libc/src/__support/FPUtil/dyadic_float.h
+++ b/libc/src/__support/FPUtil/dyadic_float.h
@@ -270,11 +270,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);
@@ -301,7 +301,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;
@@ -317,7 +317,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 59886ca44be12..05a47791deed8 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 aa0069d821d0d..2e33c37dd90f0 100644
--- a/libc/src/math/generic/CMakeLists.txt
+++ b/libc/src/math/generic/CMakeLists.txt
@@ -135,6 +135,22 @@ add_header_library(
libc.src.__support.common
)
+add_header_library(
+ range_reduction_double
+ HDRS
+ range_reduction_double.h
+ range_reduction_double_fma.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 +162,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 +192,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.h b/libc/src/math/generic/range_reduction_double.h
new file mode 100644
index 0000000000000..3cd76d722da1a
--- /dev/null
+++ b/libc/src/math/generic/range_reduction_double.h
@@ -0,0 +1,67 @@
+//===-- 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_H
+#define LLVM_LIBC_SRC_MATH_GENERIC_RANGE_REDUCTION_DOUBLE_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 {
+
+using fputil::DoubleDouble;
+
+LIBC_INLINE constexpr int FAST_PASS_EXPONENT = 23;
+
+// Digits of pi/128, generated by Sollya with:
+// > a = round(pi/128, D, RN);
+// > b = round(pi/128 - a, D, RN);
+LIBC_INLINE constexpr DoubleDouble PI_OVER_128 = {0x1.1a62633145c07p-60,
+ 0x1.921fb54442d18p-6};
+
+// 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);
+// The precisions of the parts are chosen so that:
+// 1) k * a, k * b, k * c are exact in double precision
+// 2) k * b + fractional part of (k * a) is exact in double precsion
+LIBC_INLINE constexpr double MPI_OVER_128[4] = {
+ -0x1.921fb5p-6, -0x1.110b48p-32, +0x1.ee59dap-56, -0x1.98a2e03707345p-83};
+
+LIBC_INLINE constexpr double ONE_TWENTY_EIGHT_OVER_PI_D = 0x1.45f306dc9c883p5;
+
+namespace generic {
+
+LIBC_INLINE int range_reduction_small(double x, DoubleDouble &u) {
+ double prod_hi = x * ONE_TWENTY_EIGHT_OVER_PI_D;
+ double kd = fputil::nearest_integer(prod_hi);
+ int k = static_cast<int>(kd);
+
+ // x - k * (pi/128)
+ double c = fputil::multiply_add(kd, MPI_OVER_128[0], x); // Exact
+ double y_hi = fputil::multiply_add(kd, MPI_OVER_128[1], c); // Exact
+ double y_lo = fputil::multiply_add(kd, MPI_OVER_128[2], kd * MPI_OVER_128[3]);
+ u = fputil::exact_add(y_hi, y_lo);
+
+ return k;
+}
+
+// TODO: Implement generic's range_reduction_large correctly rounded for all
+// rounding modes. The current fma's range_reduction_large only works for
+// round-to-nearest without FMA instruction.
+
+} // namespace generic
+
+} // namespace LIBC_NAMESPACE
+
+#endif // LLVM_LIBC_SRC_MATH_GENERIC_RANGE_REDUCTION_DOUBLE_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 0000000000000..dafcefe4ef72e
--- /dev/null
+++ b/libc/src/math/generic/range_reduction_double_fma.h
@@ -0,0 +1,334 @@
+//===-- 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/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 fma {
+
+using fputil::DoubleDouble;
+using Float128 = fputil::DyadicFloat<128>;
+
+LIBC_INLINE constexpr int FAST_PASS_EXPONENT = 32;
+
+// Digits of pi/128, generated by Sollya with:
+// > a = round(pi/128, D, RN);
+// > b = round(pi/128 - a, D, RN);
+LIBC_INLINE constexpr DoubleDouble PI_OVER_128 = {0x1.1a62633145c07p-60,
+ 0x1.921fb54442d18p-6};
+LIBC_INLINE constexpr Float128 PI_OVER_128_F128 = {
+ Sign::POS, -133, 0xc90f'daa2'2168'c234'c4c6'628b'80dc'1cd1_u128};
+
+// Digits of 2^(16*i) / pi, generated by Sollya with:
+// 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},
+};
+
+LIBC_INLINE int range_reduction_small(double x, DoubleDouble &u) {
+ double prod_hi = x * ONE_TWENTY_EIGHT_OVER_PI[3][0];
+ double kd = fputil::nearest_integer(prod_hi);
+ int k = static_cast<int>(static_cast<int64_t>(kd));
+
+ // 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.hi, x); // Exact
+ // u_hi + u_lo ~ (y_hi + kd*(-PI_OVER_128[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.lo, y_hi);
+ u.lo = y_hi - u.hi; // Exact;
+ u.lo = fputil::multiply_add(kd, -PI_OVER_128.lo, u.lo);
+
+ return k;
+}
+
+// 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.
+//
+// Note: this algorithm works correctly without FMA instruction for the default
+// rounding mode, round-to-nearest. The limitation is due to Veltkamp's
+// Splitting algorithm used by exact_mult: double x double -> double-double.
+LIBC_INLINE int range_reduction_large(double x, DoubleDouble &u) {
+ // |x| >= 2^32.
+ using FPBits = typename fputil::FPBits<double>;
+ FPBits xbits(x);
+
+ int x_e_m62 = xbits.get_biased_exponent() - (FPBits::EXP_BIAS + 62);
+ int idx = (x_e_m62 >> 4) + 3;
+ // Scale x down by 2^(-(16 * (idx - 2))
+ xbits.set_biased_exponent((x_e_m62 & 15) + FPBits::EXP_BIAS + 62);
+ double 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 = u.hi + u.lo
+ double y_hi = ph_lo_frac + pm_hi_frac; // Exact
+ // 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);
+ u = fputil::quick_mult(y, PI_OVER_128);
+ int k = static_cast<int>(kh) + static_cast<int>(km);
+
+ return k;
+}
+
+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);
+ Float128 u = fputil::quick_mul(y, PI_OVER_128_F128);
+
+ return u;
+}
+
+// Maybe not redo-ing most of the computation, instead getting
+// y_hi, idx, pm.lo, x_reduced from range_reduction_large.
+LIBC_INLINE Float128 range_reduction_large_f128(double x) {
+ // |x| >= 2^32.
+ using FPBits = typename fputil::FPBits<double>;
+ FPBits xbits(x);
+
+ int x_e_m62 = xbits.get_biased_exponent() - (FPBits::EXP_BIAS + 62);
+ int idx = (x_e_m62 >> 4) + 3;
+ // Scale x down by 2^(-(16 * (idx - 2))
+ xbits.set_biased_exponent((x_e_m62 & 15) + FPBits::EXP_BIAS + 62);
+ double 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 = u.hi + u.lo
+ double y_hi = ph_lo_frac + pm_hi_frac; // Exact
+ // 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);
+
+ using fputil::quick_add;
+ Float128 y =
+ quick_add(y_hi_f128, quick_add(y_lo_2, quick_add(y_lo_1, y_lo_0)));
+ Float128 u = fputil::quick_mul(y, PI_OVER_128_F128);
+
+ return u;
+}
+
+} // namespace fma
+
+} // namespace LIBC_NAMESPACE
+
+#endif // LLVM_LIBC_SRC_MATH_GENERIC_RANGE_REDUCTION_DOUBLE_FMA_H
diff --git a/libc/src/math/generic/sin.cpp b/libc/src/math/generic/sin.cpp
new file mode 100644
index 0000000000000..74be21d23d929
--- /dev/null
+++ b/libc/src/math/generic/sin.cpp
@@ -0,0 +1,567 @@
+//===-- 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 "range_reduction_double_fma.h"
+
+#ifdef LIBC_TARGET_CPU_HAS_FMA
+using LIBC_NAMESPACE::fma::range_reduction_small;
+#else
+#include "range_reduction_double.h"
+using LIBC_NAMESPACE::generic::range_reduction_small;
+#endif // LIBC_TARGET_CPU_HAS_FMA
+
+// TODO: Implement generic's range_reduction_large correctly rounded for all
+// rounding modes. The current fma's range_reduction_large only works for
+// round-to-nearest without FMA instruction.
+using LIBC_NAMESPACE::fma::range_reduction_large;
+using LIBC_NAMESPACE::fma::range_reduction_large_f128;
+using LIBC_NAMESPACE::fma::range_reduction_small_f128;
+
+#include "sincos_eval.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 {
+
+// 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},
+};
+
+#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
+constexpr double ERR = 0x1.0p-67;
+#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;
+ int k;
+
+#ifdef LIBC_TARGET_CPU_HAS_FMA
+ constexpr int SMALL_EXPONENT = 32;
+#else
+ constexpr int SMALL_EXPONENT = 23;
+#endif
+
+ if (LIBC_LIKELY(x_e < FPBits::EXP_BIAS + SMALL_EXPONENT)) {
+ // |x| < 2^32
+ 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
+ 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();
+#endif // LIBC_TARGET_CPU_HAS_FMA
+ }
+
+ // // Small range reduction.
+ k = range_reduction_small(x, y);
+ } else {
+ if (LIBC_UNLIKELY(x_e > 2 * FPBits::EXP_BIAS)) {
+ // Inf or 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(x, y);
+ }
+
+ 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 = [](int kk) -> DoubleDouble {
+ // int 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 + SMALL_EXPONENT))
+ u_f128 = range_reduction_small_f128(x);
+ else
+ u_f128 = range_reduction_large_f128(x);
+
+ 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 = [](int kk) -> Float128 {
+ int 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));
+
+ 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 0000000000000..d5db18f04a8f1
--- /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 cd129e3eefb75..882181b33b9f8 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 2c7b8aa6e8c83..0000000000000
--- 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 0171b79810d4e..8ffc06b112070 100644
--- a/libc/test/src/math/sin_test.cpp
+++ b/libc/test/src/math/sin_test.cpp
@@ -12,22 +12,102 @@
#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 int N = 50;
+ 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,
+ };
+ for (int i = 0; i < N; ++i) {
+ double x = INPUTS[i];
+ EXPECT_MPFR_MATCH(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);
+
+ // TODO: Enable these tests when non-FMA versions works correctly for other
+ // rounding modes.
+ // 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 a67d0437592d5..5c8495d486ece 100644
--- a/libc/test/src/math/smoke/CMakeLists.txt
+++ b/libc/test/src/math/smoke/CMakeLists.txt
@@ -3576,3 +3576,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 0000000000000..16ced68709ca7
--- /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));
+}
>From 6210cca06379277f8e1f696412a3a6efaea31df6 Mon Sep 17 00:00:00 2001
From: Tue Ly <lntue.h at gmail.com>
Date: Mon, 17 Jun 2024 21:24:57 +0000
Subject: [PATCH 2/3] Fix small inputs for non-FMA targets.
---
libc/src/math/generic/sin.cpp | 13 ++++++++-----
1 file changed, 8 insertions(+), 5 deletions(-)
diff --git a/libc/src/math/generic/sin.cpp b/libc/src/math/generic/sin.cpp
index 74be21d23d929..dbcb5ca07c07c 100644
--- a/libc/src/math/generic/sin.cpp
+++ b/libc/src/math/generic/sin.cpp
@@ -420,11 +420,14 @@ LLVM_LIBC_FUNCTION(double, sin, (double x)) {
#ifdef LIBC_TARGET_CPU_HAS_FMA
return fputil::multiply_add(x, -0x1.0p-54, x);
#else
- 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();
+ 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
}
>From 696f27a44f46e6766095b94ccd1adbd11cb338a0 Mon Sep 17 00:00:00 2001
From: Tue Ly <lntue.h at gmail.com>
Date: Wed, 19 Jun 2024 13:43:30 +0000
Subject: [PATCH 3/3] Address comments.
---
.../src/math/generic/range_reduction_double.h | 45 +++++++-------
.../math/generic/range_reduction_double_fma.h | 59 ++++++++++++-------
libc/src/math/generic/sin.cpp | 23 ++++----
3 files changed, 75 insertions(+), 52 deletions(-)
diff --git a/libc/src/math/generic/range_reduction_double.h b/libc/src/math/generic/range_reduction_double.h
index 3cd76d722da1a..6c0e28cf7a612 100644
--- a/libc/src/math/generic/range_reduction_double.h
+++ b/libc/src/math/generic/range_reduction_double.h
@@ -21,39 +21,42 @@ using fputil::DoubleDouble;
LIBC_INLINE constexpr int FAST_PASS_EXPONENT = 23;
-// Digits of pi/128, generated by Sollya with:
-// > a = round(pi/128, D, RN);
-// > b = round(pi/128 - a, D, RN);
-LIBC_INLINE constexpr DoubleDouble PI_OVER_128 = {0x1.1a62633145c07p-60,
- 0x1.921fb54442d18p-6};
+namespace generic {
// 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);
+// > 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 + fractional part of (k * a) is exact in double precsion
+// 2) k * b + (x - (k * a)) is exact in double precsion
LIBC_INLINE constexpr double MPI_OVER_128[4] = {
-0x1.921fb5p-6, -0x1.110b48p-32, +0x1.ee59dap-56, -0x1.98a2e03707345p-83};
-LIBC_INLINE constexpr double ONE_TWENTY_EIGHT_OVER_PI_D = 0x1.45f306dc9c883p5;
-
-namespace generic {
+LIBC_INLINE unsigned range_reduction_small(double x, DoubleDouble &u) {
+ constexpr double ONE_TWENTY_EIGHT_OVER_PI = 0x1.45f306dc9c883p5;
-LIBC_INLINE int range_reduction_small(double x, DoubleDouble &u) {
- double prod_hi = x * ONE_TWENTY_EIGHT_OVER_PI_D;
+ double prod_hi = x * ONE_TWENTY_EIGHT_OVER_PI;
double kd = fputil::nearest_integer(prod_hi);
- int k = static_cast<int>(kd);
- // x - k * (pi/128)
- double c = fputil::multiply_add(kd, MPI_OVER_128[0], x); // Exact
- double y_hi = fputil::multiply_add(kd, MPI_OVER_128[1], c); // Exact
+ // 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);
-
- return k;
+ // 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));
}
// TODO: Implement generic's range_reduction_large correctly rounded for all
diff --git a/libc/src/math/generic/range_reduction_double_fma.h b/libc/src/math/generic/range_reduction_double_fma.h
index dafcefe4ef72e..0b9ded717962f 100644
--- a/libc/src/math/generic/range_reduction_double_fma.h
+++ b/libc/src/math/generic/range_reduction_double_fma.h
@@ -29,8 +29,8 @@ LIBC_INLINE constexpr int FAST_PASS_EXPONENT = 32;
// Digits of pi/128, generated by Sollya with:
// > a = round(pi/128, D, RN);
// > b = round(pi/128 - a, D, RN);
-LIBC_INLINE constexpr DoubleDouble PI_OVER_128 = {0x1.1a62633145c07p-60,
- 0x1.921fb54442d18p-6};
+LIBC_INLINE constexpr DoubleDouble PI_OVER_128_DD = {0x1.1a62633145c07p-60,
+ 0x1.921fb54442d18p-6};
LIBC_INLINE constexpr Float128 PI_OVER_128_F128 = {
Sign::POS, -133, 0xc90f'daa2'2168'c234'c4c6'628b'80dc'1cd1_u128};
@@ -194,24 +194,33 @@ LIBC_INLINE constexpr double ONE_TWENTY_EIGHT_OVER_PI[64][4] = {
-0x1.ca8bdea7f33eep-164},
};
-LIBC_INLINE int range_reduction_small(double x, DoubleDouble &u) {
+// 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) {
double prod_hi = x * ONE_TWENTY_EIGHT_OVER_PI[3][0];
double kd = fputil::nearest_integer(prod_hi);
- int k = static_cast<int>(static_cast<int64_t>(kd));
// 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.hi, x); // Exact
- // u_hi + u_lo ~ (y_hi + kd*(-PI_OVER_128[1]))
+ 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.lo, y_hi);
+ 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.lo, u.lo);
-
- return k;
+ 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
+ return static_cast<unsigned>(static_cast<int64_t>(kd));
}
// For large range |x| >= 2^32, we use the exponent of x to find 3 double-chunks
@@ -234,15 +243,15 @@ LIBC_INLINE int range_reduction_small(double x, DoubleDouble &u) {
// Note: this algorithm works correctly without FMA instruction for the default
// rounding mode, round-to-nearest. The limitation is due to Veltkamp's
// Splitting algorithm used by exact_mult: double x double -> double-double.
-LIBC_INLINE int range_reduction_large(double x, DoubleDouble &u) {
- // |x| >= 2^32.
+LIBC_INLINE unsigned range_reduction_large(double x, DoubleDouble &u) {
using FPBits = typename fputil::FPBits<double>;
FPBits xbits(x);
int x_e_m62 = xbits.get_biased_exponent() - (FPBits::EXP_BIAS + 62);
int idx = (x_e_m62 >> 4) + 3;
- // Scale x down by 2^(-(16 * (idx - 2))
+ // 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
double x_reduced = xbits.get_val();
// x * c_hi = ph.hi + ph.lo exactly.
DoubleDouble ph =
@@ -261,10 +270,20 @@ LIBC_INLINE int range_reduction_large(double x, DoubleDouble &u) {
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);
- u = fputil::quick_mult(y, PI_OVER_128);
- int k = static_cast<int>(kh) + static_cast<int>(km);
+ // 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])
+ // <= 2 * 2^77 * 2^-103 * 2^-52
+ // = 2^-77.
+ // Hence,
+ // | {x mod pi/128} - (u.hi + u.lo) | < 2 * 2^-6 * 2^-77.
+ // = 2^-82.
+ u = fputil::quick_mult(y, PI_OVER_128_DD);
- return k;
+ return static_cast<unsigned>(static_cast<int>(kh) + static_cast<int>(km));
}
LIBC_INLINE Float128 range_reduction_small_f128(double x) {
@@ -282,12 +301,11 @@ LIBC_INLINE Float128 range_reduction_small_f128(double x) {
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);
- Float128 u = fputil::quick_mul(y, PI_OVER_128_F128);
- return u;
+ return fputil::quick_mul(y, PI_OVER_128_F128);
}
-// Maybe not redo-ing most of the computation, instead getting
+// TODO: Maybe not redo-ing most of the computation, instead getting
// y_hi, idx, pm.lo, x_reduced from range_reduction_large.
LIBC_INLINE Float128 range_reduction_large_f128(double x) {
// |x| >= 2^32.
@@ -322,9 +340,8 @@ LIBC_INLINE Float128 range_reduction_large_f128(double x) {
using fputil::quick_add;
Float128 y =
quick_add(y_hi_f128, quick_add(y_lo_2, quick_add(y_lo_1, y_lo_0)));
- Float128 u = fputil::quick_mul(y, PI_OVER_128_F128);
- return u;
+ return fputil::quick_mul(y, PI_OVER_128_F128);
}
} // namespace fma
diff --git a/libc/src/math/generic/sin.cpp b/libc/src/math/generic/sin.cpp
index dbcb5ca07c07c..801db03c32c55 100644
--- a/libc/src/math/generic/sin.cpp
+++ b/libc/src/math/generic/sin.cpp
@@ -19,6 +19,7 @@
#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"
#include "range_reduction_double_fma.h"
@@ -36,8 +37,6 @@ using LIBC_NAMESPACE::fma::range_reduction_large;
using LIBC_NAMESPACE::fma::range_reduction_large_f128;
using LIBC_NAMESPACE::fma::range_reduction_small_f128;
-#include "sincos_eval.h"
-
#if ((LIBC_MATH & LIBC_MATH_SKIP_ACCURATE_PASS) != 0)
#define LIBC_MATH_SIN_SKIP_ACCURATE_PASS
#endif
@@ -388,8 +387,10 @@ LIBC_INLINE constexpr Float128 SIN_K_PI_OVER_128_F128[65] = {
#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-67;
#endif // LIBC_TARGET_CPU_HAS_FMA
+
#endif // !LIBC_MATH_SIN_SKIP_ACCURATE_PASS
} // anonymous namespace
@@ -401,7 +402,7 @@ LLVM_LIBC_FUNCTION(double, sin, (double x)) {
uint16_t x_e = xbits.get_biased_exponent();
DoubleDouble y;
- int k;
+ unsigned k;
#ifdef LIBC_TARGET_CPU_HAS_FMA
constexpr int SMALL_EXPONENT = 32;
@@ -409,8 +410,9 @@ LLVM_LIBC_FUNCTION(double, sin, (double x)) {
constexpr int SMALL_EXPONENT = 23;
#endif
+ // |x| < 2^32 (with FMA) or |x| < 2^23 (w/o FMA)
if (LIBC_LIKELY(x_e < FPBits::EXP_BIAS + SMALL_EXPONENT)) {
- // |x| < 2^32
+ // |x| < 2^-26
if (LIBC_UNLIKELY(x_e < FPBits::EXP_BIAS - 26)) {
// Signed zeros.
if (LIBC_UNLIKELY(x == 0.0))
@@ -434,8 +436,9 @@ LLVM_LIBC_FUNCTION(double, sin, (double x)) {
// // Small range reduction.
k = range_reduction_small(x, y);
} else {
+ // Inf or NaN
if (LIBC_UNLIKELY(x_e > 2 * FPBits::EXP_BIAS)) {
- // Inf or NaN
+ // sin(+-Inf) = NaN
if (xbits.get_mantissa() == 0) {
fputil::set_errno_if_required(EDOM);
fputil::raise_except_if_required(FE_INVALID);
@@ -443,7 +446,7 @@ LLVM_LIBC_FUNCTION(double, sin, (double x)) {
return x + FPBits::quiet_nan().get_val();
}
- // // Large range reduction.
+ // Large range reduction.
k = range_reduction_large(x, y);
}
@@ -465,8 +468,8 @@ LLVM_LIBC_FUNCTION(double, sin, (double x)) {
// cos_k.lo = FPBits(FPBits(cos_k.hi).uintval() ^ cos_s).get_val();
// Use 64-entry table instead:
- // auto get_idx_dd = [](int kk) -> DoubleDouble {
- // int idx = (kk & 64) ? 64 - (kk & 63) : (kk & 63);
+ // 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;
@@ -546,8 +549,8 @@ LLVM_LIBC_FUNCTION(double, sin, (double x)) {
COS_COEFFS[2], COS_COEFFS[3], COS_COEFFS[4],
COS_COEFFS[5], COS_COEFFS[6]);
- auto get_sin_k = [](int kk) -> Float128 {
- int idx = (kk & 64) ? 64 - (kk & 63) : (kk & 63);
+ 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;
More information about the libc-commits
mailing list