[libc-commits] [libc] [llvm] [libc][math] Refactor sincosf implementation to header only (PR #177523)
via libc-commits
libc-commits at lists.llvm.org
Sat Jan 24 18:38:37 PST 2026
https://github.com/tnuha updated https://github.com/llvm/llvm-project/pull/177523
>From 2f23252a330b838f5ee4c565f61428313be5a730 Mon Sep 17 00:00:00 2001
From: ash <ashkernel02 at gmail.com>
Date: Fri, 23 Jan 2026 08:36:39 -0700
Subject: [PATCH 1/7] [libc][math] Refactor sincosf implementation to header
only in src/__support/math directory. Part of #147386
in preparation for: https://discourse.llvm.org/t/rfc-make-clang-builtin-math-functions-constexpr-with-llvm-libc-to-support-c-23-constexpr-math-functions/86450
---
libc/shared/math.h | 1 +
libc/shared/math/sincosf.h | 23 ++
libc/src/__support/math/CMakeLists.txt | 22 ++
libc/src/__support/math/sincosf.h | 223 ++++++++++++++++++
libc/src/__support/math/sincosf_utils.h | 17 +-
libc/src/math/generic/sincosf.cpp | 198 +---------------
.../llvm-project-overlay/libc/BUILD.bazel | 7 +
7 files changed, 287 insertions(+), 204 deletions(-)
create mode 100644 libc/shared/math/sincosf.h
create mode 100644 libc/src/__support/math/sincosf.h
diff --git a/libc/shared/math.h b/libc/shared/math.h
index 6b66c5e75e457..6615fc5edfe8b 100644
--- a/libc/shared/math.h
+++ b/libc/shared/math.h
@@ -78,5 +78,6 @@
#include "math/rsqrtf.h"
#include "math/rsqrtf16.h"
#include "math/sin.h"
+#include "math/sincosf.h"
#endif // LLVM_LIBC_SHARED_MATH_H
diff --git a/libc/shared/math/sincosf.h b/libc/shared/math/sincosf.h
new file mode 100644
index 0000000000000..6996e10028874
--- /dev/null
+++ b/libc/shared/math/sincosf.h
@@ -0,0 +1,23 @@
+//===-- Shared sincosf function ---------------------------------*- 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_SHARED_MATH_SINCOSF_H
+#define LLVM_LIBC_SHARED_MATH_SINCOSF_H
+
+#include "shared/libc_common.h"
+#include "src/__support/math/sincosf.h"
+
+namespace LIBC_NAMESPACE_DECL {
+namespace shared {
+
+using math::sincosf;
+
+} // namespace shared
+} // namespace LIBC_NAMESPACE_DECL
+
+#endif // LLVM_LIBC_SHARED_MATH_SINCOSF_H
diff --git a/libc/src/__support/math/CMakeLists.txt b/libc/src/__support/math/CMakeLists.txt
index 0e7a082434fb5..32d21f53dcafd 100644
--- a/libc/src/__support/math/CMakeLists.txt
+++ b/libc/src/__support/math/CMakeLists.txt
@@ -1065,6 +1065,28 @@ add_header_library(
libc.src.__support.integer_literals
)
+add_header_library(
+ sincosf
+ HDRS
+ sincosf.h
+ DEPENDS
+ .sincosf_utils
+ libc.src.__support.math.sincosf_utils
+ libc.src.__support.FPUtil.FEnvImpl
+ libc.src.__support.FPUtil.FPBits
+ libc.src.__support.FPUtil.multiply_add
+ libc.src.__support.FPUtil.rounding_mode
+ libc.src.__support.common
+ libc.src.__support.macros.config
+ libc.src.__support.macros.optimization
+ libc.src.__support.macros.properties.cpu_features
+ libc.src.__support.math.sincosf_utils
+
+ libc.src.__support.macros.properties.types
+
+ libc.src.__support.macros.properties.cpu_features
+)
+
add_header_library(
sincosf_utils
HDRS
diff --git a/libc/src/__support/math/sincosf.h b/libc/src/__support/math/sincosf.h
new file mode 100644
index 0000000000000..12d8ad8c1555c
--- /dev/null
+++ b/libc/src/__support/math/sincosf.h
@@ -0,0 +1,223 @@
+//===-- Implementation header for sincosf ------------------------*- 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___SUPPORT_MATH_SINCOSF_H
+#define LLVM_LIBC_SRC___SUPPORT_MATH_SINCOSF_H
+
+#include "sincosf_utils.h"
+#include "src/__support/FPUtil/FEnvImpl.h"
+#include "src/__support/FPUtil/FPBits.h"
+#include "src/__support/FPUtil/multiply_add.h"
+#include "src/__support/FPUtil/rounding_mode.h"
+#include "src/__support/common.h"
+#include "src/__support/macros/config.h"
+#include "src/__support/macros/optimization.h" // LIBC_UNLIKELY
+#include "src/__support/macros/properties/cpu_features.h" // LIBC_TARGET_CPU_HAS_FMA
+#include "src/__support/math/sincosf_utils.h"
+
+namespace LIBC_NAMESPACE_DECL {
+
+namespace math {
+
+#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
+// Exceptional values
+static constexpr int N_EXCEPTS = 6;
+
+static constexpr uint32_t EXCEPT_INPUTS[N_EXCEPTS] = {
+ 0x46199998, // x = 0x1.33333p13 x
+ 0x55325019, // x = 0x1.64a032p43 x
+ 0x5922aa80, // x = 0x1.4555p51 x
+ 0x5f18b878, // x = 0x1.3170fp63 x
+ 0x6115cb11, // x = 0x1.2b9622p67 x
+ 0x7beef5ef, // x = 0x1.ddebdep120 x
+};
+
+static constexpr uint32_t EXCEPT_OUTPUTS_SIN[N_EXCEPTS][4] = {
+ {0xbeb1fa5d, 0, 1, 0}, // x = 0x1.33333p13, sin(x) = -0x1.63f4bap-2 (RZ)
+ {0xbf171adf, 0, 1, 1}, // x = 0x1.64a032p43, sin(x) = -0x1.2e35bep-1 (RZ)
+ {0xbf587521, 0, 1, 1}, // x = 0x1.4555p51, sin(x) = -0x1.b0ea42p-1 (RZ)
+ {0x3dad60f6, 1, 0, 1}, // x = 0x1.3170fp63, sin(x) = 0x1.5ac1ecp-4 (RZ)
+ {0xbe7cc1e0, 0, 1, 1}, // x = 0x1.2b9622p67, sin(x) = -0x1.f983cp-3 (RZ)
+ {0xbf587d1b, 0, 1, 1}, // x = 0x1.ddebdep120, sin(x) = -0x1.b0fa36p-1 (RZ)
+};
+
+static constexpr uint32_t EXCEPT_OUTPUTS_COS[N_EXCEPTS][4] = {
+ {0xbf70090b, 0, 1, 0}, // x = 0x1.33333p13, cos(x) = -0x1.e01216p-1 (RZ)
+ {0x3f4ea5d2, 1, 0, 0}, // x = 0x1.64a032p43, cos(x) = 0x1.9d4ba4p-1 (RZ)
+ {0x3f08aebe, 1, 0, 1}, // x = 0x1.4555p51, cos(x) = 0x1.115d7cp-1 (RZ)
+ {0x3f7f14bb, 1, 0, 0}, // x = 0x1.3170fp63, cos(x) = 0x1.fe2976p-1 (RZ)
+ {0x3f78142e, 1, 0, 1}, // x = 0x1.2b9622p67, cos(x) = 0x1.f0285cp-1 (RZ)
+ {0x3f08a21c, 1, 0, 0}, // x = 0x1.ddebdep120, cos(x) = 0x1.114438p-1 (RZ)
+};
+#endif // !LIBC_MATH_HAS_SKIP_ACCURATE_PASS
+
+LIBC_INLINE static constexpr void sincosf(float x, float *sinp, float *cosp) {
+ using FPBits = typename fputil::FPBits<float>;
+ FPBits xbits(x);
+
+ uint32_t x_abs = xbits.uintval() & 0x7fff'ffffU;
+ double xd = static_cast<double>(x);
+
+ // Range reduction:
+ // For |x| >= 2^-12, we perform range reduction as follows:
+ // Find k and y such that:
+ // x = (k + y) * pi/32
+ // k is an integer
+ // |y| < 0.5
+ // For small range (|x| < 2^45 when FMA instructions are available, 2^22
+ // otherwise), this is done by performing:
+ // k = round(x * 32/pi)
+ // y = x * 32/pi - k
+ // For large range, we will omit all the higher parts of 32/pi such that the
+ // least significant bits of their full products with x are larger than 63,
+ // since:
+ // sin((k + y + 64*i) * pi/32) = sin(x + i * 2pi) = sin(x), and
+ // cos((k + y + 64*i) * pi/32) = cos(x + i * 2pi) = cos(x).
+ //
+ // When FMA instructions are not available, we store the digits of 32/pi in
+ // chunks of 28-bit precision. This will make sure that the products:
+ // x * THIRTYTWO_OVER_PI_28[i] are all exact.
+ // When FMA instructions are available, we simply store the digits of326/pi in
+ // chunks of doubles (53-bit of precision).
+ // So when multiplying by the largest values of single precision, the
+ // resulting output should be correct up to 2^(-208 + 128) ~ 2^-80. By the
+ // worst-case analysis of range reduction, |y| >= 2^-38, so this should give
+ // us more than 40 bits of accuracy. For the worst-case estimation of range
+ // reduction, see for instances:
+ // Elementary Functions by J-M. Muller, Chapter 11,
+ // Handbook of Floating-Point Arithmetic by J-M. Muller et. al.,
+ // Chapter 10.2.
+ //
+ // Once k and y are computed, we then deduce the answer by the sine and cosine
+ // of sum formulas:
+ // sin(x) = sin((k + y)*pi/32)
+ // = sin(y*pi/32) * cos(k*pi/32) + cos(y*pi/32) * sin(k*pi/32)
+ // cos(x) = cos((k + y)*pi/32)
+ // = cos(y*pi/32) * cos(k*pi/32) - sin(y*pi/32) * sin(k*pi/32)
+ // The values of sin(k*pi/32) and cos(k*pi/32) for k = 0..63 are precomputed
+ // and stored using a vector of 32 doubles. Sin(y*pi/32) and cos(y*pi/32) are
+ // computed using degree-7 and degree-6 minimax polynomials generated by
+ // Sollya respectively.
+
+ // |x| < 0x1.0p-12f
+ if (LIBC_UNLIKELY(x_abs < 0x3980'0000U)) {
+ if (LIBC_UNLIKELY(x_abs == 0U)) {
+ // For signed zeros.
+ *sinp = x;
+ *cosp = 1.0f;
+ return;
+ }
+ // When |x| < 2^-12, the relative errors of the approximations
+ // sin(x) ~ x, cos(x) ~ 1
+ // are:
+ // |sin(x) - x| / |sin(x)| < |x^3| / (6|x|)
+ // = x^2 / 6
+ // < 2^-25
+ // < epsilon(1)/2.
+ // |cos(x) - 1| < |x^2 / 2| = 2^-25 < epsilon(1)/2.
+ // So the correctly rounded values of sin(x) and cos(x) are:
+ // sin(x) = x - sign(x)*eps(x) if rounding mode = FE_TOWARDZERO,
+ // or (rounding mode = FE_UPWARD and x is
+ // negative),
+ // = x otherwise.
+ // cos(x) = 1 - eps(x) if rounding mode = FE_TOWARDZERO or FE_DOWWARD,
+ // = 1 otherwise.
+ // To simplify the rounding decision and make it more efficient and to
+ // prevent compiler to perform constant folding, we use
+ // sin(x) = fma(x, -2^-25, x),
+ // cos(x) = fma(x*0.5f, -x, 1)
+ // instead.
+ // Note: to use the formula x - 2^-25*x to decide the correct rounding, we
+ // do need fma(x, -2^-25, x) to prevent underflow caused by -2^-25*x when
+ // |x| < 2^-125. For targets without FMA instructions, we simply use
+ // double for intermediate results as it is more efficient than using an
+ // emulated version of FMA.
+#if defined(LIBC_TARGET_CPU_HAS_FMA_FLOAT)
+ *sinp = fputil::multiply_add(x, -0x1.0p-25f, x);
+ *cosp = fputil::multiply_add(FPBits(x_abs).get_val(), -0x1.0p-25f, 1.0f);
+#else
+ *sinp = static_cast<float>(fputil::multiply_add(xd, -0x1.0p-25, xd));
+ *cosp = static_cast<float>(fputil::multiply_add(
+ static_cast<double>(FPBits(x_abs).get_val()), -0x1.0p-25, 1.0));
+#endif // LIBC_TARGET_CPU_HAS_FMA_FLOAT
+ return;
+ }
+
+ // x is inf or nan.
+ if (LIBC_UNLIKELY(x_abs >= 0x7f80'0000U)) {
+ if (xbits.is_signaling_nan()) {
+ fputil::raise_except_if_required(FE_INVALID);
+ *sinp = *cosp = FPBits::quiet_nan().get_val();
+ return;
+ }
+
+ if (x_abs == 0x7f80'0000U) {
+ fputil::set_errno_if_required(EDOM);
+ fputil::raise_except_if_required(FE_INVALID);
+ }
+ *sinp = FPBits::quiet_nan().get_val();
+ *cosp = *sinp;
+ return;
+ }
+
+#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
+ // Check exceptional values.
+ for (int i = 0; i < N_EXCEPTS; ++i) {
+ if (LIBC_UNLIKELY(x_abs == EXCEPT_INPUTS[i])) {
+ uint32_t s = EXCEPT_OUTPUTS_SIN[i][0]; // FE_TOWARDZERO
+ uint32_t c = EXCEPT_OUTPUTS_COS[i][0]; // FE_TOWARDZERO
+ bool x_sign = x < 0;
+ switch (fputil::quick_get_round()) {
+ case FE_UPWARD:
+ s += x_sign ? EXCEPT_OUTPUTS_SIN[i][2] : EXCEPT_OUTPUTS_SIN[i][1];
+ c += EXCEPT_OUTPUTS_COS[i][1];
+ break;
+ case FE_DOWNWARD:
+ s += x_sign ? EXCEPT_OUTPUTS_SIN[i][1] : EXCEPT_OUTPUTS_SIN[i][2];
+ c += EXCEPT_OUTPUTS_COS[i][2];
+ break;
+ case FE_TONEAREST:
+ s += EXCEPT_OUTPUTS_SIN[i][3];
+ c += EXCEPT_OUTPUTS_COS[i][3];
+ break;
+ }
+ *sinp = x_sign ? -FPBits(s).get_val() : FPBits(s).get_val();
+ *cosp = FPBits(c).get_val();
+
+ return;
+ }
+ }
+#endif // !LIBC_MATH_HAS_SKIP_ACCURATE_PASS
+
+ // Combine the results with the sine and cosine of sum formulas:
+ // sin(x) = sin((k + y)*pi/32)
+ // = sin(y*pi/32) * cos(k*pi/32) + cos(y*pi/32) * sin(k*pi/32)
+ // = sin_y * cos_k + (1 + cosm1_y) * sin_k
+ // = sin_y * cos_k + (cosm1_y * sin_k + sin_k)
+ // cos(x) = cos((k + y)*pi/32)
+ // = cos(y*pi/32) * cos(k*pi/32) - sin(y*pi/32) * sin(k*pi/32)
+ // = cosm1_y * cos_k + sin_y * sin_k
+ // = (cosm1_y * cos_k + cos_k) + sin_y * sin_k
+ double sin_k = 0;
+ double cos_k = 0;
+ double sin_y = 0;
+ double cosm1_y = 0;
+
+ sincosf_eval(xd, x_abs, sin_k, cos_k, sin_y, cosm1_y);
+
+ *sinp = static_cast<float>(fputil::multiply_add(
+ sin_y, cos_k, fputil::multiply_add(cosm1_y, sin_k, sin_k)));
+ *cosp = static_cast<float>(fputil::multiply_add(
+ sin_y, -sin_k, fputil::multiply_add(cosm1_y, cos_k, cos_k)));
+}
+
+} // namespace math
+
+} // namespace LIBC_NAMESPACE_DECL
+
+#endif // LLVM_LIBC_SRC___SUPPORT_MATH_SINCOSF_H
diff --git a/libc/src/__support/math/sincosf_utils.h b/libc/src/__support/math/sincosf_utils.h
index ed9d9f6bda0eb..fd2f0ac27102b 100644
--- a/libc/src/__support/math/sincosf_utils.h
+++ b/libc/src/__support/math/sincosf_utils.h
@@ -6,8 +6,8 @@
//
//===----------------------------------------------------------------------===//
-#ifndef LIBC_SRC___SUPPORT_MATH_SINCOSF_UTILS_H
-#define LIBC_SRC___SUPPORT_MATH_SINCOSF_UTILS_H
+#ifndef LLVM_LIBC_SRC___SUPPORT_MATH_SINCOSF_UTILS_H
+#define LLVM_LIBC_SRC___SUPPORT_MATH_SINCOSF_UTILS_H
#include "src/__support/FPUtil/FPBits.h"
#include "src/__support/FPUtil/PolyEval.h"
@@ -35,7 +35,7 @@ namespace LIBC_NAMESPACE_DECL {
// Table is generated with Sollya as follow:
// > display = hexadecimal;
// > for k from 0 to 63 do { D(sin(k * pi/32)); };
-const double SIN_K_PI_OVER_32[64] = {
+static const double SIN_K_PI_OVER_32[64] = {
0x0.0000000000000p+0, 0x1.917a6bc29b42cp-4, 0x1.8f8b83c69a60bp-3,
0x1.294062ed59f06p-2, 0x1.87de2a6aea963p-2, 0x1.e2b5d3806f63bp-2,
0x1.1c73b39ae68c8p-1, 0x1.44cf325091dd6p-1, 0x1.6a09e667f3bcdp-1,
@@ -88,8 +88,9 @@ static LIBC_INLINE void sincosf_poly_eval(int64_t k, double y, double &sin_k,
0x1.03c1f070c2e27p-18, -0x1.55cc84bd942p-30);
}
-LIBC_INLINE void sincosf_eval(double xd, uint32_t x_abs, double &sin_k,
- double &cos_k, double &sin_y, double &cosm1_y) {
+LIBC_INLINE static void sincosf_eval(double xd, uint32_t x_abs, double &sin_k,
+ double &cos_k, double &sin_y,
+ double &cosm1_y) {
int64_t k;
double y;
@@ -113,8 +114,8 @@ static LIBC_INLINE int64_t range_reduction_sincospi(double x, double &y) {
return static_cast<int64_t>(kd);
}
-LIBC_INLINE void sincospif_eval(double xd, double &sin_k, double &cos_k,
- double &sin_y, double &cosm1_y) {
+LIBC_INLINE static void sincospif_eval(double xd, double &sin_k, double &cos_k,
+ double &sin_y, double &cosm1_y) {
double y;
int64_t k = range_reduction_sincospi(xd, y);
sincosf_poly_eval(k, y, sin_k, cos_k, sin_y, cosm1_y);
@@ -122,4 +123,4 @@ LIBC_INLINE void sincospif_eval(double xd, double &sin_k, double &cos_k,
} // namespace LIBC_NAMESPACE_DECL
-#endif // LIBC_SRC___SUPPORT_MATH_SINCOSF_UTILS_H
+#endif // LLVM_LIBC_SRC___SUPPORT_MATH_SINCOSF_UTILS_H
diff --git a/libc/src/math/generic/sincosf.cpp b/libc/src/math/generic/sincosf.cpp
index 5179c985a2240..ebe4f8da5a0d9 100644
--- a/libc/src/math/generic/sincosf.cpp
+++ b/libc/src/math/generic/sincosf.cpp
@@ -6,206 +6,12 @@
//
//===----------------------------------------------------------------------===//
+#include "src/__support/math/sincosf.h"
#include "src/math/sincosf.h"
-#include "src/__support/FPUtil/FEnvImpl.h"
-#include "src/__support/FPUtil/FPBits.h"
-#include "src/__support/FPUtil/multiply_add.h"
-#include "src/__support/FPUtil/rounding_mode.h"
-#include "src/__support/common.h"
-#include "src/__support/macros/config.h"
-#include "src/__support/macros/optimization.h" // LIBC_UNLIKELY
-#include "src/__support/macros/properties/cpu_features.h" // LIBC_TARGET_CPU_HAS_FMA
-#include "src/__support/math/sincosf_utils.h"
namespace LIBC_NAMESPACE_DECL {
-
-#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
-// Exceptional values
-static constexpr int N_EXCEPTS = 6;
-
-static constexpr uint32_t EXCEPT_INPUTS[N_EXCEPTS] = {
- 0x46199998, // x = 0x1.33333p13 x
- 0x55325019, // x = 0x1.64a032p43 x
- 0x5922aa80, // x = 0x1.4555p51 x
- 0x5f18b878, // x = 0x1.3170fp63 x
- 0x6115cb11, // x = 0x1.2b9622p67 x
- 0x7beef5ef, // x = 0x1.ddebdep120 x
-};
-
-static constexpr uint32_t EXCEPT_OUTPUTS_SIN[N_EXCEPTS][4] = {
- {0xbeb1fa5d, 0, 1, 0}, // x = 0x1.33333p13, sin(x) = -0x1.63f4bap-2 (RZ)
- {0xbf171adf, 0, 1, 1}, // x = 0x1.64a032p43, sin(x) = -0x1.2e35bep-1 (RZ)
- {0xbf587521, 0, 1, 1}, // x = 0x1.4555p51, sin(x) = -0x1.b0ea42p-1 (RZ)
- {0x3dad60f6, 1, 0, 1}, // x = 0x1.3170fp63, sin(x) = 0x1.5ac1ecp-4 (RZ)
- {0xbe7cc1e0, 0, 1, 1}, // x = 0x1.2b9622p67, sin(x) = -0x1.f983cp-3 (RZ)
- {0xbf587d1b, 0, 1, 1}, // x = 0x1.ddebdep120, sin(x) = -0x1.b0fa36p-1 (RZ)
-};
-
-static constexpr uint32_t EXCEPT_OUTPUTS_COS[N_EXCEPTS][4] = {
- {0xbf70090b, 0, 1, 0}, // x = 0x1.33333p13, cos(x) = -0x1.e01216p-1 (RZ)
- {0x3f4ea5d2, 1, 0, 0}, // x = 0x1.64a032p43, cos(x) = 0x1.9d4ba4p-1 (RZ)
- {0x3f08aebe, 1, 0, 1}, // x = 0x1.4555p51, cos(x) = 0x1.115d7cp-1 (RZ)
- {0x3f7f14bb, 1, 0, 0}, // x = 0x1.3170fp63, cos(x) = 0x1.fe2976p-1 (RZ)
- {0x3f78142e, 1, 0, 1}, // x = 0x1.2b9622p67, cos(x) = 0x1.f0285cp-1 (RZ)
- {0x3f08a21c, 1, 0, 0}, // x = 0x1.ddebdep120, cos(x) = 0x1.114438p-1 (RZ)
-};
-#endif // !LIBC_MATH_HAS_SKIP_ACCURATE_PASS
-
LLVM_LIBC_FUNCTION(void, sincosf, (float x, float *sinp, float *cosp)) {
- using FPBits = typename fputil::FPBits<float>;
- FPBits xbits(x);
-
- uint32_t x_abs = xbits.uintval() & 0x7fff'ffffU;
- double xd = static_cast<double>(x);
-
- // Range reduction:
- // For |x| >= 2^-12, we perform range reduction as follows:
- // Find k and y such that:
- // x = (k + y) * pi/32
- // k is an integer
- // |y| < 0.5
- // For small range (|x| < 2^45 when FMA instructions are available, 2^22
- // otherwise), this is done by performing:
- // k = round(x * 32/pi)
- // y = x * 32/pi - k
- // For large range, we will omit all the higher parts of 32/pi such that the
- // least significant bits of their full products with x are larger than 63,
- // since:
- // sin((k + y + 64*i) * pi/32) = sin(x + i * 2pi) = sin(x), and
- // cos((k + y + 64*i) * pi/32) = cos(x + i * 2pi) = cos(x).
- //
- // When FMA instructions are not available, we store the digits of 32/pi in
- // chunks of 28-bit precision. This will make sure that the products:
- // x * THIRTYTWO_OVER_PI_28[i] are all exact.
- // When FMA instructions are available, we simply store the digits of326/pi in
- // chunks of doubles (53-bit of precision).
- // So when multiplying by the largest values of single precision, the
- // resulting output should be correct up to 2^(-208 + 128) ~ 2^-80. By the
- // worst-case analysis of range reduction, |y| >= 2^-38, so this should give
- // us more than 40 bits of accuracy. For the worst-case estimation of range
- // reduction, see for instances:
- // Elementary Functions by J-M. Muller, Chapter 11,
- // Handbook of Floating-Point Arithmetic by J-M. Muller et. al.,
- // Chapter 10.2.
- //
- // Once k and y are computed, we then deduce the answer by the sine and cosine
- // of sum formulas:
- // sin(x) = sin((k + y)*pi/32)
- // = sin(y*pi/32) * cos(k*pi/32) + cos(y*pi/32) * sin(k*pi/32)
- // cos(x) = cos((k + y)*pi/32)
- // = cos(y*pi/32) * cos(k*pi/32) - sin(y*pi/32) * sin(k*pi/32)
- // The values of sin(k*pi/32) and cos(k*pi/32) for k = 0..63 are precomputed
- // and stored using a vector of 32 doubles. Sin(y*pi/32) and cos(y*pi/32) are
- // computed using degree-7 and degree-6 minimax polynomials generated by
- // Sollya respectively.
-
- // |x| < 0x1.0p-12f
- if (LIBC_UNLIKELY(x_abs < 0x3980'0000U)) {
- if (LIBC_UNLIKELY(x_abs == 0U)) {
- // For signed zeros.
- *sinp = x;
- *cosp = 1.0f;
- return;
- }
- // When |x| < 2^-12, the relative errors of the approximations
- // sin(x) ~ x, cos(x) ~ 1
- // are:
- // |sin(x) - x| / |sin(x)| < |x^3| / (6|x|)
- // = x^2 / 6
- // < 2^-25
- // < epsilon(1)/2.
- // |cos(x) - 1| < |x^2 / 2| = 2^-25 < epsilon(1)/2.
- // So the correctly rounded values of sin(x) and cos(x) are:
- // sin(x) = x - sign(x)*eps(x) if rounding mode = FE_TOWARDZERO,
- // or (rounding mode = FE_UPWARD and x is
- // negative),
- // = x otherwise.
- // cos(x) = 1 - eps(x) if rounding mode = FE_TOWARDZERO or FE_DOWWARD,
- // = 1 otherwise.
- // To simplify the rounding decision and make it more efficient and to
- // prevent compiler to perform constant folding, we use
- // sin(x) = fma(x, -2^-25, x),
- // cos(x) = fma(x*0.5f, -x, 1)
- // instead.
- // Note: to use the formula x - 2^-25*x to decide the correct rounding, we
- // do need fma(x, -2^-25, x) to prevent underflow caused by -2^-25*x when
- // |x| < 2^-125. For targets without FMA instructions, we simply use
- // double for intermediate results as it is more efficient than using an
- // emulated version of FMA.
-#if defined(LIBC_TARGET_CPU_HAS_FMA_FLOAT)
- *sinp = fputil::multiply_add(x, -0x1.0p-25f, x);
- *cosp = fputil::multiply_add(FPBits(x_abs).get_val(), -0x1.0p-25f, 1.0f);
-#else
- *sinp = static_cast<float>(fputil::multiply_add(xd, -0x1.0p-25, xd));
- *cosp = static_cast<float>(fputil::multiply_add(
- static_cast<double>(FPBits(x_abs).get_val()), -0x1.0p-25, 1.0));
-#endif // LIBC_TARGET_CPU_HAS_FMA_FLOAT
- return;
- }
-
- // x is inf or nan.
- if (LIBC_UNLIKELY(x_abs >= 0x7f80'0000U)) {
- if (xbits.is_signaling_nan()) {
- fputil::raise_except_if_required(FE_INVALID);
- *sinp = *cosp = FPBits::quiet_nan().get_val();
- return;
- }
-
- if (x_abs == 0x7f80'0000U) {
- fputil::set_errno_if_required(EDOM);
- fputil::raise_except_if_required(FE_INVALID);
- }
- *sinp = FPBits::quiet_nan().get_val();
- *cosp = *sinp;
- return;
- }
-
-#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
- // Check exceptional values.
- for (int i = 0; i < N_EXCEPTS; ++i) {
- if (LIBC_UNLIKELY(x_abs == EXCEPT_INPUTS[i])) {
- uint32_t s = EXCEPT_OUTPUTS_SIN[i][0]; // FE_TOWARDZERO
- uint32_t c = EXCEPT_OUTPUTS_COS[i][0]; // FE_TOWARDZERO
- bool x_sign = x < 0;
- switch (fputil::quick_get_round()) {
- case FE_UPWARD:
- s += x_sign ? EXCEPT_OUTPUTS_SIN[i][2] : EXCEPT_OUTPUTS_SIN[i][1];
- c += EXCEPT_OUTPUTS_COS[i][1];
- break;
- case FE_DOWNWARD:
- s += x_sign ? EXCEPT_OUTPUTS_SIN[i][1] : EXCEPT_OUTPUTS_SIN[i][2];
- c += EXCEPT_OUTPUTS_COS[i][2];
- break;
- case FE_TONEAREST:
- s += EXCEPT_OUTPUTS_SIN[i][3];
- c += EXCEPT_OUTPUTS_COS[i][3];
- break;
- }
- *sinp = x_sign ? -FPBits(s).get_val() : FPBits(s).get_val();
- *cosp = FPBits(c).get_val();
-
- return;
- }
- }
-#endif // !LIBC_MATH_HAS_SKIP_ACCURATE_PASS
-
- // Combine the results with the sine and cosine of sum formulas:
- // sin(x) = sin((k + y)*pi/32)
- // = sin(y*pi/32) * cos(k*pi/32) + cos(y*pi/32) * sin(k*pi/32)
- // = sin_y * cos_k + (1 + cosm1_y) * sin_k
- // = sin_y * cos_k + (cosm1_y * sin_k + sin_k)
- // cos(x) = cos((k + y)*pi/32)
- // = cos(y*pi/32) * cos(k*pi/32) - sin(y*pi/32) * sin(k*pi/32)
- // = cosm1_y * cos_k + sin_y * sin_k
- // = (cosm1_y * cos_k + cos_k) + sin_y * sin_k
- double sin_k, cos_k, sin_y, cosm1_y;
-
- sincosf_eval(xd, x_abs, sin_k, cos_k, sin_y, cosm1_y);
-
- *sinp = static_cast<float>(fputil::multiply_add(
- sin_y, cos_k, fputil::multiply_add(cosm1_y, sin_k, sin_k)));
- *cosp = static_cast<float>(fputil::multiply_add(
- sin_y, -sin_k, fputil::multiply_add(cosm1_y, cos_k, cos_k)));
+ return math::sincosf(x, sinp, cosp);
}
} // namespace LIBC_NAMESPACE_DECL
diff --git a/utils/bazel/llvm-project-overlay/libc/BUILD.bazel b/utils/bazel/llvm-project-overlay/libc/BUILD.bazel
index 7251f491bc757..3101ccad5f4d4 100644
--- a/utils/bazel/llvm-project-overlay/libc/BUILD.bazel
+++ b/utils/bazel/llvm-project-overlay/libc/BUILD.bazel
@@ -4902,6 +4902,13 @@ libc_math_function(
],
)
+libc_math_function(
+ name = "sincosf",
+ additional_deps = [
+ ":__support_math_sincosf",
+ ],
+)
+
libc_math_function(
name = "sinf",
additional_deps = [
>From c676b612fbf4694f5111c5247d142b43f91ba0b2 Mon Sep 17 00:00:00 2001
From: Ash <97464181+tnuha at users.noreply.github.com>
Date: Fri, 23 Jan 2026 22:14:06 +0000
Subject: [PATCH 2/7] Update libc/src/math/generic/sincosf.cpp
Co-authored-by: Muhammad Bassiouni <60100307+bassiounix at users.noreply.github.com>
---
libc/src/math/generic/sincosf.cpp | 2 +-
1 file changed, 1 insertion(+), 1 deletion(-)
diff --git a/libc/src/math/generic/sincosf.cpp b/libc/src/math/generic/sincosf.cpp
index ebe4f8da5a0d9..f2881cfed7a80 100644
--- a/libc/src/math/generic/sincosf.cpp
+++ b/libc/src/math/generic/sincosf.cpp
@@ -6,8 +6,8 @@
//
//===----------------------------------------------------------------------===//
-#include "src/__support/math/sincosf.h"
#include "src/math/sincosf.h"
+#include "src/__support/math/sincosf.h"
namespace LIBC_NAMESPACE_DECL {
LLVM_LIBC_FUNCTION(void, sincosf, (float x, float *sinp, float *cosp)) {
>From 2e13d16172b75c7ffba9c0fb5ae94cb7c50920ce Mon Sep 17 00:00:00 2001
From: Ash <97464181+tnuha at users.noreply.github.com>
Date: Fri, 23 Jan 2026 22:14:17 +0000
Subject: [PATCH 3/7] Update libc/src/__support/math/sincosf_utils.h
Co-authored-by: Muhammad Bassiouni <60100307+bassiounix at users.noreply.github.com>
---
libc/src/__support/math/sincosf_utils.h | 2 +-
1 file changed, 1 insertion(+), 1 deletion(-)
diff --git a/libc/src/__support/math/sincosf_utils.h b/libc/src/__support/math/sincosf_utils.h
index fd2f0ac27102b..265998fb0ec7a 100644
--- a/libc/src/__support/math/sincosf_utils.h
+++ b/libc/src/__support/math/sincosf_utils.h
@@ -35,7 +35,7 @@ namespace LIBC_NAMESPACE_DECL {
// Table is generated with Sollya as follow:
// > display = hexadecimal;
// > for k from 0 to 63 do { D(sin(k * pi/32)); };
-static const double SIN_K_PI_OVER_32[64] = {
+LIBC_INLINE_VAR const double SIN_K_PI_OVER_32[64] = {
0x0.0000000000000p+0, 0x1.917a6bc29b42cp-4, 0x1.8f8b83c69a60bp-3,
0x1.294062ed59f06p-2, 0x1.87de2a6aea963p-2, 0x1.e2b5d3806f63bp-2,
0x1.1c73b39ae68c8p-1, 0x1.44cf325091dd6p-1, 0x1.6a09e667f3bcdp-1,
>From 3a6a62ba78debeffeec677b0692fe18c1cff6359 Mon Sep 17 00:00:00 2001
From: Ash <97464181+tnuha at users.noreply.github.com>
Date: Fri, 23 Jan 2026 22:14:29 +0000
Subject: [PATCH 4/7] Update libc/src/__support/math/CMakeLists.txt
Co-authored-by: Muhammad Bassiouni <60100307+bassiounix at users.noreply.github.com>
---
libc/src/__support/math/CMakeLists.txt | 1 -
1 file changed, 1 deletion(-)
diff --git a/libc/src/__support/math/CMakeLists.txt b/libc/src/__support/math/CMakeLists.txt
index 32d21f53dcafd..a30a466f86007 100644
--- a/libc/src/__support/math/CMakeLists.txt
+++ b/libc/src/__support/math/CMakeLists.txt
@@ -1080,7 +1080,6 @@ add_header_library(
libc.src.__support.macros.config
libc.src.__support.macros.optimization
libc.src.__support.macros.properties.cpu_features
- libc.src.__support.math.sincosf_utils
libc.src.__support.macros.properties.types
>From 087114842195b1eb633ce7eb58b7f793ad2c97f8 Mon Sep 17 00:00:00 2001
From: Ash <97464181+tnuha at users.noreply.github.com>
Date: Fri, 23 Jan 2026 22:14:45 +0000
Subject: [PATCH 5/7] Update libc/src/__support/math/sincosf.h
Co-authored-by: Muhammad Bassiouni <60100307+bassiounix at users.noreply.github.com>
---
libc/src/__support/math/sincosf.h | 12 ++++++++----
1 file changed, 8 insertions(+), 4 deletions(-)
diff --git a/libc/src/__support/math/sincosf.h b/libc/src/__support/math/sincosf.h
index 12d8ad8c1555c..3f681bec09822 100644
--- a/libc/src/__support/math/sincosf.h
+++ b/libc/src/__support/math/sincosf.h
@@ -24,11 +24,13 @@ namespace LIBC_NAMESPACE_DECL {
namespace math {
+namespace sincosf_internal {
+
#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
// Exceptional values
-static constexpr int N_EXCEPTS = 6;
+LIBC_INLINE_VAR constexpr int N_EXCEPTS = 6;
-static constexpr uint32_t EXCEPT_INPUTS[N_EXCEPTS] = {
+LIBC_INLINE_VAR constexpr uint32_t EXCEPT_INPUTS[N_EXCEPTS] = {
0x46199998, // x = 0x1.33333p13 x
0x55325019, // x = 0x1.64a032p43 x
0x5922aa80, // x = 0x1.4555p51 x
@@ -37,7 +39,7 @@ static constexpr uint32_t EXCEPT_INPUTS[N_EXCEPTS] = {
0x7beef5ef, // x = 0x1.ddebdep120 x
};
-static constexpr uint32_t EXCEPT_OUTPUTS_SIN[N_EXCEPTS][4] = {
+LIBC_INLINE_VAR constexpr uint32_t EXCEPT_OUTPUTS_SIN[N_EXCEPTS][4] = {
{0xbeb1fa5d, 0, 1, 0}, // x = 0x1.33333p13, sin(x) = -0x1.63f4bap-2 (RZ)
{0xbf171adf, 0, 1, 1}, // x = 0x1.64a032p43, sin(x) = -0x1.2e35bep-1 (RZ)
{0xbf587521, 0, 1, 1}, // x = 0x1.4555p51, sin(x) = -0x1.b0ea42p-1 (RZ)
@@ -46,7 +48,7 @@ static constexpr uint32_t EXCEPT_OUTPUTS_SIN[N_EXCEPTS][4] = {
{0xbf587d1b, 0, 1, 1}, // x = 0x1.ddebdep120, sin(x) = -0x1.b0fa36p-1 (RZ)
};
-static constexpr uint32_t EXCEPT_OUTPUTS_COS[N_EXCEPTS][4] = {
+LIBC_INLINE_VAR constexpr uint32_t EXCEPT_OUTPUTS_COS[N_EXCEPTS][4] = {
{0xbf70090b, 0, 1, 0}, // x = 0x1.33333p13, cos(x) = -0x1.e01216p-1 (RZ)
{0x3f4ea5d2, 1, 0, 0}, // x = 0x1.64a032p43, cos(x) = 0x1.9d4ba4p-1 (RZ)
{0x3f08aebe, 1, 0, 1}, // x = 0x1.4555p51, cos(x) = 0x1.115d7cp-1 (RZ)
@@ -55,8 +57,10 @@ static constexpr uint32_t EXCEPT_OUTPUTS_COS[N_EXCEPTS][4] = {
{0x3f08a21c, 1, 0, 0}, // x = 0x1.ddebdep120, cos(x) = 0x1.114438p-1 (RZ)
};
#endif // !LIBC_MATH_HAS_SKIP_ACCURATE_PASS
+} // namespace sincosf_internal
LIBC_INLINE static constexpr void sincosf(float x, float *sinp, float *cosp) {
+ using namespace sincosf_internal;
using FPBits = typename fputil::FPBits<float>;
FPBits xbits(x);
>From 21c028dcae8507ebfac0d01dc76cd3dbf041d487 Mon Sep 17 00:00:00 2001
From: ash <ashkernel02 at gmail.com>
Date: Fri, 23 Jan 2026 16:02:37 -0700
Subject: [PATCH 6/7] update CMakeLists.txt
---
libc/src/math/generic/CMakeLists.txt | 11 +----------
1 file changed, 1 insertion(+), 10 deletions(-)
diff --git a/libc/src/math/generic/CMakeLists.txt b/libc/src/math/generic/CMakeLists.txt
index bf799ce41363c..e0f1924a1a08b 100644
--- a/libc/src/math/generic/CMakeLists.txt
+++ b/libc/src/math/generic/CMakeLists.txt
@@ -450,16 +450,7 @@ add_entrypoint_object(
HDRS
../sincosf.h
DEPENDS
- libc.src.__support.math.range_reduction
- libc.src.__support.math.sincosf_utils
- libc.src.errno.errno
- 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.polyeval
- libc.src.__support.FPUtil.rounding_mode
- libc.src.__support.macros.optimization
+ libc.src.__support.math.sincosf
)
add_entrypoint_object(
>From c614d21a0fdbda647388f9a71a9686b4e99c679b Mon Sep 17 00:00:00 2001
From: ash <ashkernel02 at gmail.com>
Date: Sat, 24 Jan 2026 19:38:46 -0700
Subject: [PATCH 7/7] update shared tests
---
libc/test/shared/CMakeLists.txt | 1 +
libc/test/shared/shared_math_test.cpp | 5 +++++
2 files changed, 6 insertions(+)
diff --git a/libc/test/shared/CMakeLists.txt b/libc/test/shared/CMakeLists.txt
index 3be8f992aa486..9d7b170cd67ed 100644
--- a/libc/test/shared/CMakeLists.txt
+++ b/libc/test/shared/CMakeLists.txt
@@ -74,4 +74,5 @@ add_fp_unittest(
libc.src.__support.math.rsqrtf
libc.src.__support.math.rsqrtf16
libc.src.__support.math.sin
+ libc.src.__support.math.sincosf
)
diff --git a/libc/test/shared/shared_math_test.cpp b/libc/test/shared/shared_math_test.cpp
index 655d96692f1b8..56258e620cd36 100644
--- a/libc/test/shared/shared_math_test.cpp
+++ b/libc/test/shared/shared_math_test.cpp
@@ -53,6 +53,7 @@ TEST(LlvmLibcSharedMathTest, AllFloat16) {
TEST(LlvmLibcSharedMathTest, AllFloat) {
int exponent;
+ float sin, cos;
EXPECT_FP_EQ(0x1.921fb6p+0, LIBC_NAMESPACE::shared::acosf(0.0f));
EXPECT_FP_EQ(0x0p+0f, LIBC_NAMESPACE::shared::acoshf(1.0f));
@@ -84,6 +85,10 @@ TEST(LlvmLibcSharedMathTest, AllFloat) {
EXPECT_EQ(long(0), LIBC_NAMESPACE::shared::llogbf(1.0f));
EXPECT_FP_EQ(0x0p+0f, LIBC_NAMESPACE::shared::logbf(1.0f));
EXPECT_FP_EQ(0x1p+0f, LIBC_NAMESPACE::shared::rsqrtf(1.0f));
+
+ LIBC_NAMESPACE::shared::sincosf(0, &sin, &cos);
+ ASSERT_FP_EQ(1.0, cos);
+ ASSERT_FP_EQ(0.0, sin);
}
TEST(LlvmLibcSharedMathTest, AllDouble) {
More information about the libc-commits
mailing list