[libc-commits] [libc] [llvm] [libc][math][c23] implement double-precision asinpi (PR #188158)
Mohamed Emad via libc-commits
libc-commits at lists.llvm.org
Thu Mar 26 22:11:48 PDT 2026
https://github.com/hulxv updated https://github.com/llvm/llvm-project/pull/188158
>From 9659dd9b4f1c41d72b7fbee1fbfb71bbc53bb105 Mon Sep 17 00:00:00 2001
From: hulxv <hulxxv at gmail.com>
Date: Tue, 24 Mar 2026 03:26:49 +0200
Subject: [PATCH 1/6] [libc][math][c23] implement double-precision asinpi
---
libc/config/baremetal/aarch64/entrypoints.txt | 1 +
libc/config/baremetal/arm/entrypoints.txt | 1 +
libc/config/baremetal/riscv/entrypoints.txt | 1 +
libc/config/darwin/aarch64/entrypoints.txt | 1 +
libc/config/gpu/amdgpu/entrypoints.txt | 1 +
libc/config/gpu/nvptx/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/config/linux/x86_64/entrypoints.txt | 1 +
libc/config/windows/entrypoints.txt | 1 +
libc/include/math.yaml | 12 +
libc/shared/math.h | 1 +
libc/shared/math/asinpi.h | 23 ++
libc/src/__support/math/CMakeLists.txt | 17 ++
libc/src/__support/math/asin_utils.h | 68 +++++
libc/src/__support/math/asinpi.h | 242 ++++++++++++++++++
libc/src/math/CMakeLists.txt | 1 +
libc/src/math/asinpi.h | 20 ++
libc/src/math/generic/CMakeLists.txt | 10 +
libc/src/math/generic/asinpi.cpp | 16 ++
libc/test/shared/CMakeLists.txt | 1 +
libc/test/shared/shared_math_test.cpp | 1 +
libc/test/src/math/CMakeLists.txt | 11 +
libc/test/src/math/asinpi_test.cpp | 95 +++++++
libc/test/src/math/smoke/CMakeLists.txt | 12 +
libc/test/src/math/smoke/asinpi_test.cpp | 52 ++++
.../llvm-project-overlay/libc/BUILD.bazel | 25 ++
28 files changed, 618 insertions(+)
create mode 100644 libc/shared/math/asinpi.h
create mode 100644 libc/src/__support/math/asinpi.h
create mode 100644 libc/src/math/asinpi.h
create mode 100644 libc/src/math/generic/asinpi.cpp
create mode 100644 libc/test/src/math/asinpi_test.cpp
create mode 100644 libc/test/src/math/smoke/asinpi_test.cpp
diff --git a/libc/config/baremetal/aarch64/entrypoints.txt b/libc/config/baremetal/aarch64/entrypoints.txt
index 278e1f5afde70..b6936653ef128 100644
--- a/libc/config/baremetal/aarch64/entrypoints.txt
+++ b/libc/config/baremetal/aarch64/entrypoints.txt
@@ -343,6 +343,7 @@ set(TARGET_LIBM_ENTRYPOINTS
libc.src.math.asin
libc.src.math.asinf
libc.src.math.asinhf
+ libc.src.math.asinpi
libc.src.math.asinpif
libc.src.math.atan2
libc.src.math.atan2f
diff --git a/libc/config/baremetal/arm/entrypoints.txt b/libc/config/baremetal/arm/entrypoints.txt
index 955dc643586e9..cbc054d212d0a 100644
--- a/libc/config/baremetal/arm/entrypoints.txt
+++ b/libc/config/baremetal/arm/entrypoints.txt
@@ -347,6 +347,7 @@ set(TARGET_LIBM_ENTRYPOINTS
libc.src.math.asin
libc.src.math.asinf
libc.src.math.asinhf
+ libc.src.math.asinpi
libc.src.math.asinpif
libc.src.math.atan2
libc.src.math.atan2f
diff --git a/libc/config/baremetal/riscv/entrypoints.txt b/libc/config/baremetal/riscv/entrypoints.txt
index 472792e48fa01..de93de5e5f2fc 100644
--- a/libc/config/baremetal/riscv/entrypoints.txt
+++ b/libc/config/baremetal/riscv/entrypoints.txt
@@ -342,6 +342,7 @@ set(TARGET_LIBM_ENTRYPOINTS
libc.src.math.acospif
libc.src.math.asinf
libc.src.math.asinhf
+ libc.src.math.asinpi
libc.src.math.asinpif
libc.src.math.atan2
libc.src.math.atan2f
diff --git a/libc/config/darwin/aarch64/entrypoints.txt b/libc/config/darwin/aarch64/entrypoints.txt
index 051ab45f450fa..885e30caf3b1f 100644
--- a/libc/config/darwin/aarch64/entrypoints.txt
+++ b/libc/config/darwin/aarch64/entrypoints.txt
@@ -158,6 +158,7 @@ set(TARGET_LIBM_ENTRYPOINTS
libc.src.math.asin
libc.src.math.asinf
libc.src.math.asinhf
+ libc.src.math.asinpi
libc.src.math.asinpif
libc.src.math.atan2
libc.src.math.atan2f
diff --git a/libc/config/gpu/amdgpu/entrypoints.txt b/libc/config/gpu/amdgpu/entrypoints.txt
index a65b6f0274fd8..24ade13e4012d 100644
--- a/libc/config/gpu/amdgpu/entrypoints.txt
+++ b/libc/config/gpu/amdgpu/entrypoints.txt
@@ -287,6 +287,7 @@ set(TARGET_LIBM_ENTRYPOINTS
libc.src.math.asin
libc.src.math.asinf
libc.src.math.asinhf
+ libc.src.math.asinpi
libc.src.math.asinpif
libc.src.math.atan
libc.src.math.atan2
diff --git a/libc/config/gpu/nvptx/entrypoints.txt b/libc/config/gpu/nvptx/entrypoints.txt
index fee0038c88cc0..e5fff18867025 100644
--- a/libc/config/gpu/nvptx/entrypoints.txt
+++ b/libc/config/gpu/nvptx/entrypoints.txt
@@ -287,6 +287,7 @@ set(TARGET_LIBM_ENTRYPOINTS
libc.src.math.asin
libc.src.math.asinf
libc.src.math.asinhf
+ libc.src.math.asinpi
libc.src.math.asinpif
libc.src.math.atan
libc.src.math.atan2
diff --git a/libc/config/linux/aarch64/entrypoints.txt b/libc/config/linux/aarch64/entrypoints.txt
index bdca0baf478c7..21f04e3b36710 100644
--- a/libc/config/linux/aarch64/entrypoints.txt
+++ b/libc/config/linux/aarch64/entrypoints.txt
@@ -428,6 +428,7 @@ set(TARGET_LIBM_ENTRYPOINTS
libc.src.math.asin
libc.src.math.asinf
libc.src.math.asinhf
+ libc.src.math.asinpi
libc.src.math.asinpif
libc.src.math.atan2
libc.src.math.atan2f
diff --git a/libc/config/linux/arm/entrypoints.txt b/libc/config/linux/arm/entrypoints.txt
index 3d3e77d826b06..36f4007a65cb7 100644
--- a/libc/config/linux/arm/entrypoints.txt
+++ b/libc/config/linux/arm/entrypoints.txt
@@ -258,6 +258,7 @@ set(TARGET_LIBM_ENTRYPOINTS
libc.src.math.asin
libc.src.math.asinf
libc.src.math.asinhf
+ libc.src.math.asinpi
libc.src.math.asinpif
libc.src.math.atan2
libc.src.math.atan2f
diff --git a/libc/config/linux/riscv/entrypoints.txt b/libc/config/linux/riscv/entrypoints.txt
index 13343dea76fe3..81271c9d17909 100644
--- a/libc/config/linux/riscv/entrypoints.txt
+++ b/libc/config/linux/riscv/entrypoints.txt
@@ -432,6 +432,7 @@ set(TARGET_LIBM_ENTRYPOINTS
libc.src.math.asin
libc.src.math.asinf
libc.src.math.asinhf
+ libc.src.math.asinpi
libc.src.math.asinpif
libc.src.math.atan2
libc.src.math.atan2f
diff --git a/libc/config/linux/x86_64/entrypoints.txt b/libc/config/linux/x86_64/entrypoints.txt
index 24df7e63e3ebf..ce35ebd6fbdef 100644
--- a/libc/config/linux/x86_64/entrypoints.txt
+++ b/libc/config/linux/x86_64/entrypoints.txt
@@ -477,6 +477,7 @@ set(TARGET_LIBM_ENTRYPOINTS
libc.src.math.asin
libc.src.math.asinf
libc.src.math.asinhf
+ libc.src.math.asinpi
libc.src.math.asinpif
libc.src.math.atan2
libc.src.math.atan2f
diff --git a/libc/config/windows/entrypoints.txt b/libc/config/windows/entrypoints.txt
index f777fc6c94bc2..24b4b67cedee3 100644
--- a/libc/config/windows/entrypoints.txt
+++ b/libc/config/windows/entrypoints.txt
@@ -137,6 +137,7 @@ set(TARGET_LIBM_ENTRYPOINTS
libc.src.math.asin
libc.src.math.asinf
libc.src.math.asinhf
+ libc.src.math.asinpi
libc.src.math.asinpif
libc.src.math.atan2
libc.src.math.atan2f
diff --git a/libc/include/math.yaml b/libc/include/math.yaml
index afd3ae33305c1..59ffcf50f86c9 100644
--- a/libc/include/math.yaml
+++ b/libc/include/math.yaml
@@ -79,6 +79,18 @@ functions:
arguments:
- type: _Float16
guard: LIBC_TYPES_HAS_FLOAT16
+ - name: asinpi
+ standards:
+ - stdc
+ return_type: double
+ arguments:
+ - type: double
+ - name: asinpif
+ standards:
+ - stdc
+ return_type: float
+ arguments:
+ - type: float
- name: asinpif16
standards:
- stdc
diff --git a/libc/shared/math.h b/libc/shared/math.h
index ede0ebd5371ac..e4762de92cbd0 100644
--- a/libc/shared/math.h
+++ b/libc/shared/math.h
@@ -23,6 +23,7 @@
#include "math/asinf16.h"
#include "math/asinhf.h"
#include "math/asinhf16.h"
+#include "math/asinpi.h"
#include "math/asinpif.h"
#include "math/asinpif16.h"
#include "math/atan.h"
diff --git a/libc/shared/math/asinpi.h b/libc/shared/math/asinpi.h
new file mode 100644
index 0000000000000..aed947ef00e76
--- /dev/null
+++ b/libc/shared/math/asinpi.h
@@ -0,0 +1,23 @@
+//===-- Shared asinpi header -------------------------------------*- 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_ASINPI_H
+#define LLVM_LIBC_SHARED_MATH_ASINPI_H
+
+#include "shared/libc_common.h"
+#include "src/__support/math/asinpi.h"
+
+namespace LIBC_NAMESPACE_DECL {
+namespace shared {
+
+using math::asinpi;
+
+} // namespace shared
+} // namespace LIBC_NAMESPACE_DECL
+
+#endif // LLVM_LIBC_SHARED_MATH_ASINPI_H
diff --git a/libc/src/__support/math/CMakeLists.txt b/libc/src/__support/math/CMakeLists.txt
index 0df8262cfb5f2..bc2f871218da2 100644
--- a/libc/src/__support/math/CMakeLists.txt
+++ b/libc/src/__support/math/CMakeLists.txt
@@ -187,6 +187,23 @@ add_header_library(
libc.src.__support.macros.optimization
)
+add_header_library(
+ asinpi
+ HDRS
+ asinpi.h
+ DEPENDS
+ .asin_utils
+ 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.multiply_add
+ libc.src.__support.FPUtil.polyeval
+ libc.src.__support.FPUtil.sqrt
+ libc.src.__support.macros.optimization
+ libc.src.__support.macros.properties.cpu_features
+)
+
add_header_library(
asinpif
HDRS
diff --git a/libc/src/__support/math/asin_utils.h b/libc/src/__support/math/asin_utils.h
index 2feb1485ec667..2991474f6236e 100644
--- a/libc/src/__support/math/asin_utils.h
+++ b/libc/src/__support/math/asin_utils.h
@@ -63,6 +63,45 @@ LIBC_INLINE double asin_eval(double u) {
return fputil::polyeval(u4, d0, d1, d2);
}
+// Coefficients for the polynomial approximation of asin(x)/(pi*x) on [0, 0.5].
+// Generated by Sollya:
+// > prec = 200;
+// > g = asin(x) / (pi * x);
+// > P = fpminimax(g, [|0, 2, 4, 6, 8, 10, 12, 14, 16, 18, 20, 22|],
+// > [|D...|], [0, 0.5]);
+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
+};
+
+// Evaluate P(u) where P(u) ~ asin(sqrt(u))/(pi*sqrt(u)), using Estrin's scheme.
+LIBC_INLINE double asinpi_eval(double u) {
+ double u2 = u * u;
+ double c0 = fputil::multiply_add(u, ASINPI_COEFFS[1], ASINPI_COEFFS[0]);
+ double c1 = fputil::multiply_add(u, ASINPI_COEFFS[3], ASINPI_COEFFS[2]);
+ double c2 = fputil::multiply_add(u, ASINPI_COEFFS[5], ASINPI_COEFFS[4]);
+ double c3 = fputil::multiply_add(u, ASINPI_COEFFS[7], ASINPI_COEFFS[6]);
+ double c4 = fputil::multiply_add(u, ASINPI_COEFFS[9], ASINPI_COEFFS[8]);
+ double c5 = fputil::multiply_add(u, ASINPI_COEFFS[11], ASINPI_COEFFS[10]);
+
+ double u4 = u2 * u2;
+ double d0 = fputil::multiply_add(u2, c1, c0);
+ double d1 = fputil::multiply_add(u2, c3, c2);
+ double d2 = fputil::multiply_add(u2, c5, c4);
+
+ return fputil::polyeval(u4, d0, d1, d2);
+}
+
#else
// The Taylor expansion of asin(x) around 0 is:
@@ -565,6 +604,35 @@ LIBC_INLINE constexpr Float128 asin_eval(const Float128 &u, unsigned idx) {
ASIN_COEFFS_F128[idx][14], ASIN_COEFFS_F128[idx][15]);
}
+// 1/pi as DoubleDouble (~106 bits of precision).
+LIBC_INLINE_VAR constexpr DoubleDouble ONE_OVER_PI_DD = {
+ -0x1.6b01ec5417056p-56, 0x1.45f306dc9c883p-2};
+
+// 1/pi as Float128 (~128 bits of precision).
+LIBC_INLINE_VAR constexpr Float128 ONE_OVER_PI_F128 = {
+ Sign::POS, -129, 0xa2f9836e'4e441529'fc2757d1'f534ddc0_u128};
+
+// 0.5 as Float128 (exact).
+LIBC_INLINE_VAR constexpr Float128 HALF_F128 = {
+ Sign::POS, -128, 0x80000000'00000000'00000000'00000000_u128};
+
+// Compute asin(sqrt(u))/(pi*sqrt(u)) in DoubleDouble precision by wrapping
+// asin_eval and multiplying by 1/pi.
+LIBC_INLINE DoubleDouble asinpi_eval(const DoubleDouble &u, unsigned &idx,
+ double &err) {
+ DoubleDouble p = asin_eval(u, idx, err);
+ DoubleDouble result = fputil::quick_mult(p, ONE_OVER_PI_DD);
+ // Scale error by ~1/pi and add multiplication rounding error.
+ err = err * 0x1.46p-2 + 0x1.0p-98;
+ return result;
+}
+
+// Compute asin(sqrt(u))/(pi*sqrt(u)) in Float128 precision.
+LIBC_INLINE constexpr Float128 asinpi_eval(const Float128 &u, unsigned idx) {
+ Float128 p = asin_eval(u, idx);
+ return fputil::quick_mul(p, ONE_OVER_PI_F128);
+}
+
#endif // LIBC_MATH_HAS_SKIP_ACCURATE_PASS
} // namespace asin_internal
diff --git a/libc/src/__support/math/asinpi.h b/libc/src/__support/math/asinpi.h
new file mode 100644
index 0000000000000..2c1b499deebaf
--- /dev/null
+++ b/libc/src/__support/math/asinpi.h
@@ -0,0 +1,242 @@
+//===-- Implementation header for asinpi -------------------------*- 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_ASINPI_H
+#define LLVM_LIBC_SRC___SUPPORT_MATH_ASINPI_H
+
+#include "asin_utils.h"
+#include "src/__support/FPUtil/FEnvImpl.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/sqrt.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/asin_utils.h"
+
+namespace LIBC_NAMESPACE_DECL {
+
+namespace math {
+
+LIBC_INLINE double asinpi(double x) {
+ using namespace asin_internal;
+ using FPBits = fputil::FPBits<double>;
+
+ FPBits xbits(x);
+ int x_exp = xbits.get_biased_exponent();
+
+ // |x| < 0.5.
+ if (x_exp < FPBits::EXP_BIAS - 1) {
+ // |x| < 2^-26.
+ if (LIBC_UNLIKELY(x_exp < FPBits::EXP_BIAS - 26)) {
+ // asinpi(+-0) = +-0.
+ if (LIBC_UNLIKELY(xbits.abs().uintval() == 0))
+ return x;
+ // When |x| < 2^-26, asinpi(x) ~ x/pi.
+ // The relative error of x/pi is:
+ // |asinpi(x) - x/pi| / |asinpi(x)| < x^2/6 < 2^-54.
+#if defined(LIBC_MATH_HAS_SKIP_ACCURATE_PASS)
+ return x * ASINPI_COEFFS[0];
+#elif defined(LIBC_TARGET_CPU_HAS_FMA_DOUBLE)
+ return fputil::multiply_add(x, ONE_OVER_PI_DD.hi,
+ x * ONE_OVER_PI_DD.lo);
+#else
+ double hi = x * ONE_OVER_PI_DD.hi;
+ return hi + x * ONE_OVER_PI_DD.lo;
+#endif // LIBC_MATH_HAS_SKIP_ACCURATE_PASS
+ }
+
+#ifdef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
+ return x * asinpi_eval(x * x);
+#else
+ using Float128 = fputil::DyadicFloat<128>;
+ using DoubleDouble = fputil::DoubleDouble;
+
+ unsigned idx = 0;
+ DoubleDouble x_sq = fputil::exact_mult(x, x);
+ double err = xbits.abs().get_val() * 0x1.0p-51;
+ // Polynomial approximation:
+ // p ~ asin(x)/(pi*x)
+
+ DoubleDouble p = asinpi_eval(x_sq, idx, err);
+ // asinpi(x) ~ x * p
+ DoubleDouble r0 = fputil::exact_mult(x, p.hi);
+ double r_lo = fputil::multiply_add(x, p.lo, r0.lo);
+
+ // Ziv's accuracy test.
+ double r_upper = r0.hi + (r_lo + err);
+ double r_lower = r0.hi + (r_lo - err);
+
+ if (LIBC_LIKELY(r_upper == r_lower))
+ return r_upper;
+
+ // Ziv's accuracy test failed, perform 128-bit calculation.
+
+ // Recalculate mod 1/64.
+ idx = static_cast<unsigned>(fputil::nearest_integer(x_sq.hi * 0x1.0p6));
+
+ Float128 x_f128(x);
+
+#ifdef LIBC_TARGET_CPU_HAS_FMA_DOUBLE
+ Float128 u_hi(
+ fputil::multiply_add(static_cast<double>(idx), -0x1.0p-6, x_sq.hi));
+ Float128 u = fputil::quick_add(u_hi, Float128(x_sq.lo));
+#else
+ Float128 x_sq_f128 = fputil::quick_mul(x_f128, x_f128);
+ Float128 u = fputil::quick_add(
+ x_sq_f128, Float128(static_cast<double>(idx) * (-0x1.0p-6)));
+#endif // LIBC_TARGET_CPU_HAS_FMA_DOUBLE
+
+ Float128 p_f128 = asinpi_eval(u, idx);
+ Float128 r = fputil::quick_mul(x_f128, p_f128);
+
+ return static_cast<double>(r);
+#endif // LIBC_MATH_HAS_SKIP_ACCURATE_PASS
+ }
+ // |x| >= 0.5
+
+ double x_abs = xbits.abs().get_val();
+
+ // Maintaining the sign:
+ constexpr double SIGN[2] = {1.0, -1.0};
+ double x_sign = SIGN[xbits.is_neg()];
+
+ // |x| >= 1
+ if (LIBC_UNLIKELY(x_exp >= FPBits::EXP_BIAS)) {
+ // x = +-1, asinpi(x) = +- 0.5
+ if (x_abs == 1.0) {
+ return x_sign * 0.5;
+ }
+ // |x| > 1, return NaN.
+ if (xbits.is_quiet_nan())
+ return x;
+
+ // Set domain error for non-NaN input.
+ if (!xbits.is_nan())
+ fputil::set_errno_if_required(EDOM);
+
+ fputil::raise_except_if_required(FE_INVALID);
+ return FPBits::quiet_nan().get_val();
+ }
+
+ // When |x| >= 0.5, we perform range reduction as follow:
+ //
+ // Assume further that 0.5 <= x < 1, and let:
+ // y = asin(x)
+ // Using the identity:
+ // asin(x) = pi/2 - 2 * asin( sqrt( (1 - x)/2 ) )
+ // We get:
+ // asinpi(x) = asin(x)/pi = 0.5 - 2 * asin(sqrt(u)) / pi
+ // = 0.5 - 2 * sqrt(u) * [asin(sqrt(u)) / (pi * sqrt(u))]
+ // = 0.5 - 2 * sqrt(u) * asinpi_eval(u)
+ // where u = (1 - |x|) / 2.
+
+ // u = (1 - |x|)/2
+ double u = fputil::multiply_add(x_abs, -0.5, 0.5);
+ // v_hi ~ sqrt(u).
+ double v_hi = fputil::sqrt<double>(u);
+
+#ifdef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
+ double p = asinpi_eval(u);
+ double r = x_sign * fputil::multiply_add(-2.0 * v_hi, p, 0.5);
+ return r;
+#else
+ using Float128 = fputil::DyadicFloat<128>;
+ using DoubleDouble = fputil::DoubleDouble;
+
+#ifdef LIBC_TARGET_CPU_HAS_FMA_DOUBLE
+ double h = fputil::multiply_add(v_hi, -v_hi, u);
+#else
+ DoubleDouble v_hi_sq = fputil::exact_mult(v_hi, v_hi);
+ double h = (u - v_hi_sq.hi) - v_hi_sq.lo;
+#endif // LIBC_TARGET_CPU_HAS_FMA_DOUBLE
+
+ // Scale v_lo and v_hi by 2 from the formula:
+ // vh = v_hi * 2
+ // vl = 2*v_lo = h / v_hi.
+ double vh = v_hi * 2.0;
+ double vl = h / v_hi;
+
+ // Polynomial approximation:
+ // p ~ asin(sqrt(u))/(pi*sqrt(u))
+ unsigned idx = 0;
+ double err = vh * 0x1.0p-51;
+
+ DoubleDouble p = asinpi_eval(DoubleDouble{0.0, u}, idx, err);
+
+ // Perform computations in double-double arithmetic:
+ // asinpi(x) = 0.5 - (vh + vl) * p
+ DoubleDouble r0 = fputil::quick_mult(DoubleDouble{vl, vh}, p);
+ DoubleDouble r = fputil::exact_add(0.5, -r0.hi);
+
+ double r_lo = -r0.lo + r.lo;
+
+ // Ziv's accuracy test.
+
+#ifdef LIBC_TARGET_CPU_HAS_FMA_DOUBLE
+ double r_upper = fputil::multiply_add(
+ r.hi, x_sign, fputil::multiply_add(r_lo, x_sign, err));
+ double r_lower = fputil::multiply_add(
+ r.hi, x_sign, fputil::multiply_add(r_lo, x_sign, -err));
+#else
+ r_lo *= x_sign;
+ r.hi *= x_sign;
+ double r_upper = r.hi + (r_lo + err);
+ double r_lower = r.hi + (r_lo - err);
+#endif // LIBC_TARGET_CPU_HAS_FMA_DOUBLE
+
+ if (LIBC_LIKELY(r_upper == r_lower))
+ return r_upper;
+
+ // Ziv's accuracy test failed, we redo the computations in Float128.
+ // Recalculate mod 1/64.
+ idx = static_cast<unsigned>(fputil::nearest_integer(u * 0x1.0p6));
+
+ // After the first step of Newton-Raphson approximating v = sqrt(u):
+ // sqrt(u) = v_hi + h / (sqrt(u) + v_hi)
+ // v_lo = h / (2 * v_hi)
+ // Add second-order correction:
+ // v_ll = -v_lo * (h / (4u))
+
+ // Get the rounding error of vl = 2 * v_lo ~ h / vh
+#ifdef LIBC_TARGET_CPU_HAS_FMA_DOUBLE
+ double vl_lo = fputil::multiply_add(-v_hi, vl, h) / v_hi;
+#else
+ DoubleDouble vh_vl = fputil::exact_mult(v_hi, vl);
+ double vl_lo = ((h - vh_vl.hi) - vh_vl.lo) / v_hi;
+#endif // LIBC_TARGET_CPU_HAS_FMA_DOUBLE
+ // vll = 2*v_ll = -vl * (h / (4u)).
+ double t = h * (-0.25) / u;
+ double vll = fputil::multiply_add(vl, t, vl_lo);
+ // m_v = -(v_hi + v_lo + v_ll).
+ Float128 m_v = fputil::quick_add(
+ Float128(vh), fputil::quick_add(Float128(vl), Float128(vll)));
+ m_v.sign = Sign::NEG;
+
+ // Perform computations in Float128:
+ // asinpi(x) = 0.5 - (v_hi + v_lo + vll) * P_pi(u).
+ Float128 y_f128(fputil::multiply_add(static_cast<double>(idx), -0x1.0p-6, u));
+
+ Float128 p_f128 = asinpi_eval(y_f128, idx);
+ Float128 r0_f128 = fputil::quick_mul(m_v, p_f128);
+ Float128 r_f128 = fputil::quick_add(HALF_F128, r0_f128);
+
+ if (xbits.is_neg())
+ r_f128.sign = Sign::NEG;
+
+ return static_cast<double>(r_f128);
+#endif // LIBC_MATH_HAS_SKIP_ACCURATE_PASS
+}
+
+} // namespace math
+
+} // namespace LIBC_NAMESPACE_DECL
+
+#endif // LLVM_LIBC_SRC___SUPPORT_MATH_ASINPI_H
diff --git a/libc/src/math/CMakeLists.txt b/libc/src/math/CMakeLists.txt
index e617950368994..232bf63d0c653 100644
--- a/libc/src/math/CMakeLists.txt
+++ b/libc/src/math/CMakeLists.txt
@@ -59,6 +59,7 @@ add_math_entrypoint_object(asinh)
add_math_entrypoint_object(asinhf)
add_math_entrypoint_object(asinhf16)
+add_math_entrypoint_object(asinpi)
add_math_entrypoint_object(asinpif)
add_math_entrypoint_object(asinpif16)
diff --git a/libc/src/math/asinpi.h b/libc/src/math/asinpi.h
new file mode 100644
index 0000000000000..086c3acfa8dfd
--- /dev/null
+++ b/libc/src/math/asinpi.h
@@ -0,0 +1,20 @@
+//===-- Implementation header for asinpi -------------------------*- 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_ASINPI_H
+#define LLVM_LIBC_SRC_MATH_ASINPI_H
+
+#include "src/__support/macros/config.h"
+
+namespace LIBC_NAMESPACE_DECL {
+
+double asinpi(double x);
+
+} // namespace LIBC_NAMESPACE_DECL
+
+#endif // LLVM_LIBC_SRC_MATH_ASINPI_H
diff --git a/libc/src/math/generic/CMakeLists.txt b/libc/src/math/generic/CMakeLists.txt
index 418cf85b84a20..38b21b6b9cc9e 100644
--- a/libc/src/math/generic/CMakeLists.txt
+++ b/libc/src/math/generic/CMakeLists.txt
@@ -3992,6 +3992,16 @@ add_entrypoint_object(
libc.src.__support.math.asinhf16
)
+add_entrypoint_object(
+ asinpi
+ SRCS
+ asinpi.cpp
+ HDRS
+ ../asinpi.h
+ DEPENDS
+ libc.src.__support.math.asinpi
+)
+
add_entrypoint_object(
asinpif
SRCS
diff --git a/libc/src/math/generic/asinpi.cpp b/libc/src/math/generic/asinpi.cpp
new file mode 100644
index 0000000000000..18bb61db6be35
--- /dev/null
+++ b/libc/src/math/generic/asinpi.cpp
@@ -0,0 +1,16 @@
+//===-- Double-precision asinpi 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/asinpi.h"
+#include "src/__support/math/asinpi.h"
+
+namespace LIBC_NAMESPACE_DECL {
+
+LLVM_LIBC_FUNCTION(double, asinpi, (double x)) { return math::asinpi(x); }
+
+} // namespace LIBC_NAMESPACE_DECL
diff --git a/libc/test/shared/CMakeLists.txt b/libc/test/shared/CMakeLists.txt
index c90e5687d8c33..8a04afa347c5a 100644
--- a/libc/test/shared/CMakeLists.txt
+++ b/libc/test/shared/CMakeLists.txt
@@ -20,6 +20,7 @@ add_fp_unittest(
libc.src.__support.math.asinf16
libc.src.__support.math.asinhf
libc.src.__support.math.asinhf16
+ libc.src.__support.math.asinpi
libc.src.__support.math.asinpif
libc.src.__support.math.asinpif16
libc.src.__support.math.atan
diff --git a/libc/test/shared/shared_math_test.cpp b/libc/test/shared/shared_math_test.cpp
index 17045ce5edfdb..b888c1eb4b8c0 100644
--- a/libc/test/shared/shared_math_test.cpp
+++ b/libc/test/shared/shared_math_test.cpp
@@ -235,6 +235,7 @@ TEST(LlvmLibcSharedMathTest, AllDouble) {
EXPECT_FP_EQ(0x1.921fb54442d18p+0, LIBC_NAMESPACE::shared::acos(0.0));
EXPECT_FP_EQ(0x0p+0, LIBC_NAMESPACE::shared::asin(0.0));
+ EXPECT_FP_EQ(0x0p+0, LIBC_NAMESPACE::shared::asinpi(0.0));
EXPECT_FP_EQ(0x0p+0, LIBC_NAMESPACE::shared::atan(0.0));
EXPECT_FP_EQ(0x0p+0, LIBC_NAMESPACE::shared::atan2(0.0, 0.0));
EXPECT_FP_EQ(0x0p+0, LIBC_NAMESPACE::shared::cbrt(0.0));
diff --git a/libc/test/src/math/CMakeLists.txt b/libc/test/src/math/CMakeLists.txt
index 73b5ebf5a856e..1922815481e2c 100644
--- a/libc/test/src/math/CMakeLists.txt
+++ b/libc/test/src/math/CMakeLists.txt
@@ -2476,6 +2476,17 @@ add_fp_unittest(
libc.src.math.asinf16
)
+add_fp_unittest(
+ asinpi_test
+ NEED_MPFR
+ SUITE
+ libc-math-unittests
+ SRCS
+ asinpi_test.cpp
+ DEPENDS
+ libc.src.math.asinpi
+)
+
add_fp_unittest(
asinpif_test
NEED_MPFR
diff --git a/libc/test/src/math/asinpi_test.cpp b/libc/test/src/math/asinpi_test.cpp
new file mode 100644
index 0000000000000..dcaba1b25a6e5
--- /dev/null
+++ b/libc/test/src/math/asinpi_test.cpp
@@ -0,0 +1,95 @@
+//===-- Unittests for asinpi ----------------------------------------------===//
+//
+// 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/__support/macros/optimization.h"
+#include "src/math/asinpi.h"
+#include "test/UnitTest/FPMatcher.h"
+#include "test/UnitTest/Test.h"
+#include "utils/MPFRWrapper/MPFRUtils.h"
+
+#ifdef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
+#define TOLERANCE 6
+#else
+#define TOLERANCE 0
+#endif // LIBC_MATH_HAS_SKIP_ACCURATE_PASS
+
+using LlvmLibcAsinpiTest = LIBC_NAMESPACE::testing::FPTest<double>;
+
+namespace mpfr = LIBC_NAMESPACE::testing::mpfr;
+
+using LIBC_NAMESPACE::testing::tlog;
+
+TEST_F(LlvmLibcAsinpiTest, InDoubleRange) {
+ constexpr uint64_t COUNT = 123'451;
+ uint64_t START = FPBits(0x1.0p-60).uintval();
+ uint64_t STOP = FPBits(1.0).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 = 0.0, 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 (FPBits(v).is_nan() || FPBits(v).is_inf())
+ continue;
+ libc_errno = 0;
+ double result = LIBC_NAMESPACE::asinpi(x);
+ ++cc;
+ if (FPBits(result).is_nan() || FPBits(result).is_inf())
+ continue;
+
+ ++count;
+
+ if (!TEST_MPFR_MATCH_ROUNDING_SILENTLY(mpfr::Operation::Asinpi, x,
+ result, TOLERANCE + 0.5,
+ rounding_mode)) {
+ ++fails;
+ while (!TEST_MPFR_MATCH_ROUNDING_SILENTLY(mpfr::Operation::Asinpi, x,
+ result, tol,
+ rounding_mode)) {
+ mx = x;
+ mr = result;
+
+ if (tol > 1000.0)
+ break;
+
+ tol *= 2.0;
+ }
+ }
+ }
+ if (fails) {
+ tlog << " Asinpi failed: " << fails << "/" << count << "/" << cc
+ << " tests.\n";
+ tlog << " Max ULPs is at most: " << static_cast<uint64_t>(tol)
+ << ".\n";
+ EXPECT_MPFR_MATCH(mpfr::Operation::Asinpi, mx, mr, 0.5, rounding_mode);
+ }
+ };
+
+ tlog << " Test Rounding To Nearest...\n";
+ test(mpfr::RoundingMode::Nearest);
+
+#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
+ 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);
+#endif // LIBC_MATH_HAS_SKIP_ACCURATE_PASS
+}
diff --git a/libc/test/src/math/smoke/CMakeLists.txt b/libc/test/src/math/smoke/CMakeLists.txt
index 3d52873c40bb8..893b992db71f9 100644
--- a/libc/test/src/math/smoke/CMakeLists.txt
+++ b/libc/test/src/math/smoke/CMakeLists.txt
@@ -4703,6 +4703,18 @@ add_fp_unittest(
libc.src.math.asinhf16
)
+add_fp_unittest(
+ asinpi_test
+ SUITE
+ libc-math-smoke-tests
+ SRCS
+ asinpi_test.cpp
+ DEPENDS
+ libc.hdr.fenv_macros
+ libc.src.math.asinpi
+ libc.src.__support.FPUtil.fp_bits
+)
+
add_fp_unittest(
asinpif_test
SUITE
diff --git a/libc/test/src/math/smoke/asinpi_test.cpp b/libc/test/src/math/smoke/asinpi_test.cpp
new file mode 100644
index 0000000000000..fcd3dd3e5ddc3
--- /dev/null
+++ b/libc/test/src/math/smoke/asinpi_test.cpp
@@ -0,0 +1,52 @@
+//===-- Unittests for asinpi ----------------------------------------------===//
+//
+// 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 "hdr/fenv_macros.h"
+#include "src/math/asinpi.h"
+#include "test/UnitTest/FPMatcher.h"
+#include "test/UnitTest/Test.h"
+
+using LlvmLibcAsinpiTest = LIBC_NAMESPACE::testing::FPTest<double>;
+
+TEST_F(LlvmLibcAsinpiTest, SpecialNumbers) {
+ EXPECT_FP_EQ_WITH_EXCEPTION_ALL_ROUNDING(aNaN, LIBC_NAMESPACE::asinpi(sNaN),
+ FE_INVALID);
+ EXPECT_FP_EQ_ALL_ROUNDING(aNaN, LIBC_NAMESPACE::asinpi(aNaN));
+ EXPECT_FP_EQ_ALL_ROUNDING(zero, LIBC_NAMESPACE::asinpi(zero));
+ EXPECT_FP_EQ_ALL_ROUNDING(neg_zero, LIBC_NAMESPACE::asinpi(neg_zero));
+ EXPECT_FP_EQ_ALL_ROUNDING(aNaN, LIBC_NAMESPACE::asinpi(inf));
+ EXPECT_FP_EQ_ALL_ROUNDING(aNaN, LIBC_NAMESPACE::asinpi(neg_inf));
+ EXPECT_FP_EQ_ALL_ROUNDING(aNaN, LIBC_NAMESPACE::asinpi(2.0));
+ EXPECT_FP_EQ_ALL_ROUNDING(aNaN, LIBC_NAMESPACE::asinpi(-2.0));
+ EXPECT_FP_EQ(0.5, LIBC_NAMESPACE::asinpi(1.0));
+ EXPECT_FP_EQ(-0.5, LIBC_NAMESPACE::asinpi(-1.0));
+}
+
+#ifdef LIBC_TEST_FTZ_DAZ
+
+using namespace LIBC_NAMESPACE::testing;
+
+TEST_F(LlvmLibcAsinpiTest, FTZMode) {
+ ModifyMXCSR mxcsr(FTZ);
+
+ EXPECT_TRUE(zero == LIBC_NAMESPACE::asinpi(min_denormal));
+}
+
+TEST_F(LlvmLibcAsinpiTest, DAZMode) {
+ ModifyMXCSR mxcsr(DAZ);
+
+ EXPECT_TRUE(zero == LIBC_NAMESPACE::asinpi(min_denormal));
+}
+
+TEST_F(LlvmLibcAsinpiTest, FTZDAZMode) {
+ ModifyMXCSR mxcsr(FTZ | DAZ);
+
+ EXPECT_TRUE(zero == LIBC_NAMESPACE::asinpi(min_denormal));
+}
+
+#endif
diff --git a/utils/bazel/llvm-project-overlay/libc/BUILD.bazel b/utils/bazel/llvm-project-overlay/libc/BUILD.bazel
index 31a229d94a187..17140ebd4249a 100644
--- a/utils/bazel/llvm-project-overlay/libc/BUILD.bazel
+++ b/utils/bazel/llvm-project-overlay/libc/BUILD.bazel
@@ -2850,6 +2850,24 @@ libc_support_library(
],
)
+libc_support_library(
+ name = "__support_math_asinpi",
+ hdrs = ["src/__support/math/asinpi.h"],
+ deps = [
+ ":__support_fputil_cast",
+ ":__support_fputil_except_value_utils",
+ ":__support_fputil_fenv_impl",
+ ":__support_fputil_fp_bits",
+ ":__support_fputil_multiply_add",
+ ":__support_fputil_polyeval",
+ ":__support_fputil_rounding_mode",
+ ":__support_fputil_sqrt",
+ ":__support_macros_config",
+ ":__support_macros_optimization",
+ ":__support_math_asin_utils",
+ ],
+)
+
libc_support_library(
name = "__support_math_asinpif",
hdrs = ["src/__support/math/asinpif.h"],
@@ -5967,6 +5985,13 @@ libc_math_function(
],
)
+libc_math_function(
+ name = "asinpi",
+ additional_deps = [
+ ":__support_math_asinpi",
+ ],
+)
+
libc_math_function(
name = "asinpif",
additional_deps = [
>From 1fb72ce26354ca42888ea2a92e2e0ae4578a15c3 Mon Sep 17 00:00:00 2001
From: hulxv <hulxxv at gmail.com>
Date: Tue, 24 Mar 2026 03:39:55 +0200
Subject: [PATCH 2/6] fix format
---
libc/shared/math/asinpi.h | 2 +-
libc/src/__support/math/asin_utils.h | 4 ++--
libc/src/__support/math/asinpi.h | 5 ++---
libc/src/math/asinpi.h | 2 +-
libc/src/math/generic/asinpi.cpp | 2 +-
libc/test/src/math/asinpi_test.cpp | 11 ++++-------
6 files changed, 11 insertions(+), 15 deletions(-)
diff --git a/libc/shared/math/asinpi.h b/libc/shared/math/asinpi.h
index aed947ef00e76..42ada49e0a1ae 100644
--- a/libc/shared/math/asinpi.h
+++ b/libc/shared/math/asinpi.h
@@ -1,4 +1,4 @@
-//===-- Shared asinpi header -------------------------------------*- C++ -*-===//
+//===-- Shared asinpi header ------------------------------------*- C++ -*-===//
//
// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
// See https://llvm.org/LICENSE.txt for license information.
diff --git a/libc/src/__support/math/asin_utils.h b/libc/src/__support/math/asin_utils.h
index 2991474f6236e..e39262f7d19a8 100644
--- a/libc/src/__support/math/asin_utils.h
+++ b/libc/src/__support/math/asin_utils.h
@@ -605,8 +605,8 @@ LIBC_INLINE constexpr Float128 asin_eval(const Float128 &u, unsigned idx) {
}
// 1/pi as DoubleDouble (~106 bits of precision).
-LIBC_INLINE_VAR constexpr DoubleDouble ONE_OVER_PI_DD = {
- -0x1.6b01ec5417056p-56, 0x1.45f306dc9c883p-2};
+LIBC_INLINE_VAR constexpr DoubleDouble ONE_OVER_PI_DD = {-0x1.6b01ec5417056p-56,
+ 0x1.45f306dc9c883p-2};
// 1/pi as Float128 (~128 bits of precision).
LIBC_INLINE_VAR constexpr Float128 ONE_OVER_PI_F128 = {
diff --git a/libc/src/__support/math/asinpi.h b/libc/src/__support/math/asinpi.h
index 2c1b499deebaf..e492a5e64ee7e 100644
--- a/libc/src/__support/math/asinpi.h
+++ b/libc/src/__support/math/asinpi.h
@@ -1,4 +1,4 @@
-//===-- Implementation header for asinpi -------------------------*- C++ -*-===//
+//===-- Implementation header for asinpi ------------------------*- C++ -*-===//
//
// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
// See https://llvm.org/LICENSE.txt for license information.
@@ -45,8 +45,7 @@ LIBC_INLINE double asinpi(double x) {
#if defined(LIBC_MATH_HAS_SKIP_ACCURATE_PASS)
return x * ASINPI_COEFFS[0];
#elif defined(LIBC_TARGET_CPU_HAS_FMA_DOUBLE)
- return fputil::multiply_add(x, ONE_OVER_PI_DD.hi,
- x * ONE_OVER_PI_DD.lo);
+ return fputil::multiply_add(x, ONE_OVER_PI_DD.hi, x * ONE_OVER_PI_DD.lo);
#else
double hi = x * ONE_OVER_PI_DD.hi;
return hi + x * ONE_OVER_PI_DD.lo;
diff --git a/libc/src/math/asinpi.h b/libc/src/math/asinpi.h
index 086c3acfa8dfd..7b8770b5e52a9 100644
--- a/libc/src/math/asinpi.h
+++ b/libc/src/math/asinpi.h
@@ -1,4 +1,4 @@
-//===-- Implementation header for asinpi -------------------------*- C++ -*-===//
+//===-- Implementation header for asinpi ------------------------*- C++ -*-===//
//
// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
// See https://llvm.org/LICENSE.txt for license information.
diff --git a/libc/src/math/generic/asinpi.cpp b/libc/src/math/generic/asinpi.cpp
index 18bb61db6be35..20b76ff1ffd6f 100644
--- a/libc/src/math/generic/asinpi.cpp
+++ b/libc/src/math/generic/asinpi.cpp
@@ -1,4 +1,4 @@
-//===-- Double-precision asinpi function -----------------------------------===//
+//===-- Double-precision asinpi function ----------------------------------===//
//
// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
// See https://llvm.org/LICENSE.txt for license information.
diff --git a/libc/test/src/math/asinpi_test.cpp b/libc/test/src/math/asinpi_test.cpp
index dcaba1b25a6e5..75dfbc6608eac 100644
--- a/libc/test/src/math/asinpi_test.cpp
+++ b/libc/test/src/math/asinpi_test.cpp
@@ -53,13 +53,11 @@ TEST_F(LlvmLibcAsinpiTest, InDoubleRange) {
++count;
- if (!TEST_MPFR_MATCH_ROUNDING_SILENTLY(mpfr::Operation::Asinpi, x,
- result, TOLERANCE + 0.5,
- rounding_mode)) {
+ if (!TEST_MPFR_MATCH_ROUNDING_SILENTLY(mpfr::Operation::Asinpi, x, result,
+ TOLERANCE + 0.5, rounding_mode)) {
++fails;
while (!TEST_MPFR_MATCH_ROUNDING_SILENTLY(mpfr::Operation::Asinpi, x,
- result, tol,
- rounding_mode)) {
+ result, tol, rounding_mode)) {
mx = x;
mr = result;
@@ -73,8 +71,7 @@ TEST_F(LlvmLibcAsinpiTest, InDoubleRange) {
if (fails) {
tlog << " Asinpi failed: " << fails << "/" << count << "/" << cc
<< " tests.\n";
- tlog << " Max ULPs is at most: " << static_cast<uint64_t>(tol)
- << ".\n";
+ tlog << " Max ULPs is at most: " << static_cast<uint64_t>(tol) << ".\n";
EXPECT_MPFR_MATCH(mpfr::Operation::Asinpi, mx, mr, 0.5, rounding_mode);
}
};
>From a1b0a1facf6c3782abb0bd97f232d86176d296a7 Mon Sep 17 00:00:00 2001
From: hulxv <hulxxv at gmail.com>
Date: Tue, 24 Mar 2026 07:09:33 +0200
Subject: [PATCH 3/6] fix mpfr
---
libc/src/__support/math/asinpi.h | 7 +------
1 file changed, 1 insertion(+), 6 deletions(-)
diff --git a/libc/src/__support/math/asinpi.h b/libc/src/__support/math/asinpi.h
index e492a5e64ee7e..1acf8aaea9dcf 100644
--- a/libc/src/__support/math/asinpi.h
+++ b/libc/src/__support/math/asinpi.h
@@ -42,13 +42,8 @@ LIBC_INLINE double asinpi(double x) {
// When |x| < 2^-26, asinpi(x) ~ x/pi.
// The relative error of x/pi is:
// |asinpi(x) - x/pi| / |asinpi(x)| < x^2/6 < 2^-54.
-#if defined(LIBC_MATH_HAS_SKIP_ACCURATE_PASS)
+#ifdef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
return x * ASINPI_COEFFS[0];
-#elif defined(LIBC_TARGET_CPU_HAS_FMA_DOUBLE)
- return fputil::multiply_add(x, ONE_OVER_PI_DD.hi, x * ONE_OVER_PI_DD.lo);
-#else
- double hi = x * ONE_OVER_PI_DD.hi;
- return hi + x * ONE_OVER_PI_DD.lo;
#endif // LIBC_MATH_HAS_SKIP_ACCURATE_PASS
}
>From 48703675e0f0b8f924cb42e1173877cdf7b3c01e Mon Sep 17 00:00:00 2001
From: hulxv <hulxxv at gmail.com>
Date: Wed, 25 Mar 2026 06:40:42 +0200
Subject: [PATCH 4/6] fix: raise underflow for subnormal numbers
---
libc/src/__support/math/asinpi.h | 55 ++++++++++++++++++++++++++++++++
1 file changed, 55 insertions(+)
diff --git a/libc/src/__support/math/asinpi.h b/libc/src/__support/math/asinpi.h
index 1acf8aaea9dcf..a088ae322ad79 100644
--- a/libc/src/__support/math/asinpi.h
+++ b/libc/src/__support/math/asinpi.h
@@ -53,6 +53,61 @@ LIBC_INLINE double asinpi(double x) {
using Float128 = fputil::DyadicFloat<128>;
using DoubleDouble = fputil::DoubleDouble;
+ // For |x| < 2^-511, x^2 would underflow to subnormal, raising a
+ // spurious underflow exception. Since asinpi(x) = x/pi with correction
+ // x^2/(6*pi) < 2^-1024 relative (negligible), compute x/pi directly
+ // in Float128.
+ if (LIBC_UNLIKELY(x_exp < 512)) {
+ Float128 x_f128(x);
+ Float128 r = fputil::quick_mul(x_f128, ONE_OVER_PI_F128);
+ double result = static_cast<double>(r);
+
+ // IEEE 754 "after rounding" tininess: the 53-bit unlimited-exponent
+ // result is strictly between +-2^-1022. DyadicFloat's conversion
+ // checks the *IEEE subnormal* result (52-bit at the boundary), not
+ // the 53-bit unlimited-exponent result, so we detect it here.
+ int exp_hi = r.exponent + 127 + FPBits::EXP_BIAS;
+ if (LIBC_UNLIKELY(exp_hi <= 0) && !r.mantissa.is_zero()) {
+ bool raise_underflow = true;
+ // When exp_hi == 0, a carry in 53-bit rounding can push the
+ // result to exactly 2^-1022 (not tiny). Check for this.
+ if (exp_hi == 0) {
+ constexpr unsigned SHIFT_53 = 128 - FPBits::SIG_LEN - 1;
+ using MantT = typename Float128::MantissaType;
+ MantT m53 = r.mantissa >> SHIFT_53;
+ constexpr MantT ALL_ONES_53 = (MantT(1) << (FPBits::SIG_LEN + 1)) - 1;
+ if (m53 == ALL_ONES_53) {
+ // All 53 bits set. carry happens if rounding rounds away
+ // from zero at this precision.
+ bool round_bit =
+ static_cast<bool>((r.mantissa >> (SHIFT_53 - 1)) & 1);
+ MantT sticky_mask = (MantT(1) << (SHIFT_53 - 1)) - 1;
+ bool sticky = (r.mantissa & sticky_mask) != 0;
+ bool lsb = static_cast<bool>(m53 & 1);
+ switch (fputil::quick_get_round()) {
+ case FE_TONEAREST:
+ // Carry if round_bit && (lsb || sticky) (round half to even).
+ raise_underflow = !(round_bit && (lsb || sticky));
+ break;
+ case FE_UPWARD:
+ raise_underflow = xbits.is_neg() || !(round_bit || sticky);
+ break;
+ case FE_DOWNWARD:
+ raise_underflow = !xbits.is_neg() || !(round_bit || sticky);
+ break;
+ case FE_TOWARDZERO:
+ default:
+ raise_underflow = true; // truncation never carries
+ break;
+ }
+ }
+ }
+ if (raise_underflow)
+ fputil::raise_except_if_required(FE_UNDERFLOW | FE_INEXACT);
+ }
+ return result;
+ }
+
unsigned idx = 0;
DoubleDouble x_sq = fputil::exact_mult(x, x);
double err = xbits.abs().get_val() * 0x1.0p-51;
>From 1ef36a5e8eec916f83ebf1d663f6bab8ce07a95c Mon Sep 17 00:00:00 2001
From: hulxv <hulxxv at gmail.com>
Date: Wed, 25 Mar 2026 07:13:23 +0200
Subject: [PATCH 5/6] docs: check asinpi double precision
---
libc/docs/headers/math/index.rst | 2 +-
1 file changed, 1 insertion(+), 1 deletion(-)
diff --git a/libc/docs/headers/math/index.rst b/libc/docs/headers/math/index.rst
index 7ea54fb4d8263..45166f050df93 100644
--- a/libc/docs/headers/math/index.rst
+++ b/libc/docs/headers/math/index.rst
@@ -273,7 +273,7 @@ Higher Math Functions
+-----------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+------------------------+----------------------------+
| asinh | |check| | | | |check| | | | 7.12.5.2 | F.10.2.2 |
+-----------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+------------------------+----------------------------+
-| asinpi | |check| | | | |check| | | | 7.12.4.9 | F.10.1.9 |
+| asinpi | |check| | |check| | | |check| | | | 7.12.4.9 | F.10.1.9 |
+-----------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+------------------------+----------------------------+
| atan | |check| | 1 ULP | | |check| | | | 7.12.4.3 | F.10.1.3 |
+-----------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+------------------------+----------------------------+
>From 209eff4c580b0658d4357f8051c50d87d018ac04 Mon Sep 17 00:00:00 2001
From: hulxv <hulxxv at gmail.com>
Date: Fri, 27 Mar 2026 07:11:27 +0200
Subject: [PATCH 6/6] add the estimated relative error of asinpi
---
libc/src/__support/math/asin_utils.h | 2 ++
1 file changed, 2 insertions(+)
diff --git a/libc/src/__support/math/asin_utils.h b/libc/src/__support/math/asin_utils.h
index e39262f7d19a8..faeb38cefbdd4 100644
--- a/libc/src/__support/math/asin_utils.h
+++ b/libc/src/__support/math/asin_utils.h
@@ -69,6 +69,8 @@ LIBC_INLINE double asin_eval(double u) {
// > g = asin(x) / (pi * x);
// > P = fpminimax(g, [|0, 2, 4, 6, 8, 10, 12, 14, 16, 18, 20, 22|],
// > [|D...|], [0, 0.5]);
+//
+// Error : 0x1.a53f84eafa3ea69bb81b6c52b3278872083fca2c757bd778acp-54 ~= 2^-54
LIBC_INLINE_VAR constexpr double ASINPI_COEFFS[12] = {
0x1.45f306dc9c881p-2, // x^0
0x1.b2995e7b7e756p-5, // x^2
More information about the libc-commits
mailing list