1 //===-- Single-precision acos 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/acosf.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 19 #include "inv_trigf_utils.h" 20 21 namespace LIBC_NAMESPACE_DECL { 22 23 static constexpr size_t N_EXCEPTS = 4; 24 25 // Exceptional values when |x| <= 0.5 26 static constexpr fputil::ExceptValues<float, N_EXCEPTS> ACOSF_EXCEPTS = {{ 27 // (inputs, RZ output, RU offset, RD offset, RN offset) 28 // x = 0x1.110b46p-26, acosf(x) = 0x1.921fb4p0 (RZ) 29 {0x328885a3, 0x3fc90fda, 1, 0, 1}, 30 // x = -0x1.110b46p-26, acosf(x) = 0x1.921fb4p0 (RZ) 31 {0xb28885a3, 0x3fc90fda, 1, 0, 1}, 32 // x = 0x1.04c444p-12, acosf(x) = 0x1.920f68p0 (RZ) 33 {0x39826222, 0x3fc907b4, 1, 0, 1}, 34 // x = -0x1.04c444p-12, acosf(x) = 0x1.923p0 (RZ) 35 {0xb9826222, 0x3fc91800, 1, 0, 1}, 36 }}; 37 38 LLVM_LIBC_FUNCTION(float, acosf, (float x)) { 39 using FPBits = typename fputil::FPBits<float>; 40 41 FPBits xbits(x); 42 uint32_t x_uint = xbits.uintval(); 43 uint32_t x_abs = xbits.uintval() & 0x7fff'ffffU; 44 uint32_t x_sign = x_uint >> 31; 45 46 // |x| <= 0.5 47 if (LIBC_UNLIKELY(x_abs <= 0x3f00'0000U)) { 48 // |x| < 0x1p-10 49 if (LIBC_UNLIKELY(x_abs < 0x3a80'0000U)) { 50 // When |x| < 2^-10, we use the following approximation: 51 // acos(x) = pi/2 - asin(x) 52 // ~ pi/2 - x - x^3 / 6 53 54 // Check for exceptional values 55 if (auto r = ACOSF_EXCEPTS.lookup(x_uint); LIBC_UNLIKELY(r.has_value())) 56 return r.value(); 57 58 double xd = static_cast<double>(x); 59 return static_cast<float>(fputil::multiply_add( 60 -0x1.5555555555555p-3 * xd, xd * xd, M_MATH_PI_2 - xd)); 61 } 62 63 // For |x| <= 0.5, we approximate acosf(x) by: 64 // acos(x) = pi/2 - asin(x) = pi/2 - x * P(x^2) 65 // Where P(X^2) = Q(X) is a degree-20 minimax even polynomial approximating 66 // asin(x)/x on [0, 0.5] generated by Sollya with: 67 // > Q = fpminimax(asin(x)/x, [|0, 2, 4, 6, 8, 10, 12, 14, 16, 18, 20|], 68 // [|1, D...|], [0, 0.5]); 69 double xd = static_cast<double>(x); 70 double xsq = xd * xd; 71 double x3 = xd * xsq; 72 double r = asin_eval(xsq); 73 return static_cast<float>(fputil::multiply_add(-x3, r, M_MATH_PI_2 - xd)); 74 } 75 76 // |x| >= 1, return 0, 2pi, or NaNs. 77 if (LIBC_UNLIKELY(x_abs >= 0x3f80'0000U)) { 78 if (x_abs == 0x3f80'0000U) 79 return x_sign ? /* x == -1.0f */ fputil::round_result_slightly_down( 80 0x1.921fb6p+1f) 81 : /* x == 1.0f */ 0.0f; 82 83 if (x_abs <= 0x7f80'0000U) { 84 fputil::set_errno_if_required(EDOM); 85 fputil::raise_except_if_required(FE_INVALID); 86 } 87 return x + FPBits::quiet_nan().get_val(); 88 } 89 90 // When 0.5 < |x| < 1, we perform range reduction as follow: 91 // 92 // Assume further that 0.5 < x <= 1, and let: 93 // y = acos(x) 94 // We use the double angle formula: 95 // x = cos(y) = 1 - 2 sin^2(y/2) 96 // So: 97 // sin(y/2) = sqrt( (1 - x)/2 ) 98 // And hence: 99 // y = 2 * asin( sqrt( (1 - x)/2 ) ) 100 // Let u = (1 - x)/2, then 101 // acos(x) = 2 * asin( sqrt(u) ) 102 // Moreover, since 0.5 < x <= 1, 103 // 0 <= u < 1/4, and 0 <= sqrt(u) < 0.5, 104 // And hence we can reuse the same polynomial approximation of asin(x) when 105 // |x| <= 0.5: 106 // acos(x) ~ 2 * sqrt(u) * P(u). 107 // 108 // When -1 < x <= -0.5, we use the identity: 109 // acos(x) = pi - acos(-x) 110 // which is reduced to the postive case. 111 112 xbits.set_sign(Sign::POS); 113 double xd = static_cast<double>(xbits.get_val()); 114 double u = fputil::multiply_add(-0.5, xd, 0.5); 115 double cv = 2 * fputil::sqrt<double>(u); 116 117 double r3 = asin_eval(u); 118 double r = fputil::multiply_add(cv * u, r3, cv); 119 return static_cast<float>(x_sign ? M_MATH_PI - r : r); 120 } 121 122 } // namespace LIBC_NAMESPACE_DECL 123