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

source code of libc/src/math/generic/asinf.cpp