176 lines · c
1//===-- Implementation header for asinf -------------------------*- C++ -*-===//2//3// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.4// See https://llvm.org/LICENSE.txt for license information.5// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception6//7//===----------------------------------------------------------------------===//8 9#ifndef LLVM_LIBC_SRC___SUPPORT_MATH_ASINF_H10#define LLVM_LIBC_SRC___SUPPORT_MATH_ASINF_H11 12#include "inv_trigf_utils.h"13#include "src/__support/FPUtil/FEnvImpl.h"14#include "src/__support/FPUtil/FPBits.h"15#include "src/__support/FPUtil/except_value_utils.h"16#include "src/__support/FPUtil/multiply_add.h"17#include "src/__support/FPUtil/sqrt.h"18#include "src/__support/macros/config.h"19#include "src/__support/macros/optimization.h" // LIBC_UNLIKELY20#include "src/__support/macros/properties/cpu_features.h" // LIBC_TARGET_CPU_HAS_FMA21 22namespace LIBC_NAMESPACE_DECL {23 24namespace math {25 26LIBC_INLINE static constexpr float asinf(float x) {27 using namespace inv_trigf_utils_internal;28 29#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS30 constexpr size_t N_EXCEPTS = 2;31 32 // Exceptional values when |x| <= 0.533 constexpr fputil::ExceptValues<float, N_EXCEPTS> ASINF_EXCEPTS_LO = {{34 // (inputs, RZ output, RU offset, RD offset, RN offset)35 // x = 0x1.137f0cp-5, asinf(x) = 0x1.138c58p-5 (RZ)36 {0x3d09bf86, 0x3d09c62c, 1, 0, 1},37 // x = 0x1.cbf43cp-4, asinf(x) = 0x1.cced1cp-4 (RZ)38 {0x3de5fa1e, 0x3de6768e, 1, 0, 0},39 }};40 41 // Exceptional values when 0.5 < |x| <= 142 constexpr fputil::ExceptValues<float, N_EXCEPTS> ASINF_EXCEPTS_HI = {{43 // (inputs, RZ output, RU offset, RD offset, RN offset)44 // x = 0x1.107434p-1, asinf(x) = 0x1.1f4b64p-1 (RZ)45 {0x3f083a1a, 0x3f0fa5b2, 1, 0, 0},46 // x = 0x1.ee836cp-1, asinf(x) = 0x1.4f0654p0 (RZ)47 {0x3f7741b6, 0x3fa7832a, 1, 0, 0},48 }};49#endif // !LIBC_MATH_HAS_SKIP_ACCURATE_PASS50 51 using namespace inv_trigf_utils_internal;52 using FPBits = typename fputil::FPBits<float>;53 54 FPBits xbits(x);55 uint32_t x_uint = xbits.uintval();56 uint32_t x_abs = xbits.uintval() & 0x7fff'ffffU;57 constexpr double SIGN[2] = {1.0, -1.0};58 uint32_t x_sign = x_uint >> 31;59 60 // |x| <= 0.5-ish61 if (x_abs < 0x3f04'471dU) {62 // |x| < 0x1.d12edp-1263 if (LIBC_UNLIKELY(x_abs < 0x39e8'9768U)) {64 // When |x| < 2^-12, the relative error of the approximation asin(x) ~ x65 // is:66 // |asin(x) - x| / |asin(x)| < |x^3| / (6|x|)67 // = x^2 / 668 // < 2^-2569 // < epsilon(1)/2.70 // So the correctly rounded values of asin(x) are:71 // = x + sign(x)*eps(x) if rounding mode = FE_TOWARDZERO,72 // or (rounding mode = FE_UPWARD and x is73 // negative),74 // = x otherwise.75 // To simplify the rounding decision and make it more efficient, we use76 // fma(x, 2^-25, x) instead.77 // An exhaustive test shows that this formula work correctly for all78 // rounding modes up to |x| < 0x1.d12edp-12.79 // Note: to use the formula x + 2^-25*x to decide the correct rounding, we80 // do need fma(x, 2^-25, x) to prevent underflow caused by 2^-25*x when81 // |x| < 2^-125. For targets without FMA instructions, we simply use82 // double for intermediate results as it is more efficient than using an83 // emulated version of FMA.84#if defined(LIBC_TARGET_CPU_HAS_FMA_FLOAT)85 return fputil::multiply_add(x, 0x1.0p-25f, x);86#else87 double xd = static_cast<double>(x);88 return static_cast<float>(fputil::multiply_add(xd, 0x1.0p-25, xd));89#endif // LIBC_TARGET_CPU_HAS_FMA_FLOAT90 }91 92#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS93 // Check for exceptional values94 if (auto r = ASINF_EXCEPTS_LO.lookup_odd(x_abs, x_sign);95 LIBC_UNLIKELY(r.has_value()))96 return r.value();97#endif // !LIBC_MATH_HAS_SKIP_ACCURATE_PASS98 99 // For |x| <= 0.5, we approximate asinf(x) by:100 // asin(x) = x * P(x^2)101 // Where P(X^2) = Q(X) is a degree-20 minimax even polynomial approximating102 // asin(x)/x on [0, 0.5] generated by Sollya with:103 // > Q = fpminimax(asin(x)/x, [|0, 2, 4, 6, 8, 10, 12, 14, 16, 18, 20|],104 // [|1, D...|], [0, 0.5]);105 // An exhaustive test shows that this approximation works well up to a106 // little more than 0.5.107 double xd = static_cast<double>(x);108 double xsq = xd * xd;109 double x3 = xd * xsq;110 double r = asin_eval(xsq);111 return static_cast<float>(fputil::multiply_add(x3, r, xd));112 }113 114 // |x| > 1, return NaNs.115 if (LIBC_UNLIKELY(x_abs > 0x3f80'0000U)) {116 if (xbits.is_signaling_nan()) {117 fputil::raise_except_if_required(FE_INVALID);118 return FPBits::quiet_nan().get_val();119 }120 121 if (x_abs <= 0x7f80'0000U) {122 fputil::set_errno_if_required(EDOM);123 fputil::raise_except_if_required(FE_INVALID);124 }125 126 return FPBits::quiet_nan().get_val();127 }128 129#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS130 // Check for exceptional values131 if (auto r = ASINF_EXCEPTS_HI.lookup_odd(x_abs, x_sign);132 LIBC_UNLIKELY(r.has_value()))133 return r.value();134#endif // !LIBC_MATH_HAS_SKIP_ACCURATE_PASS135 136 // When |x| > 0.5, we perform range reduction as follow:137 //138 // Assume further that 0.5 < x <= 1, and let:139 // y = asin(x)140 // We will use the double angle formula:141 // cos(2y) = 1 - 2 sin^2(y)142 // and the complement angle identity:143 // x = sin(y) = cos(pi/2 - y)144 // = 1 - 2 sin^2 (pi/4 - y/2)145 // So:146 // sin(pi/4 - y/2) = sqrt( (1 - x)/2 )147 // And hence:148 // pi/4 - y/2 = asin( sqrt( (1 - x)/2 ) )149 // Equivalently:150 // asin(x) = y = pi/2 - 2 * asin( sqrt( (1 - x)/2 ) )151 // Let u = (1 - x)/2, then:152 // asin(x) = pi/2 - 2 * asin( sqrt(u) )153 // Moreover, since 0.5 < x <= 1:154 // 0 <= u < 1/4, and 0 <= sqrt(u) < 0.5,155 // And hence we can reuse the same polynomial approximation of asin(x) when156 // |x| <= 0.5:157 // asin(x) ~ pi/2 - 2 * sqrt(u) * P(u),158 159 xbits.set_sign(Sign::POS);160 double sign = SIGN[x_sign];161 double xd = static_cast<double>(xbits.get_val());162 double u = fputil::multiply_add(-0.5, xd, 0.5);163 double c1 = sign * (-2 * fputil::sqrt<double>(u));164 double c2 = fputil::multiply_add(sign, M_MATH_PI_2, c1);165 double c3 = c1 * u;166 167 double r = asin_eval(u);168 return static_cast<float>(fputil::multiply_add(c3, r, c2));169}170 171} // namespace math172 173} // namespace LIBC_NAMESPACE_DECL174 175#endif // LLVM_LIBC_SRC___SUPPORT_MATH_ASINF_H176