Skip to content

Commit 85213f2

Browse files
bassiounixbassiounix
authored
[libc][math] Refactor asinf16 implementation to header-only in src/__support/math folder. (#150800)
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 f3bcaea commit 85213f2

File tree

9 files changed

+208
-137
lines changed

9 files changed

+208
-137
lines changed

‎libc/shared/math.h

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -19,6 +19,7 @@
1919
#include "math/acospif16.h"
2020
#include "math/asin.h"
2121
#include "math/asinf.h"
22+
#include "math/asinf16.h"
2223
#include "math/erff.h"
2324
#include "math/exp.h"
2425
#include "math/exp10.h"

‎libc/shared/math/asinf16.h

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

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

Lines changed: 14 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -156,6 +156,20 @@ add_header_library(
156156
libc.src.__support.macros.properties.cpu_features
157157
)
158158

159+
add_header_library(
160+
asinf16
161+
HDRS
162+
asinf16.h
163+
DEPENDS
164+
libc.src.__support.FPUtil.fenv_impl
165+
libc.src.__support.FPUtil.fp_bits
166+
libc.src.__support.FPUtil.polyeval
167+
libc.src.__support.FPUtil.cast
168+
libc.src.__support.FPUtil.multiply_add
169+
libc.src.__support.FPUtil.sqrt
170+
libc.src.__support.macros.optimization
171+
)
172+
159173
add_header_library(
160174
erff
161175
HDRS

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

Lines changed: 146 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,146 @@
1+
//===-- Implementation header for asinf16 -----------------------*- 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_ASINF16_H
10+
#define LLVM_LIBC_SRC___SUPPORT_MATH_ASINF16_H
11+
12+
#include "include/llvm-libc-macros/float16-macros.h"
13+
14+
#ifdef LIBC_TYPES_HAS_FLOAT16
15+
16+
#include "src/__support/FPUtil/FEnvImpl.h"
17+
#include "src/__support/FPUtil/FPBits.h"
18+
#include "src/__support/FPUtil/PolyEval.h"
19+
#include "src/__support/FPUtil/cast.h"
20+
#include "src/__support/FPUtil/multiply_add.h"
21+
#include "src/__support/FPUtil/sqrt.h"
22+
#include "src/__support/macros/optimization.h"
23+
24+
namespace LIBC_NAMESPACE_DECL {
25+
26+
namespace math {
27+
28+
LIBC_INLINE static constexpr float16 asinf16(float16 x) {
29+
30+
// Generated by Sollya using the following command:
31+
// > round(pi/2, D, RN);
32+
constexpr float PI_2 = 0x1.921fb54442d18p0f;
33+
34+
using FPBits = fputil::FPBits<float16>;
35+
FPBits xbits(x);
36+
37+
uint16_t x_u = xbits.uintval();
38+
uint16_t x_abs = x_u & 0x7fff;
39+
float xf = x;
40+
41+
// |x| > 0x1p0, |x| > 1, or x is NaN.
42+
if (LIBC_UNLIKELY(x_abs > 0x3c00)) {
43+
// asinf16(NaN) = NaN
44+
if (xbits.is_nan()) {
45+
if (xbits.is_signaling_nan()) {
46+
fputil::raise_except_if_required(FE_INVALID);
47+
return FPBits::quiet_nan().get_val();
48+
}
49+
50+
return x;
51+
}
52+
53+
// 1 < |x| <= +/-inf
54+
fputil::raise_except_if_required(FE_INVALID);
55+
fputil::set_errno_if_required(EDOM);
56+
57+
return FPBits::quiet_nan().get_val();
58+
}
59+
60+
float xsq = xf * xf;
61+
62+
// |x| <= 0x1p-1, |x| <= 0.5
63+
if (x_abs <= 0x3800) {
64+
// asinf16(+/-0) = +/-0
65+
if (LIBC_UNLIKELY(x_abs == 0))
66+
return x;
67+
68+
// Exhaustive tests show that,
69+
// for |x| <= 0x1.878p-9, when:
70+
// x > 0, and rounding upward, or
71+
// x < 0, and rounding downward, then,
72+
// asin(x) = x * 2^-11 + x
73+
// else, in other rounding modes,
74+
// asin(x) = x
75+
if (LIBC_UNLIKELY(x_abs <= 0x1a1e)) {
76+
int rounding = fputil::quick_get_round();
77+
78+
if ((xbits.is_pos() && rounding == FE_UPWARD) ||
79+
(xbits.is_neg() && rounding == FE_DOWNWARD))
80+
return fputil::cast<float16>(fputil::multiply_add(xf, 0x1.0p-11f, xf));
81+
return x;
82+
}
83+
84+
// Degree-6 minimax odd polynomial of asin(x) generated by Sollya with:
85+
// > P = fpminimax(asin(x)/x, [|0, 2, 4, 6, 8|], [|SG...|], [0, 0.5]);
86+
float result =
87+
fputil::polyeval(xsq, 0x1.000002p0f, 0x1.554c2ap-3f, 0x1.3541ccp-4f,
88+
0x1.43b2d6p-5f, 0x1.a0d73ep-5f);
89+
return fputil::cast<float16>(xf * result);
90+
}
91+
92+
// When |x| > 0.5, assume that 0.5 < |x| <= 1,
93+
//
94+
// Step-by-step range-reduction proof:
95+
// 1: Let y = asin(x), such that, x = sin(y)
96+
// 2: From complimentary angle identity:
97+
// x = sin(y) = cos(pi/2 - y)
98+
// 3: Let z = pi/2 - y, such that x = cos(z)
99+
// 4: From double angle formula; cos(2A) = 1 - sin^2(A):
100+
// z = 2A, z/2 = A
101+
// cos(z) = 1 - 2 * sin^2(z/2)
102+
// 5: Make sin(z/2) subject of the formula:
103+
// sin(z/2) = sqrt((1 - cos(z))/2)
104+
// 6: Recall [3]; x = cos(z). Therefore:
105+
// sin(z/2) = sqrt((1 - x)/2)
106+
// 7: Let u = (1 - x)/2
107+
// 8: Therefore:
108+
// asin(sqrt(u)) = z/2
109+
// 2 * asin(sqrt(u)) = z
110+
// 9: Recall [3], z = pi/2 - y. Therefore:
111+
// y = pi/2 - z
112+
// y = pi/2 - 2 * asin(sqrt(u))
113+
// 10: Recall [1], y = asin(x). Therefore:
114+
// asin(x) = pi/2 - 2 * asin(sqrt(u))
115+
//
116+
// WHY?
117+
// 11: Recall [7], u = (1 - x)/2
118+
// 12: Since 0.5 < x <= 1, therefore:
119+
// 0 <= u <= 0.25 and 0 <= sqrt(u) <= 0.5
120+
//
121+
// Hence, we can reuse the same [0, 0.5] domain polynomial approximation for
122+
// Step [10] as `sqrt(u)` is in range.
123+
124+
// 0x1p-1 < |x| <= 0x1p0, 0.5 < |x| <= 1.0
125+
float xf_abs = (xf < 0 ? -xf : xf);
126+
float sign = (xbits.uintval() >> 15 == 1 ? -1.0 : 1.0);
127+
float u = fputil::multiply_add(-0.5f, xf_abs, 0.5f);
128+
float u_sqrt = fputil::sqrt<float>(u);
129+
130+
// Degree-6 minimax odd polynomial of asin(x) generated by Sollya with:
131+
// > P = fpminimax(asin(x)/x, [|0, 2, 4, 6, 8|], [|SG...|], [0, 0.5]);
132+
float asin_sqrt_u =
133+
u_sqrt * fputil::polyeval(u, 0x1.000002p0f, 0x1.554c2ap-3f,
134+
0x1.3541ccp-4f, 0x1.43b2d6p-5f, 0x1.a0d73ep-5f);
135+
136+
return fputil::cast<float16>(sign *
137+
fputil::multiply_add(-2.0f, asin_sqrt_u, PI_2));
138+
}
139+
140+
} // namespace math
141+
142+
} // namespace LIBC_NAMESPACE_DECL
143+
144+
#endif // LIBC_TYPES_HAS_FLOAT16
145+
146+
#endif // LLVM_LIBC_SRC___SUPPORT_MATH_ASINF16_H

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

Lines changed: 1 addition & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -3968,16 +3968,7 @@ add_entrypoint_object(
39683968
HDRS
39693969
../asinf16.h
39703970
DEPENDS
3971-
libc.hdr.errno_macros
3972-
libc.hdr.fenv_macros
3973-
libc.src.__support.FPUtil.cast
3974-
libc.src.__support.FPUtil.fenv_impl
3975-
libc.src.__support.FPUtil.fp_bits
3976-
libc.src.__support.FPUtil.multiply_add
3977-
libc.src.__support.FPUtil.polyeval
3978-
libc.src.__support.FPUtil.sqrt
3979-
libc.src.__support.macros.optimization
3980-
libc.src.__support.macros.properties.types
3971+
libc.src.__support.math.asinf16
39813972
)
39823973

39833974
add_entrypoint_object(

‎libc/src/math/generic/asinf16.cpp

Lines changed: 2 additions & 119 deletions
Original file line numberDiff line numberDiff line change
@@ -7,127 +7,10 @@
77
//===----------------------------------------------------------------------===//
88

99
#include "src/math/asinf16.h"
10-
#include "hdr/errno_macros.h"
11-
#include "hdr/fenv_macros.h"
12-
#include "src/__support/FPUtil/FEnvImpl.h"
13-
#include "src/__support/FPUtil/FPBits.h"
14-
#include "src/__support/FPUtil/PolyEval.h"
15-
#include "src/__support/FPUtil/cast.h"
16-
#include "src/__support/FPUtil/multiply_add.h"
17-
#include "src/__support/FPUtil/sqrt.h"
18-
#include "src/__support/macros/optimization.h"
10+
#include "src/__support/math/asinf16.h"
1911

2012
namespace LIBC_NAMESPACE_DECL {
2113

22-
// Generated by Sollya using the following command:
23-
// > round(pi/2, D, RN);
24-
static constexpr float PI_2 = 0x1.921fb54442d18p0f;
25-
26-
LLVM_LIBC_FUNCTION(float16, asinf16, (float16 x)) {
27-
using FPBits = fputil::FPBits<float16>;
28-
FPBits xbits(x);
29-
30-
uint16_t x_u = xbits.uintval();
31-
uint16_t x_abs = x_u & 0x7fff;
32-
float xf = x;
33-
34-
// |x| > 0x1p0, |x| > 1, or x is NaN.
35-
if (LIBC_UNLIKELY(x_abs > 0x3c00)) {
36-
// asinf16(NaN) = NaN
37-
if (xbits.is_nan()) {
38-
if (xbits.is_signaling_nan()) {
39-
fputil::raise_except_if_required(FE_INVALID);
40-
return FPBits::quiet_nan().get_val();
41-
}
42-
43-
return x;
44-
}
45-
46-
// 1 < |x| <= +/-inf
47-
fputil::raise_except_if_required(FE_INVALID);
48-
fputil::set_errno_if_required(EDOM);
49-
50-
return FPBits::quiet_nan().get_val();
51-
}
52-
53-
float xsq = xf * xf;
54-
55-
// |x| <= 0x1p-1, |x| <= 0.5
56-
if (x_abs <= 0x3800) {
57-
// asinf16(+/-0) = +/-0
58-
if (LIBC_UNLIKELY(x_abs == 0))
59-
return x;
60-
61-
// Exhaustive tests show that,
62-
// for |x| <= 0x1.878p-9, when:
63-
// x > 0, and rounding upward, or
64-
// x < 0, and rounding downward, then,
65-
// asin(x) = x * 2^-11 + x
66-
// else, in other rounding modes,
67-
// asin(x) = x
68-
if (LIBC_UNLIKELY(x_abs <= 0x1a1e)) {
69-
int rounding = fputil::quick_get_round();
70-
71-
if ((xbits.is_pos() && rounding == FE_UPWARD) ||
72-
(xbits.is_neg() && rounding == FE_DOWNWARD))
73-
return fputil::cast<float16>(fputil::multiply_add(xf, 0x1.0p-11f, xf));
74-
return x;
75-
}
76-
77-
// Degree-6 minimax odd polynomial of asin(x) generated by Sollya with:
78-
// > P = fpminimax(asin(x)/x, [|0, 2, 4, 6, 8|], [|SG...|], [0, 0.5]);
79-
float result =
80-
fputil::polyeval(xsq, 0x1.000002p0f, 0x1.554c2ap-3f, 0x1.3541ccp-4f,
81-
0x1.43b2d6p-5f, 0x1.a0d73ep-5f);
82-
return fputil::cast<float16>(xf * result);
83-
}
84-
85-
// When |x| > 0.5, assume that 0.5 < |x| <= 1,
86-
//
87-
// Step-by-step range-reduction proof:
88-
// 1: Let y = asin(x), such that, x = sin(y)
89-
// 2: From complimentary angle identity:
90-
// x = sin(y) = cos(pi/2 - y)
91-
// 3: Let z = pi/2 - y, such that x = cos(z)
92-
// 4: From double angle formula; cos(2A) = 1 - sin^2(A):
93-
// z = 2A, z/2 = A
94-
// cos(z) = 1 - 2 * sin^2(z/2)
95-
// 5: Make sin(z/2) subject of the formula:
96-
// sin(z/2) = sqrt((1 - cos(z))/2)
97-
// 6: Recall [3]; x = cos(z). Therefore:
98-
// sin(z/2) = sqrt((1 - x)/2)
99-
// 7: Let u = (1 - x)/2
100-
// 8: Therefore:
101-
// asin(sqrt(u)) = z/2
102-
// 2 * asin(sqrt(u)) = z
103-
// 9: Recall [3], z = pi/2 - y. Therefore:
104-
// y = pi/2 - z
105-
// y = pi/2 - 2 * asin(sqrt(u))
106-
// 10: Recall [1], y = asin(x). Therefore:
107-
// asin(x) = pi/2 - 2 * asin(sqrt(u))
108-
//
109-
// WHY?
110-
// 11: Recall [7], u = (1 - x)/2
111-
// 12: Since 0.5 < x <= 1, therefore:
112-
// 0 <= u <= 0.25 and 0 <= sqrt(u) <= 0.5
113-
//
114-
// Hence, we can reuse the same [0, 0.5] domain polynomial approximation for
115-
// Step [10] as `sqrt(u)` is in range.
116-
117-
// 0x1p-1 < |x| <= 0x1p0, 0.5 < |x| <= 1.0
118-
float xf_abs = (xf < 0 ? -xf : xf);
119-
float sign = (xbits.uintval() >> 15 == 1 ? -1.0 : 1.0);
120-
float u = fputil::multiply_add(-0.5f, xf_abs, 0.5f);
121-
float u_sqrt = fputil::sqrt<float>(u);
122-
123-
// Degree-6 minimax odd polynomial of asin(x) generated by Sollya with:
124-
// > P = fpminimax(asin(x)/x, [|0, 2, 4, 6, 8|], [|SG...|], [0, 0.5]);
125-
float asin_sqrt_u =
126-
u_sqrt * fputil::polyeval(u, 0x1.000002p0f, 0x1.554c2ap-3f,
127-
0x1.3541ccp-4f, 0x1.43b2d6p-5f, 0x1.a0d73ep-5f);
128-
129-
return fputil::cast<float16>(sign *
130-
fputil::multiply_add(-2.0f, asin_sqrt_u, PI_2));
131-
}
14+
LLVM_LIBC_FUNCTION(float16, asinf16, (float16 x)) { return math::asinf16(x); }
13215

13316
} // namespace LIBC_NAMESPACE_DECL

‎libc/test/shared/CMakeLists.txt

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -15,6 +15,7 @@ add_fp_unittest(
1515
libc.src.__support.math.acospif16
1616
libc.src.__support.math.asin
1717
libc.src.__support.math.asinf
18+
libc.src.__support.math.asinf16
1819
libc.src.__support.math.erff
1920
libc.src.__support.math.exp
2021
libc.src.__support.math.exp10

‎libc/test/shared/shared_math_test.cpp

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -17,6 +17,7 @@ TEST(LlvmLibcSharedMathTest, AllFloat16) {
1717

1818
EXPECT_FP_EQ(0x0p+0f16, LIBC_NAMESPACE::shared::acoshf16(1.0f16));
1919
EXPECT_FP_EQ(0x0p+0f16, LIBC_NAMESPACE::shared::acospif16(1.0f16));
20+
EXPECT_FP_EQ(0x0p+0f16, LIBC_NAMESPACE::shared::asinf16(0.0f16));
2021

2122
EXPECT_FP_EQ(0x1p+0f16, LIBC_NAMESPACE::shared::exp10f16(0.0f16));
2223

0 commit comments

Comments
 (0)