Skip to content

Commit 205f58b

Browse files
bassiounixbassiounix
committed
[libc][math] Refactor acospif16 implementation to header-only in src/__support/math folder.
1 parent 558975c commit 205f58b

File tree

7 files changed

+220
-129
lines changed

7 files changed

+220
-129
lines changed

‎libc/shared/math.h

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -16,6 +16,7 @@
1616
#include "math/acosf16.h"
1717
#include "math/acoshf.h"
1818
#include "math/acoshf16.h"
19+
#include "math/acospif16.h"
1920
#include "math/erff.h"
2021
#include "math/exp.h"
2122
#include "math/exp10.h"

‎libc/shared/math/acospif16.h

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

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

Lines changed: 15 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -95,6 +95,21 @@ add_header_library(
9595
libc.src.__support.macros.optimization
9696
)
9797

98+
add_header_library(
99+
acospif16
100+
HDRS
101+
acospif16.h
102+
DEPENDS
103+
libc.src.__support.FPUtil.cast
104+
libc.src.__support.FPUtil.fenv_impl
105+
libc.src.__support.FPUtil.fp_bits
106+
libc.src.__support.FPUtil.multiply_add
107+
libc.src.__support.FPUtil.polyeval
108+
libc.src.__support.FPUtil.sqrt
109+
libc.src.__support.macros.optimization
110+
libc.src.__support.macros.properties.types
111+
)
112+
98113
add_header_library(
99114
asin_utils
100115
HDRS

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

Lines changed: 147 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,147 @@
1+
//===-- Implementation header for acospif16 ---------------------*- 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_ACOSPIF16_H
10+
#define LLVM_LIBC_SRC___SUPPORT_MATH_ACOSPIF16_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+
static constexpr float16 acospif16(float16 x) {
29+
using FPBits = fputil::FPBits<float16>;
30+
FPBits xbits(x);
31+
32+
uint16_t x_u = xbits.uintval();
33+
uint16_t x_abs = x_u & 0x7fff;
34+
uint16_t x_sign = x_u >> 15;
35+
36+
// |x| > 0x1p0, |x| > 1, or x is NaN.
37+
if (LIBC_UNLIKELY(x_abs > 0x3c00)) {
38+
// acospif16(NaN) = NaN
39+
if (xbits.is_nan()) {
40+
if (xbits.is_signaling_nan()) {
41+
fputil::raise_except_if_required(FE_INVALID);
42+
return FPBits::quiet_nan().get_val();
43+
}
44+
45+
return x;
46+
}
47+
48+
// 1 < |x| <= +inf
49+
fputil::raise_except_if_required(FE_INVALID);
50+
fputil::set_errno_if_required(EDOM);
51+
52+
return FPBits::quiet_nan().get_val();
53+
}
54+
55+
// |x| == 0x1p0, x is 1 or -1
56+
// if x is (-)1, return 1
57+
// if x is (+)1, return 0
58+
if (LIBC_UNLIKELY(x_abs == 0x3c00))
59+
return fputil::cast<float16>(x_sign ? 1.0f : 0.0f);
60+
61+
float xf = x;
62+
float xsq = xf * xf;
63+
64+
// Degree-6 minimax polynomial coefficients of asin(x) generated by Sollya
65+
// with: > P = fpminimax(asin(x)/(pi * x), [|0, 2, 4, 6, 8|], [|SG...|], [0,
66+
// 0.5]);
67+
constexpr float POLY_COEFFS[5] = {0x1.45f308p-2f, 0x1.b2900cp-5f,
68+
0x1.897e36p-6f, 0x1.9efafcp-7f,
69+
0x1.06d884p-6f};
70+
// |x| <= 0x1p-1, |x| <= 0.5
71+
if (x_abs <= 0x3800) {
72+
// if x is 0, return 0.5
73+
if (LIBC_UNLIKELY(x_abs == 0))
74+
return fputil::cast<float16>(0.5f);
75+
76+
// Note that: acos(x) = pi/2 + asin(-x) = pi/2 - asin(x), then
77+
// acospi(x) = 0.5 - asin(x)/pi
78+
float interm =
79+
fputil::polyeval(xsq, POLY_COEFFS[0], POLY_COEFFS[1], POLY_COEFFS[2],
80+
POLY_COEFFS[3], POLY_COEFFS[4]);
81+
82+
return fputil::cast<float16>(fputil::multiply_add(-xf, interm, 0.5f));
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 - 2 * 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+
// 11: Recall that: acos(x) = pi/2 + asin(-x) = pi/2 - asin(x)
109+
// Therefore:
110+
// acos(x) = pi/2 - (pi/2 - 2 * asin(sqrt(u)))
111+
// acos(x) = 2 * asin(sqrt(u))
112+
// acospi(x) = 2 * (asin(sqrt(u)) / pi)
113+
//
114+
// THE RANGE REDUCTION, HOW?
115+
// 12: Recall [7], u = (1 - x)/2
116+
// 13: Since 0.5 < x <= 1, therefore:
117+
// 0 <= u <= 0.25 and 0 <= sqrt(u) <= 0.5
118+
//
119+
// Hence, we can reuse the same [0, 0.5] domain polynomial approximation for
120+
// Step [11] as `sqrt(u)` is in range.
121+
// When -1 < x <= -0.5, the identity:
122+
// acos(x) = pi - acos(-x)
123+
// acospi(x) = 1 - acos(-x)/pi
124+
// allows us to compute for the negative x value (lhs)
125+
// with a positive x value instead (rhs).
126+
127+
float xf_abs = (xf < 0 ? -xf : xf);
128+
float u = fputil::multiply_add(-0.5f, xf_abs, 0.5f);
129+
float sqrt_u = fputil::sqrt<float>(u);
130+
131+
float asin_sqrt_u =
132+
sqrt_u * fputil::polyeval(u, POLY_COEFFS[0], POLY_COEFFS[1],
133+
POLY_COEFFS[2], POLY_COEFFS[3], POLY_COEFFS[4]);
134+
135+
// Same as acos(x), but devided the expression with pi
136+
return fputil::cast<float16>(
137+
x_sign ? fputil::multiply_add(-2.0f, asin_sqrt_u, 1.0f)
138+
: 2.0f * asin_sqrt_u);
139+
}
140+
141+
} // namespace math
142+
143+
} // namespace LIBC_NAMESPACE_DECL
144+
145+
#endif // LIBC_TYPES_HAS_FLOAT16
146+
147+
#endif // LLVM_LIBC_SRC___SUPPORT_MATH_ACOSPIF16_H

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

