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

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