1//===-- Implementation header for asinpi ------------------------*- 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_ASINPI_H
10#define LLVM_LIBC_SRC___SUPPORT_MATH_ASINPI_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 asinpi(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 // asinpi(+-0) = +-0.
40 if (LIBC_UNLIKELY(xbits.abs().uintval() == 0))
41 return x;
42 // When |x| < 2^-26, asinpi(x) ~ x/pi.
43 // The relative error of x/pi is:
44 // |asinpi(x) - x/pi| / |asinpi(x)| < x^2/6 < 2^-54.
45#ifdef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
46 return x * asin_internal::ASINPIF_COEFFS[0];
47#endif // LIBC_MATH_HAS_SKIP_ACCURATE_PASS
48 }
49
50#ifdef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
51 double xsq = x * x;
52 return x * fputil::multiply_add(xsq, asin_internal::asinpi_eval(xsq),
53 asin_internal::ASINPIF_COEFFS[0]);
54#else
55 using DFloat128 = fputil::DyadicFloat<128>;
56 using DoubleDouble = fputil::DoubleDouble;
57
58 // For |x| < 2^-511, x^2 would underflow to subnormal, raising a
59 // spurious underflow exception. Since asinpi(x) = x/pi with correction
60 // x^2/(6*pi) < 2^-1024 relative (negligible), compute x/pi directly
61 // in DFloat128.
62 if (LIBC_UNLIKELY(x_exp < 512)) {
63 DFloat128 x_f128(x);
64 DFloat128 r = fputil::quick_mul(a: x_f128, b: ONE_OVER_PI_F128);
65 double result = static_cast<double>(r);
66
67 // IEEE 754 "after rounding" tininess: the 53-bit unlimited-exponent
68 // result is strictly between +-2^-1022. DyadicFloat's conversion
69 // checks the *IEEE subnormal* result (52-bit at the boundary), not
70 // the 53-bit unlimited-exponent result, so we detect it here.
71 int exp_hi = r.exponent + 127 + FPBits::EXP_BIAS;
72 if (LIBC_UNLIKELY(exp_hi <= 0) && !r.mantissa.is_zero()) {
73 bool raise_underflow = true;
74 // When exp_hi == 0, a carry in 53-bit rounding can push the
75 // result to exactly 2^-1022 (not tiny). Check for this.
76 if (exp_hi == 0) {
77 constexpr unsigned SHIFT_53 = 128 - FPBits::SIG_LEN - 1;
78 using MantT = typename DFloat128::MantissaType;
79 MantT m53 = r.mantissa >> SHIFT_53;
80 constexpr MantT ALL_ONES_53 = (MantT(1) << (FPBits::SIG_LEN + 1)) - 1;
81 if (m53 == ALL_ONES_53) {
82 // All 53 bits set. carry happens if rounding rounds away
83 // from zero at this precision.
84 bool round_bit =
85 static_cast<bool>((r.mantissa >> (SHIFT_53 - 1)) & 1);
86 MantT sticky_mask = (MantT(1) << (SHIFT_53 - 1)) - 1;
87 bool sticky = (r.mantissa & sticky_mask) != 0;
88 bool lsb = static_cast<bool>(m53 & 1);
89#ifdef LIBC_MATH_HAS_ASSUME_ROUND_NEAREST_ONLY
90 // Carry if round_bit && (lsb || sticky) (round half to even).
91 raise_underflow = !(round_bit && (lsb || sticky));
92#else
93 switch (fputil::quick_get_round()) {
94 case FE_TONEAREST:
95 // Carry if round_bit && (lsb || sticky) (round half to even).
96 raise_underflow = !(round_bit && (lsb || sticky));
97 break;
98 case FE_UPWARD:
99 raise_underflow = xbits.is_neg() || !(round_bit || sticky);
100 break;
101 case FE_DOWNWARD:
102 raise_underflow = !xbits.is_neg() || !(round_bit || sticky);
103 break;
104 case FE_TOWARDZERO:
105 default:
106 raise_underflow = true; // truncation never carries
107 break;
108 }
109#endif // LIBC_MATH_HAS_ASSUME_ROUND_NEAREST_ONLY
110 }
111 }
112 if (raise_underflow)
113 fputil::raise_except_if_required(FE_UNDERFLOW | FE_INEXACT);
114 }
115 return result;
116 }
117
118 unsigned idx = 0;
119 DoubleDouble x_sq = fputil::exact_mult(a: x, b: x);
120 double err = xbits.abs().get_val() * 0x1.0p-51;
121 // Polynomial approximation:
122 // p ~ asin(x)/(pi*x)
123
124 DoubleDouble p = asinpi_eval(u: x_sq, idx, err);
125 // asinpi(x) ~ x * p
126 DoubleDouble r0 = fputil::exact_mult(a: x, b: p.hi);
127 double r_lo = fputil::multiply_add(x, y: p.lo, z: r0.lo);
128
129 // Ziv's accuracy test.
130 double r_upper = r0.hi + (r_lo + err);
131 double r_lower = r0.hi + (r_lo - err);
132
133 if (LIBC_LIKELY(r_upper == r_lower))
134 return r_upper;
135
136 // Ziv's accuracy test failed, perform 128-bit calculation.
137
138 // Recalculate mod 1/64.
139 idx = static_cast<unsigned>(fputil::nearest_integer(x: x_sq.hi * 0x1.0p6));
140
141 DFloat128 x_f128(x);
142
143#ifdef LIBC_TARGET_CPU_HAS_FMA_DOUBLE
144 DFloat128 u_hi(
145 fputil::multiply_add(static_cast<double>(idx), -0x1.0p-6, x_sq.hi));
146 DFloat128 u = fputil::quick_add(u_hi, DFloat128(x_sq.lo));
147#else
148 DFloat128 x_sq_f128 = fputil::quick_mul(a: x_f128, b: x_f128);
149 DFloat128 u = fputil::quick_add(
150 a: x_sq_f128, b: DFloat128(static_cast<double>(idx) * (-0x1.0p-6)));
151#endif // LIBC_TARGET_CPU_HAS_FMA_DOUBLE
152
153 DFloat128 p_f128 = asinpi_eval(u, idx);
154 DFloat128 r = fputil::quick_mul(a: x_f128, b: p_f128);
155
156 return static_cast<double>(r);
157#endif // LIBC_MATH_HAS_SKIP_ACCURATE_PASS
158 }
159 // |x| >= 0.5
160
161 double x_abs = xbits.abs().get_val();
162
163 // Maintaining the sign:
164 constexpr double SIGN[2] = {1.0, -1.0};
165 double x_sign = SIGN[xbits.is_neg()];
166
167 // |x| >= 1
168 if (LIBC_UNLIKELY(x_exp >= FPBits::EXP_BIAS)) {
169 // x = +-1, asinpi(x) = +- 0.5
170 if (x_abs == 1.0) {
171 return x_sign * 0.5;
172 }
173 // |x| > 1, return NaN.
174 if (xbits.is_quiet_nan())
175 return x;
176
177 // Set domain error for non-NaN input.
178 if (!xbits.is_nan())
179 fputil::set_errno_if_required(EDOM);
180
181 fputil::raise_except_if_required(FE_INVALID);
182 return FPBits::quiet_nan().get_val();
183 }
184
185 // When |x| >= 0.5, we perform range reduction as follow:
186 //
187 // Assume further that 0.5 <= x < 1, and let:
188 // y = asin(x)
189 // Using the identity:
190 // asin(x) = pi/2 - 2 * asin( sqrt( (1 - x)/2 ) )
191 // We get:
192 // asinpi(x) = asin(x)/pi = 0.5 - 2 * asin(sqrt(u)) / pi
193 // = 0.5 - 2 * sqrt(u) * [asin(sqrt(u)) / (pi * sqrt(u))]
194 // = 0.5 - 2 * sqrt(u) * asinpi_eval(u)
195 // where u = (1 - |x|) / 2.
196
197 // u = (1 - |x|)/2
198 double u = fputil::multiply_add(x: x_abs, y: -0.5, z: 0.5);
199 // v_hi ~ sqrt(u).
200 double v_hi = fputil::sqrt<double>(x: u);
201
202#ifdef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
203 double neg2_v = -2.0 * v_hi;
204 double r = x_sign * fputil::multiply_add(
205 neg2_v * u, asin_internal::asinpi_eval(u),
206 0.5 + neg2_v * asin_internal::ASINPIF_COEFFS[0]);
207 return r;
208#else
209 using DFloat128 = fputil::DyadicFloat<128>;
210 using DoubleDouble = fputil::DoubleDouble;
211
212#ifdef LIBC_TARGET_CPU_HAS_FMA_DOUBLE
213 double h = fputil::multiply_add(v_hi, -v_hi, u);
214#else
215 DoubleDouble v_hi_sq = fputil::exact_mult(a: v_hi, b: v_hi);
216 double h = (u - v_hi_sq.hi) - v_hi_sq.lo;
217#endif // LIBC_TARGET_CPU_HAS_FMA_DOUBLE
218
219 // Scale v_lo and v_hi by 2 from the formula:
220 // vh = v_hi * 2
221 // vl = 2*v_lo = h / v_hi.
222 double vh = v_hi * 2.0;
223 double vl = h / v_hi;
224
225 // Polynomial approximation:
226 // p ~ asin(sqrt(u))/(pi*sqrt(u))
227 unsigned idx = 0;
228 double err = vh * 0x1.0p-51;
229
230 DoubleDouble p = asinpi_eval(u: DoubleDouble{.lo: 0.0, .hi: u}, idx, err);
231
232 // Perform computations in double-double arithmetic:
233 // asinpi(x) = 0.5 - (vh + vl) * p
234 DoubleDouble r0 = fputil::quick_mult(a: DoubleDouble{.lo: vl, .hi: vh}, b: p);
235 DoubleDouble r = fputil::exact_add(a: 0.5, b: -r0.hi);
236
237 double r_lo = -r0.lo + r.lo;
238
239 // Ziv's accuracy test.
240
241#ifdef LIBC_TARGET_CPU_HAS_FMA_DOUBLE
242 double r_upper = fputil::multiply_add(
243 r.hi, x_sign, fputil::multiply_add(r_lo, x_sign, err));
244 double r_lower = fputil::multiply_add(
245 r.hi, x_sign, fputil::multiply_add(r_lo, x_sign, -err));
246#else
247 r_lo *= x_sign;
248 r.hi *= x_sign;
249 double r_upper = r.hi + (r_lo + err);
250 double r_lower = r.hi + (r_lo - err);
251#endif // LIBC_TARGET_CPU_HAS_FMA_DOUBLE
252
253 if (LIBC_LIKELY(r_upper == r_lower))
254 return r_upper;
255
256 // Ziv's accuracy test failed, we redo the computations in DFloat128.
257 // Recalculate mod 1/64.
258 idx = static_cast<unsigned>(fputil::nearest_integer(x: u * 0x1.0p6));
259
260 // After the first step of Newton-Raphson approximating v = sqrt(u):
261 // sqrt(u) = v_hi + h / (sqrt(u) + v_hi)
262 // v_lo = h / (2 * v_hi)
263 // Add second-order correction:
264 // v_ll = -v_lo * (h / (4u))
265
266 // Get the rounding error of vl = 2 * v_lo ~ h / vh
267#ifdef LIBC_TARGET_CPU_HAS_FMA_DOUBLE
268 double vl_lo = fputil::multiply_add(-v_hi, vl, h) / v_hi;
269#else
270 DoubleDouble vh_vl = fputil::exact_mult(a: v_hi, b: vl);
271 double vl_lo = ((h - vh_vl.hi) - vh_vl.lo) / v_hi;
272#endif // LIBC_TARGET_CPU_HAS_FMA_DOUBLE
273 // vll = 2*v_ll = -vl * (h / (4u)).
274 double t = h * (-0.25) / u;
275 double vll = fputil::multiply_add(x: vl, y: t, z: vl_lo);
276 // m_v = -(v_hi + v_lo + v_ll).
277 DFloat128 m_v = fputil::quick_add(
278 a: DFloat128(vh), b: fputil::quick_add(a: DFloat128(vl), b: DFloat128(vll)));
279 m_v.sign = Sign::NEG;
280
281 // Perform computations in DFloat128:
282 // asinpi(x) = 0.5 - (v_hi + v_lo + vll) * P_pi(u).
283 DFloat128 y_f128(
284 fputil::multiply_add(x: static_cast<double>(idx), y: -0x1.0p-6, z: u));
285
286 DFloat128 p_f128 = asinpi_eval(u: y_f128, idx);
287 DFloat128 r0_f128 = fputil::quick_mul(a: m_v, b: p_f128);
288 DFloat128 r_f128 = fputil::quick_add(a: HALF_F128, b: r0_f128);
289
290 if (xbits.is_neg())
291 r_f128.sign = Sign::NEG;
292
293 return static_cast<double>(r_f128);
294#endif // LIBC_MATH_HAS_SKIP_ACCURATE_PASS
295}
296
297} // namespace math
298
299} // namespace LIBC_NAMESPACE_DECL
300
301#endif // LLVM_LIBC_SRC___SUPPORT_MATH_ASINPI_H
302