Lines changed: 2 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -4045,16 +4045,8 @@ add_entrypoint_object(
40454045
HDRS
40464046
../acospif16.h
40474047
DEPENDS
4048-
libc.hdr.errno_macros
4049-
libc.hdr.fenv_macros
4050-
libc.src.__support.FPUtil.cast
4051-
libc.src.__support.FPUtil.fenv_impl
4052-
libc.src.__support.FPUtil.fp_bits
4053-
libc.src.__support.FPUtil.multiply_add
4054-
libc.src.__support.FPUtil.polyeval
4055-
libc.src.__support.FPUtil.sqrt
4056-
libc.src.__support.macros.optimization
4057-
libc.src.__support.macros.properties.types
4048+
libc.src.__support.math.acospif16
4049+
libc.src.errno.errno
40584050
)
40594051

40604052
add_header_library(

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

Lines changed: 3 additions & 119 deletions
Original file line numberDiff line numberDiff line change
@@ -7,128 +7,12 @@
77
//===----------------------------------------------------------------------===//
88

99
#include "src/math/acospif16.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/acospif16.h"
1911

2012
namespace LIBC_NAMESPACE_DECL {
2113

2214
LLVM_LIBC_FUNCTION(float16, acospif16, (float16 x)) {
23-
using FPBits = fputil::FPBits<float16>;
24-
FPBits xbits(x);
25-
26-
uint16_t x_u = xbits.uintval();
27-
uint16_t x_abs = x_u & 0x7fff;
28-
uint16_t x_sign = x_u >> 15;
29-
30-
// |x| > 0x1p0, |x| > 1, or x is NaN.
31-
if (LIBC_UNLIKELY(x_abs > 0x3c00)) {
32-
// acospif16(NaN) = NaN
33-
if (xbits.is_nan()) {
34-
if (xbits.is_signaling_nan()) {
35-
fputil::raise_except_if_required(FE_INVALID);
36-
return FPBits::quiet_nan().get_val();
37-
}
38-
39-
return x;
40-
}
41-
42-
// 1 < |x| <= +inf
43-
fputil::raise_except_if_required(FE_INVALID);
44-
fputil::set_errno_if_required(EDOM);
45-
46-
return FPBits::quiet_nan().get_val();
47-
}
48-
49-
// |x| == 0x1p0, x is 1 or -1
50-
// if x is (-)1, return 1
51-
// if x is (+)1, return 0
52-
if (LIBC_UNLIKELY(x_abs == 0x3c00))
53-
return fputil::cast<float16>(x_sign ? 1.0f : 0.0f);
54-
55-
float xf = x;
56-
float xsq = xf * xf;
57-
58-
// Degree-6 minimax polynomial coefficients of asin(x) generated by Sollya
59-
// with: > P = fpminimax(asin(x)/(pi * x), [|0, 2, 4, 6, 8|], [|SG...|], [0,
60-
// 0.5]);
61-
constexpr float POLY_COEFFS[5] = {0x1.45f308p-2f, 0x1.b2900cp-5f,
62-
0x1.897e36p-6f, 0x1.9efafcp-7f,
63-
0x1.06d884p-6f};
64-
// |x| <= 0x1p-1, |x| <= 0.5
65-
if (x_abs <= 0x3800) {
66-
// if x is 0, return 0.5
67-
if (LIBC_UNLIKELY(x_abs == 0))
68-
return fputil::cast<float16>(0.5f);
69-
70-
// Note that: acos(x) = pi/2 + asin(-x) = pi/2 - asin(x), then
71-
// acospi(x) = 0.5 - asin(x)/pi
72-
float interm =
73-
fputil::polyeval(xsq, POLY_COEFFS[0], POLY_COEFFS[1], POLY_COEFFS[2],
74-
POLY_COEFFS[3], POLY_COEFFS[4]);
75-
76-
return fputil::cast<float16>(fputil::multiply_add(-xf, interm, 0.5f));
77-
}
78-
79-
// When |x| > 0.5, assume that 0.5 < |x| <= 1
80-
//
81-
// Step-by-step range-reduction proof:
82-
// 1: Let y = asin(x), such that, x = sin(y)
83-
// 2: From complimentary angle identity:
84-
// x = sin(y) = cos(pi/2 - y)
85-
// 3: Let z = pi/2 - y, such that x = cos(z)
86-
// 4: From double angle formula; cos(2A) = 1 - 2 * sin^2(A):
87-
// z = 2A, z/2 = A
88-
// cos(z) = 1 - 2 * sin^2(z/2)
89-
// 5: Make sin(z/2) subject of the formula:
90-
// sin(z/2) = sqrt((1 - cos(z))/2)
91-
// 6: Recall [3]; x = cos(z). Therefore:
92-
// sin(z/2) = sqrt((1 - x)/2)
93-
// 7: Let u = (1 - x)/2
94-
// 8: Therefore:
95-
// asin(sqrt(u)) = z/2
96-
// 2 * asin(sqrt(u)) = z
97-
// 9: Recall [3]; z = pi/2 - y. Therefore:
98-
// y = pi/2 - z
99-
// y = pi/2 - 2 * asin(sqrt(u))
100-
// 10: Recall [1], y = asin(x). Therefore:
101-
// asin(x) = pi/2 - 2 * asin(sqrt(u))
102-
// 11: Recall that: acos(x) = pi/2 + asin(-x) = pi/2 - asin(x)
103-
// Therefore:
104-
// acos(x) = pi/2 - (pi/2 - 2 * asin(sqrt(u)))
105-
// acos(x) = 2 * asin(sqrt(u))
106-
// acospi(x) = 2 * (asin(sqrt(u)) / pi)
107-
//
108-
// THE RANGE REDUCTION, HOW?
109-
// 12: Recall [7], u = (1 - x)/2
110-
// 13: Since 0.5 < x <= 1, therefore:
111-
// 0 <= u <= 0.25 and 0 <= sqrt(u) <= 0.5
112-
//
113-
// Hence, we can reuse the same [0, 0.5] domain polynomial approximation for
114-
// Step [11] as `sqrt(u)` is in range.
115-
// When -1 < x <= -0.5, the identity:
116-
// acos(x) = pi - acos(-x)
117-
// acospi(x) = 1 - acos(-x)/pi
118-
// allows us to compute for the negative x value (lhs)
119-
// with a positive x value instead (rhs).
120-
121-
float xf_abs = (xf < 0 ? -xf : xf);
122-
float u = fputil::multiply_add(-0.5f, xf_abs, 0.5f);
123-
float sqrt_u = fputil::sqrt<float>(u);
124-
125-
float asin_sqrt_u =
126-
sqrt_u * fputil::polyeval(u, POLY_COEFFS[0], POLY_COEFFS[1],
127-
POLY_COEFFS[2], POLY_COEFFS[3], POLY_COEFFS[4]);
128-
129-
// Same as acos(x), but devided the expression with pi
130-
return fputil::cast<float16>(
131-
x_sign ? fputil::multiply_add(-2.0f, asin_sqrt_u, 1.0f)
132-
: 2.0f * asin_sqrt_u);
15+
return math::acospif16(x);
13316
}
17+
13418
} // namespace LIBC_NAMESPACE_DECL

‎utils/bazel/llvm-project-overlay/libc/BUILD.bazel

Lines changed: 23 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -2158,6 +2158,21 @@ libc_support_library(
21582158
],
21592159
)
21602160

2161+
libc_support_library(
2162+
name = "__support_math_acospif16",
2163+
hdrs = ["src/__support/math/acospif16.h"],
2164+
deps = [
2165+
":__support_fputil_cast",
2166+
":__support_fputil_fenv_impl",
2167+
":__support_fputil_fp_bits",
2168+
":__support_fputil_multiply_add",
2169+
":__support_fputil_polyeval",
2170+
":__support_fputil_sqrt",
2171+
":__support_macros_optimization",
2172+
":__support_macros_properties_types",
2173+
],
2174+
)
2175+
21612176
libc_support_library(
21622177
name = "__support_math_asin_utils",
21632178
hdrs = ["src/__support/math/asin_utils.h"],
@@ -2707,6 +2722,14 @@ libc_math_function(
27072722
],
27082723
)
27092724

2725+
libc_math_function(
2726+
name = "acospif16",
2727+
additional_deps = [
2728+
":__support_math_acospif16",
2729+
":errno",
2730+
],
2731+
)
2732+
27102733
libc_math_function(
27112734
name = "asinf",
27122735
additional_deps = [

0 commit comments

Comments
 (0)