[libc-commits] [libc] b47ea96 - [libc][math] correct the output of `asinpif` and `acospi` (#185544)
via libc-commits
libc-commits at lists.llvm.org
Mon Mar 9 21:54:24 PDT 2026
Author: Mohamed Emad
Date: 2026-03-10T00:54:20-04:00
New Revision: b47ea968bbd40125e0e6387eb2e9c8ca37acb7b8
URL: https://github.com/llvm/llvm-project/commit/b47ea968bbd40125e0e6387eb2e9c8ca37acb7b8
DIFF: https://github.com/llvm/llvm-project/commit/b47ea968bbd40125e0e6387eb2e9c8ca37acb7b8.diff
LOG: [libc][math] correct the output of `asinpif` and `acospi` (#185544)
Currently, we have accuracy issues and some points fail in the asinpif
exhaustive test. This change fixes it by increasing the degree of the
used polynomial
```
-- Testing for FE_TONEAREST in range [0x0, 0x7f800000) --
Failed to match Func(x) against LIBC_NAMESPACE::testing::mpfr::get_mpfr_matcher<Op>( x, Func(x), 0.5, rounding).
Match value not within tolerance value of MPFR result:
Input decimal: 0.00000011348398487598387873731553554534912109375000
Input bits: 0x33F3B47B = (S: 0, E: 0x0067, M: 0x0073B47B)
Match decimal: 0.00000003612307253320068411994725465774536132812500
Match bits: 0x331B25BD = (S: 0, E: 0x0066, M: 0x001B25BD)
MPFR result: 0.00000003612307608591436292044818401336669921875000
MPFR rounded: 0x331B25BE = (S: 0, E: 0x0066, M: 0x001B25BE)
ULP error: 1.00000000000000000000000000000000000000000000000000
Test failed for 1 inputs in range: 871366656 to 872415232 [0x33f00000, 0x34000000), [0x1.ep-24, 0x1p-23)
Match value not within tolerance value of MPFR result:
Input decimal: 0.00000356030955117603298276662826538085937500000000
Input bits: 0x366EEDA7 = (S: 0, E: 0x006C, M: 0x006EEDA7)
Match decimal: 0.00000113328167117288103327155113220214843750000000
Match bits: 0x35981B43 = (S: 0, E: 0x006B, M: 0x00181B43)
MPFR result: 0.00000113328178485971875488758087158203125000000000
MPFR rounded: 0x35981B44 = (S: 0, E: 0x006B, M: 0x00181B44)
ULP error: 1.00000000000000000000000000000000000000000000000000
Test failed for 1 inputs in range: 912261120 to 913309696 [0x36600000, 0x36700000), [0x1.cp-19, 0x1.ep-19)
Match value not within tolerance value of MPFR result:
Input decimal: -0.00000011348398487598387873731553554534912109375000
Input bits: 0xB3F3B47B = (S: 1, E: 0x0067, M: 0x0073B47B)
Match decimal: -0.00000003612307253320068411994725465774536132812500
Match bits: 0xB31B25BD = (S: 1, E: 0x0066, M: 0x001B25BD)
MPFR result: -0.00000003612307608591436292044818401336669921875000
MPFR rounded: 0xB31B25BE = (S: 1, E: 0x0066, M: 0x001B25BE)
ULP error: 1.00000000000000000000000000000000000000000000000000
Test failed for 1 inputs in range: 3018850304 to 3019898880 [0xb3f00000, 0xb4000000), [-0x1.ep-24, -0x1p-23)
Match value not within tolerance value of MPFR result:
Input decimal: -0.00000356030955117603298276662826538085937500000000
Input bits: 0xB66EEDA7 = (S: 1, E: 0x006C, M: 0x006EEDA7)
Match decimal: -0.00000113328167117288103327155113220214843750000000
Match bits: 0xB5981B43 = (S: 1, E: 0x006B, M: 0x00181B43)
MPFR result: -0.00000113328178485971875488758087158203125000000000
MPFR rounded: 0xB5981B44 = (S: 1, E: 0x006B, M: 0x00181B44)
ULP error: 1.00000000000000000000000000000000000000000000000000
Test failed for 1 inputs in range: 3059744768 to 3060793344 [0xb6600000, 0xb6700000), [-0x1.cp-19, -0x1.ep-19)
```
Added:
Modified:
libc/src/__support/math/acospif.h
libc/src/__support/math/inv_trigf_utils.h
libc/test/src/math/exhaustive/asinpif_test.cpp
libc/test/src/math/smoke/acospif_test.cpp
Removed:
################################################################################
diff --git a/libc/src/__support/math/acospif.h b/libc/src/__support/math/acospif.h
index 15781f037f233..acf2b9748975b 100644
--- a/libc/src/__support/math/acospif.h
+++ b/libc/src/__support/math/acospif.h
@@ -30,18 +30,20 @@ LIBC_INLINE float acospif(float x) {
auto signed_result = [is_neg](auto r) -> auto { return is_neg ? -r : r; };
- if (LIBC_UNLIKELY(x_abs > 1.0)) {
+ if (LIBC_UNLIKELY(x_abs >= 1.0)) {
if (xbits.is_nan()) {
if (xbits.is_signaling_nan()) {
fputil::raise_except_if_required(FE_INVALID);
return FPBits::quiet_nan().get_val();
}
return x;
+ } else if (LIBC_UNLIKELY(x_abs == 1.0)) {
+ return is_neg ? 1.0f : 0.0f;
+ } else {
+ fputil::raise_except_if_required(FE_INVALID);
+ fputil::set_errno_if_required(EDOM);
+ return FPBits::quiet_nan().get_val();
}
-
- fputil::raise_except_if_required(FE_INVALID);
- fputil::set_errno_if_required(EDOM);
- return FPBits::quiet_nan().get_val();
}
// acospif(x) = 1/2 - asinpif(x)
diff --git a/libc/src/__support/math/inv_trigf_utils.h b/libc/src/__support/math/inv_trigf_utils.h
index 54a94f572cda3..83fe9f7f6201c 100644
--- a/libc/src/__support/math/inv_trigf_utils.h
+++ b/libc/src/__support/math/inv_trigf_utils.h
@@ -188,58 +188,63 @@ LIBC_INLINE double asin_eval(double xsq) {
// > [|D...|], [0, 0.5]);
// > for i from 0 to degree(P) do coeff(P, i);
// > print("Error:", dirtyinfnorm(P - g, [1e-30; 0.25]));
-// Error : 0x1.a53f84eafa3ea69bb81b6c52b3278872083fca2c757bd778acp-54 ~= 2^-54
+// Error : 0x1.6b01ec54170565911f924eb53361de37df00d74e2a10a21d5p-56 ~ 2^−55.496
//
// Non-zero coefficients (even powers only):
-LIBC_INLINE_VAR constexpr double ASINPI_COEFFS[12] = {
- 0x1.45f306dc9c881p-2, // x^0
- 0x1.b2995e7b7e756p-5, // x^2
- 0x1.8723a1d12f828p-6, // x^4
- 0x1.d1a45564b9545p-7, // x^6
- 0x1.3ce4ceaa0e1e9p-7, // x^8
- 0x1.d2c305898ea13p-8, // x^10
- 0x1.692212e27a5f9p-8, // x^12
- 0x1.2b22cc744d25bp-8, // x^14
- 0x1.8427b864479ffp-9, // x^16
- 0x1.815522d7a2bf1p-8, // x^18
- -0x1.f6df98438aef4p-9, // x^20
- 0x1.4b50c2eb13708p-7 // x^22
+constexpr double ASINPI_COEFFS[13] = {
+ 0x1.45f306dc9c883p-2, // x^0
+ 0x1.b2995e7b7af0fp-5, // x^2
+ 0x1.8723a1d61d2e9p-6, // x^4
+ 0x1.d1a4529a30a69p-7, // x^6
+ 0x1.3ce53861f8f1fp-7, // x^8
+ 0x1.d2b076c914efep-8, // x^10
+ 0x1.6a2b36f9aed68p-8, // x^12
+ 0x1.21604ae2879a2p-8, // x^14
+ 0x1.ff0549b4fd0d6p-9, // x^16
+ 0x1.035d343508f72p-9, // x^18
+ 0x1.a7b91f72b1592p-8, // x^20
+ -0x1.6a3fb073e97aep-8, // x^22
+ 0x1.547a51d51664ap-7 // x^24
};
// Evaluates P1(v2) = c1 + c2*v2 + c3*v2^2 + ... + c12*v2^11 (tail of P
// without c0) using Estrin's scheme for instruction-level parallelism.
//
-// The tail polynomial has 12 coefficients ASINPI_COEFFS[1..11] in powers of
-// v2:
-// P1(v2) = c1 + c2*v2 + c3*v2^2 + c4*v2^3 + ... + c11*v2^10
+// The tail polynomial has 12 coefficients (ASINPI_COEFFS[1] through
+// ASINPI_COEFFS[12]) in powers of v2:
+// P1(v2) = c1 + c2*v2 + c3*v2^2 + c4*v2^3 + ... + c12*v2^11
//
// Estrin pairs them bottom-up:
// Level 0 (6 pairs, using v2):
// p0 = c1 + c2*v2 p1 = c3 + c4*v2
// p2 = c5 + c6*v2 p3 = c7 + c8*v2
-// p4 = c9 + c10*v2 p5 = c11
+// p4 = c9 + c10*v2 p5 = c11 + c12*v2
// Level 1 (3 pairs, using v4):
// q0 = p0 + p1*v4 q1 = p2 + p3*v4
// q2 = p4 + p5*v4
// Level 2 (using v8):
// r0 = q0 + q1*v8 r1 = q2
-// result = q0 + q1*v8 + q1*v16
+// Level 3 (using v16):
+// result = r0 + r1*v16
LIBC_INLINE double asinpi_eval(double v2) {
double v4 = v2 * v2;
double v8 = v4 * v4;
+ double v16 = v8 * v8;
double p0 = fputil::multiply_add(v2, ASINPI_COEFFS[2], ASINPI_COEFFS[1]);
double p1 = fputil::multiply_add(v2, ASINPI_COEFFS[4], ASINPI_COEFFS[3]);
double p2 = fputil::multiply_add(v2, ASINPI_COEFFS[6], ASINPI_COEFFS[5]);
double p3 = fputil::multiply_add(v2, ASINPI_COEFFS[8], ASINPI_COEFFS[7]);
double p4 = fputil::multiply_add(v2, ASINPI_COEFFS[10], ASINPI_COEFFS[9]);
- double p5 = ASINPI_COEFFS[11];
+ double p5 = fputil::multiply_add(v2, ASINPI_COEFFS[12], ASINPI_COEFFS[11]);
double q0 = fputil::multiply_add(v4, p1, p0);
double q1 = fputil::multiply_add(v4, p3, p2);
double q2 = fputil::multiply_add(v4, p5, p4);
- return fputil::polyeval(v8, q0, q1, q2);
+ double r0 = fputil::multiply_add(v8, q1, q0);
+
+ return fputil::multiply_add(v16, q2, r0);
}
} // namespace inv_trigf_utils_internal
diff --git a/libc/test/src/math/exhaustive/asinpif_test.cpp b/libc/test/src/math/exhaustive/asinpif_test.cpp
index b327e350a0563..f4688b546da11 100644
--- a/libc/test/src/math/exhaustive/asinpif_test.cpp
+++ b/libc/test/src/math/exhaustive/asinpif_test.cpp
@@ -12,7 +12,7 @@
namespace mpfr = LIBC_NAMESPACE::testing::mpfr;
-using LlvmLibcAsinfExhaustiveTest =
+using LlvmLibcAsinpifExhaustiveTest =
LlvmLibcUnaryOpExhaustiveMathTest<float, mpfr::Operation::Asinpi,
LIBC_NAMESPACE::asinpif>;
@@ -20,7 +20,7 @@ using LlvmLibcAsinfExhaustiveTest =
static constexpr uint32_t POS_START = 0x0000'0000U;
static constexpr uint32_t POS_STOP = 0x7f80'0000U;
-TEST_F(LlvmLibcAsinfExhaustiveTest, PostiveRange) {
+TEST_F(LlvmLibcAsinpifExhaustiveTest, PostiveRange) {
test_full_range_all_roundings(POS_START, POS_STOP);
}
@@ -28,6 +28,6 @@ TEST_F(LlvmLibcAsinfExhaustiveTest, PostiveRange) {
static constexpr uint32_t NEG_START = 0xb000'0000U;
static constexpr uint32_t NEG_STOP = 0xff80'0000U;
-TEST_F(LlvmLibcAsinfExhaustiveTest, NegativeRange) {
+TEST_F(LlvmLibcAsinpifExhaustiveTest, NegativeRange) {
test_full_range_all_roundings(NEG_START, NEG_STOP);
}
diff --git a/libc/test/src/math/smoke/acospif_test.cpp b/libc/test/src/math/smoke/acospif_test.cpp
index ccee9a25de5e3..c3795764c6b8d 100644
--- a/libc/test/src/math/smoke/acospif_test.cpp
+++ b/libc/test/src/math/smoke/acospif_test.cpp
@@ -31,10 +31,10 @@ TEST_F(LlvmLibcAcospifTest, SpecialNumbers) {
EXPECT_FP_EQ(0.5f, LIBC_NAMESPACE::acospif(-0.0f));
EXPECT_MATH_ERRNO(0);
// acospif(1) = 0
- EXPECT_FP_EQ(0.0f, LIBC_NAMESPACE::acospif(1.0f));
+ EXPECT_FP_EQ_ALL_ROUNDING(0.0f, LIBC_NAMESPACE::acospif(1.0f));
EXPECT_MATH_ERRNO(0);
// acospif(-1) = 1
- EXPECT_FP_EQ(1.0f, LIBC_NAMESPACE::acospif(-1.0f));
+ EXPECT_FP_EQ_ALL_ROUNDING(1.0f, LIBC_NAMESPACE::acospif(-1.0f));
EXPECT_MATH_ERRNO(0);
EXPECT_FP_EQ_ALL_ROUNDING(aNaN, LIBC_NAMESPACE::acospif(inf));
More information about the libc-commits
mailing list