1//===-- Implementation header for asin --------------------------*- 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-exception
6//
7//===----------------------------------------------------------------------===//
8
9#ifndef LLVM_LIBC_SRC___SUPPORT_MATH_ASIN_H
10#define LLVM_LIBC_SRC___SUPPORT_MATH_ASIN_H
11
12#include "asin_utils.h"
13#include "src/__support/FPUtil/FEnvImpl.h"
14#include "src/__support/FPUtil/FPBits.h"
15#include "src/__support/FPUtil/double_double.h"
16#include "src/__support/FPUtil/dyadic_float.h"
17#include "src/__support/FPUtil/multiply_add.h"
18#include "src/__support/FPUtil/sqrt.h"
19#include "src/__support/macros/config.h"
20#include "src/__support/macros/optimization.h" // LIBC_UNLIKELY
21#include "src/__support/macros/properties/cpu_features.h" // LIBC_TARGET_CPU_HAS_FMA
22#include "src/__support/math/asin_utils.h"
23
24namespace LIBC_NAMESPACE_DECL {
25
26namespace math {
27
28LIBC_INLINE double asin(double x) {
29 using namespace asin_internal;
30 using FPBits = fputil::FPBits<double>;
31
32 FPBits xbits(x);
33 int x_exp = xbits.get_biased_exponent();
34
35 // |x| < 0.5.
36 if (x_exp < FPBits::EXP_BIAS - 1) {
37 // |x| < 2^-26.
38 if (LIBC_UNLIKELY(x_exp < FPBits::EXP_BIAS - 26)) {
39 // When |x| < 2^-26, the relative error of the approximation asin(x) ~ x
40 // is:
41 // |asin(x) - x| / |asin(x)| < |x^3| / (6|x|)
42 // = x^2 / 6
43 // < 2^-54
44 // < epsilon(1)/2.
45 // So the correctly rounded values of asin(x) are:
46 // = x + sign(x)*eps(x) if rounding mode = FE_TOWARDZERO,
47 // or (rounding mode = FE_UPWARD and x is
48 // negative),
49 // = x otherwise.
50 // To simplify the rounding decision and make it more efficient, we use
51 // fma(x, 2^-54, x) instead.
52 // Note: to use the formula x + 2^-54*x to decide the correct rounding, we
53 // do need fma(x, 2^-54, x) to prevent underflow caused by 2^-54*x when
54 // |x| < 2^-1022. For targets without FMA instructions, when x is close to
55 // denormal range, we normalize x,
56#if defined(LIBC_MATH_HAS_SKIP_ACCURATE_PASS)
57 return x;
58#elif defined(LIBC_TARGET_CPU_HAS_FMA_DOUBLE)
59 return fputil::multiply_add(x, 0x1.0p-54, x);
60#else
61 if (xbits.abs().uintval() == 0)
62 return x;
63 // Get sign(x) * min_normal.
64 FPBits eps_bits = FPBits::min_normal();
65 eps_bits.set_sign(xbits.sign());
66 double eps = eps_bits.get_val();
67 double normalize_const = (x_exp == 0) ? eps : 0.0;
68 double scaled_normal =
69 fputil::multiply_add(x: x + normalize_const, y: 0x1.0p54, z: eps);
70 return fputil::multiply_add(x: scaled_normal, y: 0x1.0p-54, z: -normalize_const);
71#endif // LIBC_MATH_HAS_SKIP_ACCURATE_PASS
72 }
73
74#ifdef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
75 double xsq = x * x;
76 return fputil::multiply_add(x * xsq, asin_eval(xsq), x);
77#else
78 using DFloat128 = fputil::DyadicFloat<128>;
79 using DoubleDouble = fputil::DoubleDouble;
80
81 unsigned idx = 0;
82 DoubleDouble x_sq = fputil::exact_mult(a: x, b: x);
83 double err = xbits.abs().get_val() * 0x1.0p-51;
84 // Polynomial approximation:
85 // p ~ asin(x)/x
86
87 DoubleDouble p = asin_eval(u: x_sq, idx, err);
88 // asin(x) ~ x * (ASIN_COEFFS[idx][0] + p)
89 DoubleDouble r0 = fputil::exact_mult(a: x, b: p.hi);
90 double r_lo = fputil::multiply_add(x, y: p.lo, z: r0.lo);
91
92 // Ziv's accuracy test.
93
94 double r_upper = r0.hi + (r_lo + err);
95 double r_lower = r0.hi + (r_lo - err);
96
97 if (LIBC_LIKELY(r_upper == r_lower))
98 return r_upper;
99
100 // Ziv's accuracy test failed, perform 128-bit calculation.
101
102 // Recalculate mod 1/64.
103 idx = static_cast<unsigned>(fputil::nearest_integer(x: x_sq.hi * 0x1.0p6));
104
105 // Get x^2 - idx/64 exactly. When FMA is available, double-double
106 // multiplication will be correct for all rounding modes. Otherwise we use
107 // DFloat128 directly.
108 DFloat128 x_f128(x);
109
110#ifdef LIBC_TARGET_CPU_HAS_FMA_DOUBLE
111 // u = x^2 - idx/64
112 DFloat128 u_hi(
113 fputil::multiply_add(static_cast<double>(idx), -0x1.0p-6, x_sq.hi));
114 DFloat128 u = fputil::quick_add(u_hi, DFloat128(x_sq.lo));
115#else
116 DFloat128 x_sq_f128 = fputil::quick_mul(a: x_f128, b: x_f128);
117 DFloat128 u = fputil::quick_add(
118 a: x_sq_f128, b: DFloat128(static_cast<double>(idx) * (-0x1.0p-6)));
119#endif // LIBC_TARGET_CPU_HAS_FMA_DOUBLE
120
121 DFloat128 p_f128 = asin_eval(u, idx);
122 DFloat128 r = fputil::quick_mul(a: x_f128, b: p_f128);
123
124 return static_cast<double>(r);
125#endif // LIBC_MATH_HAS_SKIP_ACCURATE_PASS
126 }
127 // |x| >= 0.5
128
129 double x_abs = xbits.abs().get_val();
130
131 // Maintaining the sign:
132 constexpr double SIGN[2] = {1.0, -1.0};
133 double x_sign = SIGN[xbits.is_neg()];
134
135 // |x| >= 1
136 if (LIBC_UNLIKELY(x_exp >= FPBits::EXP_BIAS)) {
137 // x = +-1, asin(x) = +- pi/2
138 if (x_abs == 1.0) {
139 // return +- pi/2
140 return fputil::multiply_add(x: x_sign, y: PI_OVER_TWO.hi,
141 z: x_sign * PI_OVER_TWO.lo);
142 }
143 // |x| > 1, return NaN.
144 if (xbits.is_quiet_nan())
145 return x;
146
147 // Set domain error for non-NaN input.
148 if (!xbits.is_nan())
149 fputil::set_errno_if_required(EDOM);
150
151 fputil::raise_except_if_required(FE_INVALID);
152 return FPBits::quiet_nan().get_val();
153 }
154
155 // When |x| >= 0.5, we perform range reduction as follow:
156 //
157 // Assume further that 0.5 <= x < 1, and let:
158 // y = asin(x)
159 // We will use the double angle formula:
160 // cos(2y) = 1 - 2 sin^2(y)
161 // and the complement angle identity:
162 // x = sin(y) = cos(pi/2 - y)
163 // = 1 - 2 sin^2 (pi/4 - y/2)
164 // So:
165 // sin(pi/4 - y/2) = sqrt( (1 - x)/2 )
166 // And hence:
167 // pi/4 - y/2 = asin( sqrt( (1 - x)/2 ) )
168 // Equivalently:
169 // asin(x) = y = pi/2 - 2 * asin( sqrt( (1 - x)/2 ) )
170 // Let u = (1 - x)/2, then:
171 // asin(x) = pi/2 - 2 * asin( sqrt(u) )
172 // Moreover, since 0.5 <= x < 1:
173 // 0 < u <= 1/4, and 0 < sqrt(u) <= 0.5,
174 // And hence we can reuse the same polynomial approximation of asin(x) when
175 // |x| <= 0.5:
176 // asin(x) ~ pi/2 - 2 * sqrt(u) * P(u),
177
178 // u = (1 - |x|)/2
179 double u = fputil::multiply_add(x: x_abs, y: -0.5, z: 0.5);
180 // v_hi + v_lo ~ sqrt(u).
181 // Let:
182 // h = u - v_hi^2 = (sqrt(u) - v_hi) * (sqrt(u) + v_hi)
183 // Then:
184 // sqrt(u) = v_hi + h / (sqrt(u) + v_hi)
185 // ~ v_hi + h / (2 * v_hi)
186 // So we can use:
187 // v_lo = h / (2 * v_hi).
188 // Then,
189 // asin(x) ~ pi/2 - 2*(v_hi + v_lo) * P(u)
190 double v_hi = fputil::sqrt<double>(x: u);
191
192#ifdef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
193 double neg2_v = -2.0 * v_hi;
194 double r = x_sign * fputil::multiply_add(neg2_v * u, asin_eval(u),
195 PI_OVER_TWO.hi + neg2_v);
196 return r;
197#else
198
199#ifdef LIBC_TARGET_CPU_HAS_FMA_DOUBLE
200 double h = fputil::multiply_add(v_hi, -v_hi, u);
201#else
202 DoubleDouble v_hi_sq = fputil::exact_mult(a: v_hi, b: v_hi);
203 double h = (u - v_hi_sq.hi) - v_hi_sq.lo;
204#endif // LIBC_TARGET_CPU_HAS_FMA_DOUBLE
205
206 // Scale v_lo and v_hi by 2 from the formula:
207 // vh = v_hi * 2
208 // vl = 2*v_lo = h / v_hi.
209 double vh = v_hi * 2.0;
210 double vl = h / v_hi;
211
212 // Polynomial approximation:
213 // p ~ asin(sqrt(u))/sqrt(u)
214 unsigned idx = 0;
215 double err = vh * 0x1.0p-51;
216
217 DoubleDouble p = asin_eval(u: DoubleDouble{.lo: 0.0, .hi: u}, idx, err);
218
219 // Perform computations in double-double arithmetic:
220 // asin(x) = pi/2 - (v_hi + v_lo) * (ASIN_COEFFS[idx][0] + p)
221 DoubleDouble r0 = fputil::quick_mult(a: DoubleDouble{.lo: vl, .hi: vh}, b: p);
222 DoubleDouble r = fputil::exact_add(a: PI_OVER_TWO.hi, b: -r0.hi);
223
224 double r_lo = PI_OVER_TWO.lo - r0.lo + r.lo;
225
226 // Ziv's accuracy test.
227
228#ifdef LIBC_TARGET_CPU_HAS_FMA_DOUBLE
229 double r_upper = fputil::multiply_add(
230 r.hi, x_sign, fputil::multiply_add(r_lo, x_sign, err));
231 double r_lower = fputil::multiply_add(
232 r.hi, x_sign, fputil::multiply_add(r_lo, x_sign, -err));
233#else
234 r_lo *= x_sign;
235 r.hi *= x_sign;
236 double r_upper = r.hi + (r_lo + err);
237 double r_lower = r.hi + (r_lo - err);
238#endif // LIBC_TARGET_CPU_HAS_FMA_DOUBLE
239
240 if (LIBC_LIKELY(r_upper == r_lower))
241 return r_upper;
242
243 // Ziv's accuracy test failed, we redo the computations in DFloat128.
244 // Recalculate mod 1/64.
245 idx = static_cast<unsigned>(fputil::nearest_integer(x: u * 0x1.0p6));
246
247 // After the first step of Newton-Raphson approximating v = sqrt(u), we have
248 // that:
249 // sqrt(u) = v_hi + h / (sqrt(u) + v_hi)
250 // v_lo = h / (2 * v_hi)
251 // With error:
252 // sqrt(u) - (v_hi + v_lo) = h * ( 1/(sqrt(u) + v_hi) - 1/(2*v_hi) )
253 // = -h^2 / (2*v * (sqrt(u) + v)^2).
254 // Since:
255 // (sqrt(u) + v_hi)^2 ~ (2sqrt(u))^2 = 4u,
256 // we can add another correction term to (v_hi + v_lo) that is:
257 // v_ll = -h^2 / (2*v_hi * 4u)
258 // = -v_lo * (h / 4u)
259 // = -vl * (h / 8u),
260 // making the errors:
261 // sqrt(u) - (v_hi + v_lo + v_ll) = O(h^3)
262 // well beyond 128-bit precision needed.
263
264 // Get the rounding error of vl = 2 * v_lo ~ h / vh
265 // Get full product of vh * vl
266#ifdef LIBC_TARGET_CPU_HAS_FMA_DOUBLE
267 double vl_lo = fputil::multiply_add(-v_hi, vl, h) / v_hi;
268#else
269 DoubleDouble vh_vl = fputil::exact_mult(a: v_hi, b: vl);
270 double vl_lo = ((h - vh_vl.hi) - vh_vl.lo) / v_hi;
271#endif // LIBC_TARGET_CPU_HAS_FMA_DOUBLE
272 // vll = 2*v_ll = -vl * (h / (4u)).
273 double t = h * (-0.25) / u;
274 double vll = fputil::multiply_add(x: vl, y: t, z: vl_lo);
275 // m_v = -(v_hi + v_lo + v_ll).
276 DFloat128 m_v = fputil::quick_add(
277 a: DFloat128(vh), b: fputil::quick_add(a: DFloat128(vl), b: DFloat128(vll)));
278 m_v.sign = Sign::NEG;
279
280 // Perform computations in DFloat128:
281 // asin(x) = pi/2 - (v_hi + v_lo + vll) * P(u).
282 DFloat128 y_f128(
283 fputil::multiply_add(x: static_cast<double>(idx), y: -0x1.0p-6, z: u));
284
285 DFloat128 p_f128 = asin_eval(u: y_f128, idx);
286 DFloat128 r0_f128 = fputil::quick_mul(a: m_v, b: p_f128);
287 DFloat128 r_f128 = fputil::quick_add(a: PI_OVER_TWO_F128, b: r0_f128);
288
289 if (xbits.is_neg())
290 r_f128.sign = Sign::NEG;
291
292 return static_cast<double>(r_f128);
293#endif // LIBC_MATH_HAS_SKIP_ACCURATE_PASS
294}
295
296} // namespace math
297
298} // namespace LIBC_NAMESPACE_DECL
299
300#endif // LLVM_LIBC_SRC___SUPPORT_MATH_ASIN_H
301