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