1 //===-- Single-precision asin function ------------------------------------===// 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-exception 6 // 7 //===----------------------------------------------------------------------===// 8 9 #include "src/math/asinf.h" 10 #include "src/__support/FPUtil/FEnvImpl.h" 11 #include "src/__support/FPUtil/FPBits.h" 12 #include "src/__support/FPUtil/PolyEval.h" 13 #include "src/__support/FPUtil/except_value_utils.h" 14 #include "src/__support/FPUtil/multiply_add.h" 15 #include "src/__support/FPUtil/sqrt.h" 16 #include "src/__support/macros/config.h" 17 #include "src/__support/macros/optimization.h" // LIBC_UNLIKELY 18 #include "src/__support/macros/properties/cpu_features.h" // LIBC_TARGET_CPU_HAS_FMA 19 20 #include "inv_trigf_utils.h" 21 22 namespace LIBC_NAMESPACE_DECL { 23 24 static constexpr size_t N_EXCEPTS = 2; 25 26 // Exceptional values when |x| <= 0.5 27 static constexpr fputil::ExceptValues<float, N_EXCEPTS> ASINF_EXCEPTS_LO = {{ 28 // (inputs, RZ output, RU offset, RD offset, RN offset) 29 // x = 0x1.137f0cp-5, asinf(x) = 0x1.138c58p-5 (RZ) 30 {0x3d09bf86, 0x3d09c62c, 1, 0, 1}, 31 // x = 0x1.cbf43cp-4, asinf(x) = 0x1.cced1cp-4 (RZ) 32 {0x3de5fa1e, 0x3de6768e, 1, 0, 0}, 33 }}; 34 35 // Exceptional values when 0.5 < |x| <= 1 36 static constexpr fputil::ExceptValues<float, N_EXCEPTS> ASINF_EXCEPTS_HI = {{ 37 // (inputs, RZ output, RU offset, RD offset, RN offset) 38 // x = 0x1.107434p-1, asinf(x) = 0x1.1f4b64p-1 (RZ) 39 {0x3f083a1a, 0x3f0fa5b2, 1, 0, 0}, 40 // x = 0x1.ee836cp-1, asinf(x) = 0x1.4f0654p0 (RZ) 41 {0x3f7741b6, 0x3fa7832a, 1, 0, 0}, 42 }}; 43 44 LLVM_LIBC_FUNCTION(float, asinf, (float x)) { 45 using FPBits = typename fputil::FPBits<float>; 46 47 FPBits xbits(x); 48 uint32_t x_uint = xbits.uintval(); 49 uint32_t x_abs = xbits.uintval() & 0x7fff'ffffU; 50 constexpr double SIGN[2] = {1.0, -1.0}; 51 uint32_t x_sign = x_uint >> 31; 52 53 // |x| <= 0.5-ish 54 if (x_abs < 0x3f04'471dU) { 55 // |x| < 0x1.d12edp-12 56 if (LIBC_UNLIKELY(x_abs < 0x39e8'9768U)) { 57 // When |x| < 2^-12, the relative error of the approximation asin(x) ~ x 58 // is: 59 // |asin(x) - x| / |asin(x)| < |x^3| / (6|x|) 60 // = x^2 / 6 61 // < 2^-25 62 // < epsilon(1)/2. 63 // So the correctly rounded values of asin(x) are: 64 // = x + sign(x)*eps(x) if rounding mode = FE_TOWARDZERO, 65 // or (rounding mode = FE_UPWARD and x is 66 // negative), 67 // = x otherwise. 68 // To simplify the rounding decision and make it more efficient, we use 69 // fma(x, 2^-25, x) instead. 70 // An exhaustive test shows that this formula work correctly for all 71 // rounding modes up to |x| < 0x1.d12edp-12. 72 // Note: to use the formula x + 2^-25*x to decide the correct rounding, we 73 // do need fma(x, 2^-25, x) to prevent underflow caused by 2^-25*x when 74 // |x| < 2^-125. For targets without FMA instructions, we simply use 75 // double for intermediate results as it is more efficient than using an 76 // emulated version of FMA. 77 #if defined(LIBC_TARGET_CPU_HAS_FMA) 78 return fputil::multiply_add(x, 0x1.0p-25f, x); 79 #else 80 double xd = static_cast<double>(x); 81 return static_cast<float>(fputil::multiply_add(xd, 0x1.0p-25, xd)); 82 #endif // LIBC_TARGET_CPU_HAS_FMA 83 } 84 85 // Check for exceptional values 86 if (auto r = ASINF_EXCEPTS_LO.lookup_odd(x_abs, x_sign); 87 LIBC_UNLIKELY(r.has_value())) 88 return r.value(); 89 90 // For |x| <= 0.5, we approximate asinf(x) by: 91 // asin(x) = x * P(x^2) 92 // Where P(X^2) = Q(X) is a degree-20 minimax even polynomial approximating 93 // asin(x)/x on [0, 0.5] generated by Sollya with: 94 // > Q = fpminimax(asin(x)/x, [|0, 2, 4, 6, 8, 10, 12, 14, 16, 18, 20|], 95 // [|1, D...|], [0, 0.5]); 96 // An exhaustive test shows that this approximation works well up to a 97 // little more than 0.5. 98 double xd = static_cast<double>(x); 99 double xsq = xd * xd; 100 double x3 = xd * xsq; 101 double r = asin_eval(xsq); 102 return static_cast<float>(fputil::multiply_add(x3, r, xd)); 103 } 104 105 // |x| > 1, return NaNs. 106 if (LIBC_UNLIKELY(x_abs > 0x3f80'0000U)) { 107 if (x_abs <= 0x7f80'0000U) { 108 fputil::set_errno_if_required(EDOM); 109 fputil::raise_except_if_required(FE_INVALID); 110 } 111 return FPBits::quiet_nan().get_val(); 112 } 113 114 // Check for exceptional values 115 if (auto r = ASINF_EXCEPTS_HI.lookup_odd(x_abs, x_sign); 116 LIBC_UNLIKELY(r.has_value())) 117 return r.value(); 118 119 // When |x| > 0.5, we perform range reduction as follow: 120 // 121 // Assume further that 0.5 < x <= 1, and let: 122 // y = asin(x) 123 // We will use the double angle formula: 124 // cos(2y) = 1 - 2 sin^2(y) 125 // and the complement angle identity: 126 // x = sin(y) = cos(pi/2 - y) 127 // = 1 - 2 sin^2 (pi/4 - y/2) 128 // So: 129 // sin(pi/4 - y/2) = sqrt( (1 - x)/2 ) 130 // And hence: 131 // pi/4 - y/2 = asin( sqrt( (1 - x)/2 ) ) 132 // Equivalently: 133 // asin(x) = y = pi/2 - 2 * asin( sqrt( (1 - x)/2 ) ) 134 // Let u = (1 - x)/2, then: 135 // asin(x) = pi/2 - 2 * asin( sqrt(u) ) 136 // Moreover, since 0.5 < x <= 1: 137 // 0 <= u < 1/4, and 0 <= sqrt(u) < 0.5, 138 // And hence we can reuse the same polynomial approximation of asin(x) when 139 // |x| <= 0.5: 140 // asin(x) ~ pi/2 - 2 * sqrt(u) * P(u), 141 142 xbits.set_sign(Sign::POS); 143 double sign = SIGN[x_sign]; 144 double xd = static_cast<double>(xbits.get_val()); 145 double u = fputil::multiply_add(-0.5, xd, 0.5); 146 double c1 = sign * (-2 * fputil::sqrt<double>(u)); 147 double c2 = fputil::multiply_add(sign, M_MATH_PI_2, c1); 148 double c3 = c1 * u; 149 150 double r = asin_eval(u); 151 return static_cast<float>(fputil::multiply_add(c3, r, c2)); 152 } 153 154 } // namespace LIBC_NAMESPACE_DECL 155