Skip to content

Commit bc19aed

Browse files
bassiounixbassiounix
authored
[libc][math] Refactor asin implementation to header-only in src/__support/math folder. (#148578)
Part of #147386 in preparation for: https://discourse.llvm.org/t/rfc-make-clang-builtin-math-functions-constexpr-with-llvm-libc-to-support-c-23-constexpr-math-functions/86450
1 parent 3f3d779 commit bc19aed

File tree

9 files changed

+367
-27
lines changed

9 files changed

+367
-27
lines changed

‎libc/shared/math.h

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -17,6 +17,7 @@
1717
#include "math/acoshf.h"
1818
#include "math/acoshf16.h"
1919
#include "math/acospif16.h"
20+
#include "math/asin.h"
2021
#include "math/erff.h"
2122
#include "math/exp.h"
2223
#include "math/exp10.h"

‎libc/shared/math/asin.h

Lines changed: 23 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,23 @@
1+
//===-- Shared asin function ------------------------------------*- 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_SHARED_MATH_ASIN_H
10+
#define LLVM_LIBC_SHARED_MATH_ASIN_H
11+
12+
#include "shared/libc_common.h"
13+
#include "src/__support/math/asin.h"
14+
15+
namespace LIBC_NAMESPACE_DECL {
16+
namespace shared {
17+
18+
using math::asin;
19+
20+
} // namespace shared
21+
} // namespace LIBC_NAMESPACE_DECL
22+
23+
#endif // LLVM_LIBC_SHARED_MATH_ASIN_H

‎libc/src/__support/math/CMakeLists.txt

Lines changed: 17 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -4,7 +4,6 @@ add_header_library(
44
acos.h
55
DEPENDS
66
.asin_utils
7-
libc.src.__support.math.asin_utils
87
libc.src.__support.FPUtil.double_double
98
libc.src.__support.FPUtil.dyadic_float
109
libc.src.__support.FPUtil.fenv_impl
@@ -124,6 +123,23 @@ add_header_library(
124123
libc.src.__support.macros.optimization
125124
)
126125

126+
add_header_library(
127+
asin
128+
HDRS
129+
asin.h
130+
DEPENDS
131+
.asin_utils
132+
libc.src.__support.FPUtil.double_double
133+
libc.src.__support.FPUtil.dyadic_float
134+
libc.src.__support.FPUtil.fenv_impl
135+
libc.src.__support.FPUtil.fp_bits
136+
libc.src.__support.FPUtil.multiply_add
137+
libc.src.__support.FPUtil.polyeval
138+
libc.src.__support.FPUtil.sqrt
139+
libc.src.__support.macros.optimization
140+
libc.src.__support.macros.properties.cpu_features
141+
)
142+
127143
add_header_library(
128144
erff
129145
HDRS

‎libc/src/__support/math/asin.h

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

‎libc/src/math/generic/CMakeLists.txt

Lines changed: 1 addition & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -3993,16 +3993,7 @@ add_entrypoint_object(
39933993
HDRS
39943994
../asin.h
39953995
DEPENDS
3996-
libc.src.__support.math.asin_utils
3997-
libc.src.__support.FPUtil.double_double
3998-
libc.src.__support.FPUtil.dyadic_float
3999-
libc.src.__support.FPUtil.fenv_impl
4000-
libc.src.__support.FPUtil.fp_bits
4001-
libc.src.__support.FPUtil.multiply_add
4002-
libc.src.__support.FPUtil.polyeval
4003-
libc.src.__support.FPUtil.sqrt
4004-
libc.src.__support.macros.optimization
4005-
libc.src.__support.macros.properties.cpu_features
3996+
libc.src.__support.math.asin
40063997
)
40073998

40083999
add_entrypoint_object(

0 commit comments

Comments
 (0)