Skip to content

Commit b95a303

Browse files
bassiounixmahesh-attarde
bassiounix
authored and
mahesh-attarde
committed
[libc][math] Refactor asinf implementation to header-only in src/__support/math folder. (llvm#150697)
Part of llvm#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 076275f commit b95a303

File tree

9 files changed

+237
-168
lines changed

9 files changed

+237
-168
lines changed

‎libc/shared/math.h

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -18,6 +18,7 @@
1818
#include "math/acoshf16.h"
1919
#include "math/acospif16.h"
2020
#include "math/asin.h"
21+
#include "math/asinf.h"
2122
#include "math/erff.h"
2223
#include "math/exp.h"
2324
#include "math/exp10.h"

‎libc/shared/math/asinf.h

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

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

Lines changed: 16 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -140,6 +140,22 @@ add_header_library(
140140
libc.src.__support.macros.properties.cpu_features
141141
)
142142

143+
add_header_library(
144+
asinf
145+
HDRS
146+
asinf.h
147+
DEPENDS
148+
.inv_trigf_utils
149+
libc.src.__support.FPUtil.fenv_impl
150+
libc.src.__support.FPUtil.fp_bits
151+
libc.src.__support.FPUtil.except_value_utils
152+
libc.src.__support.FPUtil.multiply_add
153+
libc.src.__support.FPUtil.sqrt
154+
libc.src.__support.macros.config
155+
libc.src.__support.macros.optimization
156+
libc.src.__support.macros.properties.cpu_features
157+
)
158+
143159
add_header_library(
144160
erff
145161
HDRS

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

Lines changed: 175 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,175 @@
1+
//===-- Implementation header for asinf -------------------------*- 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_ASINF_H
10+
#define LLVM_LIBC_SRC___SUPPORT_MATH_ASINF_H
11+
12+
#include "inv_trigf_utils.h"
13+
#include "src/__support/FPUtil/FEnvImpl.h"
14+
#include "src/__support/FPUtil/FPBits.h"
15+
#include "src/__support/FPUtil/except_value_utils.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+
22+
namespace LIBC_NAMESPACE_DECL {
23+
24+
namespace math {
25+
26+
LIBC_INLINE static constexpr float asinf(float x) {
27+
using namespace inv_trigf_utils_internal;
28+
29+
#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
30+
constexpr size_t N_EXCEPTS = 2;
31+
32+
// Exceptional values when |x| <= 0.5
33+
constexpr fputil::ExceptValues<float, N_EXCEPTS> ASINF_EXCEPTS_LO = {{
34+
// (inputs, RZ output, RU offset, RD offset, RN offset)
35+
// x = 0x1.137f0cp-5, asinf(x) = 0x1.138c58p-5 (RZ)
36+
{0x3d09bf86, 0x3d09c62c, 1, 0, 1},
37+
// x = 0x1.cbf43cp-4, asinf(x) = 0x1.cced1cp-4 (RZ)
38+
{0x3de5fa1e, 0x3de6768e, 1, 0, 0},
39+
}};
40+
41+
// Exceptional values when 0.5 < |x| <= 1
42+
constexpr fputil::ExceptValues<float, N_EXCEPTS> ASINF_EXCEPTS_HI = {{
43+
// (inputs, RZ output, RU offset, RD offset, RN offset)
44+
// x = 0x1.107434p-1, asinf(x) = 0x1.1f4b64p-1 (RZ)
45+
{0x3f083a1a, 0x3f0fa5b2, 1, 0, 0},
46+
// x = 0x1.ee836cp-1, asinf(x) = 0x1.4f0654p0 (RZ)
47+
{0x3f7741b6, 0x3fa7832a, 1, 0, 0},
48+
}};
49+
#endif // !LIBC_MATH_HAS_SKIP_ACCURATE_PASS
50+
51+
using namespace inv_trigf_utils_internal;
52+
using FPBits = typename fputil::FPBits<float>;
53+
54+
FPBits xbits(x);
55+
uint32_t x_uint = xbits.uintval();
56+
uint32_t x_abs = xbits.uintval() & 0x7fff'ffffU;
57+
constexpr double SIGN[2] = {1.0, -1.0};
58+
uint32_t x_sign = x_uint >> 31;
59+
60+
// |x| <= 0.5-ish
61+
if (x_abs < 0x3f04'471dU) {
62+
// |x| < 0x1.d12edp-12
63+
if (LIBC_UNLIKELY(x_abs < 0x39e8'9768U)) {
64+
// When |x| < 2^-12, the relative error of the approximation asin(x) ~ x
65+
// is:
66+
// |asin(x) - x| / |asin(x)| < |x^3| / (6|x|)
67+
// = x^2 / 6
68+
// < 2^-25
69+
// < epsilon(1)/2.
70+
// So the correctly rounded values of asin(x) are:
71+
// = x + sign(x)*eps(x) if rounding mode = FE_TOWARDZERO,
72+
// or (rounding mode = FE_UPWARD and x is
73+
// negative),
74+
// = x otherwise.
75+
// To simplify the rounding decision and make it more efficient, we use
76+
// fma(x, 2^-25, x) instead.
77+
// An exhaustive test shows that this formula work correctly for all
78+
// rounding modes up to |x| < 0x1.d12edp-12.
79+
// Note: to use the formula x + 2^-25*x to decide the correct rounding, we
80+
// do need fma(x, 2^-25, x) to prevent underflow caused by 2^-25*x when
81+
// |x| < 2^-125. For targets without FMA instructions, we simply use
82+
// double for intermediate results as it is more efficient than using an
83+
// emulated version of FMA.
84+
#if defined(LIBC_TARGET_CPU_HAS_FMA_FLOAT)
85+
return fputil::multiply_add(x, 0x1.0p-25f, x);
86+
#else
87+
double xd = static_cast<double>(x);
88+
return static_cast<float>(fputil::multiply_add(xd, 0x1.0p-25, xd));
89+
#endif // LIBC_TARGET_CPU_HAS_FMA_FLOAT
90+
}
91+
92+
#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
93+
// Check for exceptional values
94+
if (auto r = ASINF_EXCEPTS_LO.lookup_odd(x_abs, x_sign);
95+
LIBC_UNLIKELY(r.has_value()))
96+
return r.value();
97+
#endif // !LIBC_MATH_HAS_SKIP_ACCURATE_PASS
98+
99+
// For |x| <= 0.5, we approximate asinf(x) by:
100+
// asin(x) = x * P(x^2)
101+
// Where P(X^2) = Q(X) is a degree-20 minimax even polynomial approximating
102+
// asin(x)/x on [0, 0.5] generated by Sollya with:
103+
// > Q = fpminimax(asin(x)/x, [|0, 2, 4, 6, 8, 10, 12, 14, 16, 18, 20|],
104+
// [|1, D...|], [0, 0.5]);
105+
// An exhaustive test shows that this approximation works well up to a
106+
// little more than 0.5.
107+
double xd = static_cast<double>(x);
108+
double xsq = xd * xd;
109+
double x3 = xd * xsq;
110+
double r = asin_eval(xsq);
111+
return static_cast<float>(fputil::multiply_add(x3, r, xd));
112+
}
113+
114+
// |x| > 1, return NaNs.
115+
if (LIBC_UNLIKELY(x_abs > 0x3f80'0000U)) {
116+
if (xbits.is_signaling_nan()) {
117+
fputil::raise_except_if_required(FE_INVALID);
118+
return FPBits::quiet_nan().get_val();
119+
}
120+
121+
if (x_abs <= 0x7f80'0000U) {
122+
fputil::set_errno_if_required(EDOM);
123+
fputil::raise_except_if_required(FE_INVALID);
124+
}
125+
126+
return FPBits::quiet_nan().get_val();
127+
}
128+
129+
#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
130+
// Check for exceptional values
131+
if (auto r = ASINF_EXCEPTS_HI.lookup_odd(x_abs, x_sign);
132+
LIBC_UNLIKELY(r.has_value()))
133+
return r.value();
134+
#endif // !LIBC_MATH_HAS_SKIP_ACCURATE_PASS
135+
136+
// When |x| > 0.5, we perform range reduction as follow:
137+
//
138+
// Assume further that 0.5 < x <= 1, and let:
139+
// y = asin(x)
140+
// We will use the double angle formula:
141+
// cos(2y) = 1 - 2 sin^2(y)
142+
// and the complement angle identity:
143+
// x = sin(y) = cos(pi/2 - y)
144+
// = 1 - 2 sin^2 (pi/4 - y/2)
145+
// So:
146+
// sin(pi/4 - y/2) = sqrt( (1 - x)/2 )
147+
// And hence:
148+
// pi/4 - y/2 = asin( sqrt( (1 - x)/2 ) )
149+
// Equivalently:
150+
// asin(x) = y = pi/2 - 2 * asin( sqrt( (1 - x)/2 ) )
151+
// Let u = (1 - x)/2, then:
152+
// asin(x) = pi/2 - 2 * asin( sqrt(u) )
153+
// Moreover, since 0.5 < x <= 1:
154+
// 0 <= u < 1/4, and 0 <= sqrt(u) < 0.5,
155+
// And hence we can reuse the same polynomial approximation of asin(x) when
156+
// |x| <= 0.5:
157+
// asin(x) ~ pi/2 - 2 * sqrt(u) * P(u),
158+
159+
xbits.set_sign(Sign::POS);
160+
double sign = SIGN[x_sign];
161+
double xd = static_cast<double>(xbits.get_val());
162+
double u = fputil::multiply_add(-0.5, xd, 0.5);
163+
double c1 = sign * (-2 * fputil::sqrt<double>(u));
164+
double c2 = fputil::multiply_add(sign, M_MATH_PI_2, c1);
165+
double c3 = c1 * u;
166+
167+
double r = asin_eval(u);
168+
return static_cast<float>(fputil::multiply_add(c3, r, c2));
169+
}
170+
171+
} // namespace math
172+
173+
} // namespace LIBC_NAMESPACE_DECL
174+
175+
#endif // LLVM_LIBC_SRC___SUPPORT_MATH_ASINF_H

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

Lines changed: 1 addition & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -3958,13 +3958,7 @@ add_entrypoint_object(
39583958
HDRS
39593959
../asinf.h
39603960
DEPENDS
3961-
libc.src.__support.FPUtil.except_value_utils
3962-
libc.src.__support.FPUtil.fp_bits
3963-
libc.src.__support.FPUtil.multiply_add
3964-
libc.src.__support.FPUtil.polyeval
3965-
libc.src.__support.FPUtil.sqrt
3966-
libc.src.__support.macros.optimization
3967-
libc.src.__support.math.inv_trigf_utils
3961+
libc.src.__support.math.asinf
39683962
)
39693963

39703964
add_entrypoint_object(

0 commit comments

Comments
 (0)