diff --git a/libc/config/baremetal/aarch64/entrypoints.txt b/libc/config/baremetal/aarch64/entrypoints.txt index 19f3fb189de7..ed84bdff8ce1 100644 --- a/libc/config/baremetal/aarch64/entrypoints.txt +++ b/libc/config/baremetal/aarch64/entrypoints.txt @@ -349,6 +349,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 902f6ed44179..43ba426b8ca0 100644 --- a/libc/config/baremetal/arm/entrypoints.txt +++ b/libc/config/baremetal/arm/entrypoints.txt @@ -359,6 +359,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 5c2b5c24e158..79fbd1e4b0de 100644 --- a/libc/config/baremetal/riscv/entrypoints.txt +++ b/libc/config/baremetal/riscv/entrypoints.txt @@ -355,6 +355,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 8e0d2fb01966..71881d23c7a0 100644 --- a/libc/config/darwin/aarch64/entrypoints.txt +++ b/libc/config/darwin/aarch64/entrypoints.txt @@ -168,6 +168,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 028d57b52116..e327d9ac6f3d 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 feae60547bed..efc784d23e30 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 49288ce0699a..4311ba624e6d 100644 --- a/libc/config/linux/aarch64/entrypoints.txt +++ b/libc/config/linux/aarch64/entrypoints.txt @@ -443,6 +443,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 9fd45329ad0f..a10dc0c82e36 100644 --- a/libc/config/linux/arm/entrypoints.txt +++ b/libc/config/linux/arm/entrypoints.txt @@ -268,6 +268,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 2f5a742a272d..17efa2e368ed 100644 --- a/libc/config/linux/riscv/entrypoints.txt +++ b/libc/config/linux/riscv/entrypoints.txt @@ -447,6 +447,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 2d7de065d811..eeaccc863dbc 100644 --- a/libc/config/linux/x86_64/entrypoints.txt +++ b/libc/config/linux/x86_64/entrypoints.txt @@ -497,6 +497,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 ee8072a3465e..e594f6bbb1fd 100644 --- a/libc/config/windows/entrypoints.txt +++ b/libc/config/windows/entrypoints.txt @@ -148,6 +148,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/docs/headers/math/index.rst b/libc/docs/headers/math/index.rst index 88b284d1e63f..f9614c6d19d2 100644 --- a/libc/docs/headers/math/index.rst +++ b/libc/docs/headers/math/index.rst @@ -269,7 +269,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| | | |check| | 7.12.4.3 | F.10.1.3 | +-----------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+------------------------+----------------------------+ diff --git a/libc/include/math.yaml b/libc/include/math.yaml index dbb8479b1d0d..3984a665b234 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 71da66d37bab..08a4c994283c 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 000000000000..42ada49e0a1a --- /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 187ff4f124b2..10329830a004 100644 --- a/libc/src/__support/math/CMakeLists.txt +++ b/libc/src/__support/math/CMakeLists.txt @@ -189,6 +189,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 2feb1485ec66..4f64ad529171 100644 --- a/libc/src/__support/math/asin_utils.h +++ b/libc/src/__support/math/asin_utils.h @@ -63,6 +63,47 @@ 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]); +// > print("Error:", dirtyinfnorm(P - g, [1e-30; 0.25])); +// Error : 0x1.a53f84eafa3ea69bb81b6c52b3278872083fca2c757bd778acp-54 ~= 2^-54 +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 +606,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 000000000000..a088ae322ad7 --- /dev/null +++ b/libc/src/__support/math/asinpi.h @@ -0,0 +1,291 @@ +//===-- 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; + + 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. +#ifdef LIBC_MATH_HAS_SKIP_ACCURATE_PASS + return x * ASINPI_COEFFS[0]; +#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; + + // 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(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((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(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; + // 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(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(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(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(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(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(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(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(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 c7ded04aad7c..ebba3f04b8d5 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 000000000000..7b8770b5e52a --- /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 ad05fd83cadb..30f33749d165 100644 --- a/libc/src/math/generic/CMakeLists.txt +++ b/libc/src/math/generic/CMakeLists.txt @@ -4037,6 +4037,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 000000000000..20b76ff1ffd6 --- /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 b807ea9ac602..c3795d366a58 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 31b14bbe204c..b52abc6d8e99 100644 --- a/libc/test/shared/shared_math_test.cpp +++ b/libc/test/shared/shared_math_test.cpp @@ -242,6 +242,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 e157b0463bbd..2fd77f964e4d 100644 --- a/libc/test/src/math/CMakeLists.txt +++ b/libc/test/src/math/CMakeLists.txt @@ -2523,6 +2523,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 000000000000..75dfbc6608ea --- /dev/null +++ b/libc/test/src/math/asinpi_test.cpp @@ -0,0 +1,92 @@ +//===-- 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; + +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(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 544425a39d34..79f2ea480637 100644 --- a/libc/test/src/math/smoke/CMakeLists.txt +++ b/libc/test/src/math/smoke/CMakeLists.txt @@ -4757,6 +4757,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 000000000000..fcd3dd3e5ddc --- /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; + +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 29792cec1149..a69ebcdb63fa 100644 --- a/utils/bazel/llvm-project-overlay/libc/BUILD.bazel +++ b/utils/bazel/llvm-project-overlay/libc/BUILD.bazel @@ -3189,6 +3189,25 @@ libc_support_library( ], ) +libc_support_library( + name = "__support_math_asinpi", + hdrs = ["src/__support/math/asinpi.h"], + deps = [ + ":__support_fputil_cast", + ":__support_fputil_double_double", + ":__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"], @@ -6566,6 +6585,13 @@ libc_math_function( ], ) +libc_math_function( + name = "asinpi", + additional_deps = [ + ":__support_math_asinpi", + ], +) + libc_math_function( name = "asinpif", additional_deps = [