[libc][math][c23] implement double-precision asinpi (#188158)
Implement the double precision version of the asinpi c23 math function
This commit is contained in:
parent
2e6e36b173
commit
1bb03026b6
@ -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
|
||||
|
||||
@ -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
|
||||
|
||||
@ -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
|
||||
|
||||
@ -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
|
||||
|
||||
@ -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
|
||||
|
||||
@ -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
|
||||
|
||||
@ -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
|
||||
|
||||
@ -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
|
||||
|
||||
@ -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
|
||||
|
||||
@ -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
|
||||
|
||||
@ -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
|
||||
|
||||
@ -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 |
|
||||
+-----------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+------------------------+----------------------------+
|
||||
|
||||
@ -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
|
||||
|
||||
@ -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"
|
||||
|
||||
23
libc/shared/math/asinpi.h
Normal file
23
libc/shared/math/asinpi.h
Normal file
@ -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
|
||||
@ -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
|
||||
|
||||
@ -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
|
||||
|
||||
291
libc/src/__support/math/asinpi.h
Normal file
291
libc/src/__support/math/asinpi.h
Normal file
@ -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<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.
|
||||
#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<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;
|
||||
// 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
|
||||
@ -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)
|
||||
|
||||
|
||||
20
libc/src/math/asinpi.h
Normal file
20
libc/src/math/asinpi.h
Normal file
@ -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
|
||||
@ -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
|
||||
|
||||
16
libc/src/math/generic/asinpi.cpp
Normal file
16
libc/src/math/generic/asinpi.cpp
Normal file
@ -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
|
||||
@ -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
|
||||
|
||||
@ -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));
|
||||
|
||||
@ -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
|
||||
|
||||
92
libc/test/src/math/asinpi_test.cpp
Normal file
92
libc/test/src/math/asinpi_test.cpp
Normal file
@ -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<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
|
||||
}
|
||||
@ -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
|
||||
|
||||
52
libc/test/src/math/smoke/asinpi_test.cpp
Normal file
52
libc/test/src/math/smoke/asinpi_test.cpp
Normal file
@ -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
|
||||
@ -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 = [
|
||||
|
||||
Loading…
x
Reference in New Issue
Block a user