aboutsummaryrefslogtreecommitdiff
path: root/libc/src/math/generic/asinf.cpp
diff options
context:
space:
mode:
Diffstat (limited to 'libc/src/math/generic/asinf.cpp')
-rw-r--r--libc/src/math/generic/asinf.cpp154
1 files changed, 2 insertions, 152 deletions
diff --git a/libc/src/math/generic/asinf.cpp b/libc/src/math/generic/asinf.cpp
index 12383bf..9c6766f 100644
--- a/libc/src/math/generic/asinf.cpp
+++ b/libc/src/math/generic/asinf.cpp
@@ -7,160 +7,10 @@
//===----------------------------------------------------------------------===//
#include "src/math/asinf.h"
-#include "src/__support/FPUtil/FEnvImpl.h"
-#include "src/__support/FPUtil/FPBits.h"
-#include "src/__support/FPUtil/PolyEval.h"
-#include "src/__support/FPUtil/except_value_utils.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 "inv_trigf_utils.h"
+#include "src/__support/math/asinf.h"
namespace LIBC_NAMESPACE_DECL {
-#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
-static constexpr size_t N_EXCEPTS = 2;
-
-// Exceptional values when |x| <= 0.5
-static constexpr fputil::ExceptValues<float, N_EXCEPTS> ASINF_EXCEPTS_LO = {{
- // (inputs, RZ output, RU offset, RD offset, RN offset)
- // x = 0x1.137f0cp-5, asinf(x) = 0x1.138c58p-5 (RZ)
- {0x3d09bf86, 0x3d09c62c, 1, 0, 1},
- // x = 0x1.cbf43cp-4, asinf(x) = 0x1.cced1cp-4 (RZ)
- {0x3de5fa1e, 0x3de6768e, 1, 0, 0},
-}};
-
-// Exceptional values when 0.5 < |x| <= 1
-static constexpr fputil::ExceptValues<float, N_EXCEPTS> ASINF_EXCEPTS_HI = {{
- // (inputs, RZ output, RU offset, RD offset, RN offset)
- // x = 0x1.107434p-1, asinf(x) = 0x1.1f4b64p-1 (RZ)
- {0x3f083a1a, 0x3f0fa5b2, 1, 0, 0},
- // x = 0x1.ee836cp-1, asinf(x) = 0x1.4f0654p0 (RZ)
- {0x3f7741b6, 0x3fa7832a, 1, 0, 0},
-}};
-#endif // !LIBC_MATH_HAS_SKIP_ACCURATE_PASS
-
-LLVM_LIBC_FUNCTION(float, asinf, (float x)) {
- using FPBits = typename fputil::FPBits<float>;
-
- FPBits xbits(x);
- uint32_t x_uint = xbits.uintval();
- uint32_t x_abs = xbits.uintval() & 0x7fff'ffffU;
- constexpr double SIGN[2] = {1.0, -1.0};
- uint32_t x_sign = x_uint >> 31;
-
- // |x| <= 0.5-ish
- if (x_abs < 0x3f04'471dU) {
- // |x| < 0x1.d12edp-12
- if (LIBC_UNLIKELY(x_abs < 0x39e8'9768U)) {
- // When |x| < 2^-12, the relative error of the approximation asin(x) ~ x
- // is:
- // |asin(x) - x| / |asin(x)| < |x^3| / (6|x|)
- // = x^2 / 6
- // < 2^-25
- // < epsilon(1)/2.
- // So the correctly rounded values of asin(x) are:
- // = x + sign(x)*eps(x) if rounding mode = FE_TOWARDZERO,
- // or (rounding mode = FE_UPWARD and x is
- // negative),
- // = x otherwise.
- // To simplify the rounding decision and make it more efficient, we use
- // fma(x, 2^-25, x) instead.
- // An exhaustive test shows that this formula work correctly for all
- // rounding modes up to |x| < 0x1.d12edp-12.
- // Note: to use the formula x + 2^-25*x to decide the correct rounding, we
- // do need fma(x, 2^-25, x) to prevent underflow caused by 2^-25*x when
- // |x| < 2^-125. For targets without FMA instructions, we simply use
- // double for intermediate results as it is more efficient than using an
- // emulated version of FMA.
-#if defined(LIBC_TARGET_CPU_HAS_FMA_FLOAT)
- return fputil::multiply_add(x, 0x1.0p-25f, x);
-#else
- double xd = static_cast<double>(x);
- return static_cast<float>(fputil::multiply_add(xd, 0x1.0p-25, xd));
-#endif // LIBC_TARGET_CPU_HAS_FMA_FLOAT
- }
-
-#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
- // Check for exceptional values
- if (auto r = ASINF_EXCEPTS_LO.lookup_odd(x_abs, x_sign);
- LIBC_UNLIKELY(r.has_value()))
- return r.value();
-#endif // !LIBC_MATH_HAS_SKIP_ACCURATE_PASS
-
- // For |x| <= 0.5, we approximate asinf(x) by:
- // asin(x) = x * P(x^2)
- // Where P(X^2) = Q(X) is a degree-20 minimax even polynomial approximating
- // asin(x)/x on [0, 0.5] generated by Sollya with:
- // > Q = fpminimax(asin(x)/x, [|0, 2, 4, 6, 8, 10, 12, 14, 16, 18, 20|],
- // [|1, D...|], [0, 0.5]);
- // An exhaustive test shows that this approximation works well up to a
- // little more than 0.5.
- double xd = static_cast<double>(x);
- double xsq = xd * xd;
- double x3 = xd * xsq;
- double r = asin_eval(xsq);
- return static_cast<float>(fputil::multiply_add(x3, r, xd));
- }
-
- // |x| > 1, return NaNs.
- if (LIBC_UNLIKELY(x_abs > 0x3f80'0000U)) {
- if (xbits.is_signaling_nan()) {
- fputil::raise_except_if_required(FE_INVALID);
- return FPBits::quiet_nan().get_val();
- }
-
- if (x_abs <= 0x7f80'0000U) {
- fputil::set_errno_if_required(EDOM);
- fputil::raise_except_if_required(FE_INVALID);
- }
-
- return FPBits::quiet_nan().get_val();
- }
-
-#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
- // Check for exceptional values
- if (auto r = ASINF_EXCEPTS_HI.lookup_odd(x_abs, x_sign);
- LIBC_UNLIKELY(r.has_value()))
- return r.value();
-#endif // !LIBC_MATH_HAS_SKIP_ACCURATE_PASS
-
- // When |x| > 0.5, we perform range reduction as follow:
- //
- // Assume further that 0.5 < x <= 1, and let:
- // y = asin(x)
- // We will use the double angle formula:
- // cos(2y) = 1 - 2 sin^2(y)
- // and the complement angle identity:
- // x = sin(y) = cos(pi/2 - y)
- // = 1 - 2 sin^2 (pi/4 - y/2)
- // So:
- // sin(pi/4 - y/2) = sqrt( (1 - x)/2 )
- // And hence:
- // pi/4 - y/2 = asin( sqrt( (1 - x)/2 ) )
- // Equivalently:
- // asin(x) = y = pi/2 - 2 * asin( sqrt( (1 - x)/2 ) )
- // Let u = (1 - x)/2, then:
- // asin(x) = pi/2 - 2 * asin( sqrt(u) )
- // Moreover, since 0.5 < x <= 1:
- // 0 <= u < 1/4, and 0 <= sqrt(u) < 0.5,
- // And hence we can reuse the same polynomial approximation of asin(x) when
- // |x| <= 0.5:
- // asin(x) ~ pi/2 - 2 * sqrt(u) * P(u),
-
- xbits.set_sign(Sign::POS);
- double sign = SIGN[x_sign];
- double xd = static_cast<double>(xbits.get_val());
- double u = fputil::multiply_add(-0.5, xd, 0.5);
- double c1 = sign * (-2 * fputil::sqrt<double>(u));
- double c2 = fputil::multiply_add(sign, M_MATH_PI_2, c1);
- double c3 = c1 * u;
-
- double r = asin_eval(u);
- return static_cast<float>(fputil::multiply_add(c3, r, c2));
-}
+LLVM_LIBC_FUNCTION(float, asinf, (float x)) { return math::asinf(x); }
} // namespace LIBC_NAMESPACE_DECL
generated by cgit v1.1 at 2025-07-27 06:59:35 +0000