Skip to content

Commit f2a476d

Browse files
bassiounixbassiounix
authored
[libc][math] Refactor atanf implementation to header-only in src/__support/math folder. (#150854)
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 3c3523c commit f2a476d

File tree

9 files changed

+188
-123
lines changed

9 files changed

+188
-123
lines changed

‎libc/shared/math.h‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -23,6 +23,7 @@
2323
#include "math/asinhf.h"
2424
#include "math/asinhf16.h"
2525
#include "math/atan.h"
26+
#include "math/atanf.h"
2627
#include "math/erff.h"
2728
#include "math/exp.h"
2829
#include "math/exp10.h"

‎libc/shared/math/atanf.h‎

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

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

Lines changed: 15 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -199,6 +199,21 @@ DEPENDS
199199
libc.src.__support.macros.optimization
200200
)
201201

202+
add_header_library(
203+
atanf
204+
HDRS
205+
atanf.h
206+
DEPENDS
207+
.inv_trigf_utils
208+
libc.src.__support.FPUtil.except_value_utils
209+
libc.src.__support.FPUtil.fp_bits
210+
libc.src.__support.FPUtil.multiply_add
211+
libc.src.__support.FPUtil.nearest_integer
212+
libc.src.__support.FPUtil.polyeval
213+
libc.src.__support.FPUtil.rounding_mode
214+
libc.src.__support.macros.optimization
215+
)
216+
202217
add_header_library(
203218
asinf
204219
HDRS

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

Lines changed: 129 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,129 @@
1+
//===-- Implementation header for atanf -------------------------*- 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_ATANF_H
10+
#define LLVM_LIBC_SRC___SUPPORT_MATH_ATANF_H
11+
12+
#include "inv_trigf_utils.h"
13+
#include "src/__support/FPUtil/FPBits.h"
14+
#include "src/__support/FPUtil/PolyEval.h"
15+
#include "src/__support/FPUtil/except_value_utils.h"
16+
#include "src/__support/FPUtil/multiply_add.h"
17+
#include "src/__support/FPUtil/nearest_integer.h"
18+
#include "src/__support/macros/config.h"
19+
#include "src/__support/macros/optimization.h" // LIBC_UNLIKELY
20+
21+
namespace LIBC_NAMESPACE_DECL {
22+
23+
namespace math {
24+
25+
LIBC_INLINE static constexpr float atanf(float x) {
26+
using namespace inv_trigf_utils_internal;
27+
using FPBits = typename fputil::FPBits<float>;
28+
29+
constexpr double FINAL_SIGN[2] = {1.0, -1.0};
30+
constexpr double SIGNED_PI_OVER_2[2] = {0x1.921fb54442d18p0,
31+
-0x1.921fb54442d18p0};
32+
33+
FPBits x_bits(x);
34+
Sign sign = x_bits.sign();
35+
x_bits.set_sign(Sign::POS);
36+
uint32_t x_abs = x_bits.uintval();
37+
38+
// x is inf or nan, |x| < 2^-4 or |x|= > 16.
39+
if (LIBC_UNLIKELY(x_abs <= 0x3d80'0000U || x_abs >= 0x4180'0000U)) {
40+
double x_d = static_cast<double>(x);
41+
double const_term = 0.0;
42+
if (LIBC_UNLIKELY(x_abs >= 0x4180'0000)) {
43+
// atan(+-Inf) = +-pi/2.
44+
if (x_bits.is_inf()) {
45+
volatile double sign_pi_over_2 = SIGNED_PI_OVER_2[sign.is_neg()];
46+
return static_cast<float>(sign_pi_over_2);
47+
}
48+
if (x_bits.is_nan())
49+
return x;
50+
// x >= 16
51+
x_d = -1.0 / x_d;
52+
const_term = SIGNED_PI_OVER_2[sign.is_neg()];
53+
}
54+
// 0 <= x < 1/16;
55+
if (LIBC_UNLIKELY(x_bits.is_zero()))
56+
return x;
57+
// x <= 2^-12;
58+
if (LIBC_UNLIKELY(x_abs < 0x3980'0000)) {
59+
#if defined(LIBC_TARGET_CPU_HAS_FMA_FLOAT)
60+
return fputil::multiply_add(x, -0x1.0p-25f, x);
61+
#else
62+
double x_d = static_cast<double>(x);
63+
return static_cast<float>(fputil::multiply_add(x_d, -0x1.0p-25, x_d));
64+
#endif // LIBC_TARGET_CPU_HAS_FMA_FLOAT
65+
}
66+
// Use Taylor polynomial:
67+
// atan(x) ~ x * (1 - x^2 / 3 + x^4 / 5 - x^6 / 7 + x^8 / 9 - x^10 / 11).
68+
constexpr double ATAN_TAYLOR[6] = {
69+
0x1.0000000000000p+0, -0x1.5555555555555p-2, 0x1.999999999999ap-3,
70+
-0x1.2492492492492p-3, 0x1.c71c71c71c71cp-4, -0x1.745d1745d1746p-4,
71+
};
72+
double x2 = x_d * x_d;
73+
double x4 = x2 * x2;
74+
double c0 = fputil::multiply_add(x2, ATAN_TAYLOR[1], ATAN_TAYLOR[0]);
75+
double c1 = fputil::multiply_add(x2, ATAN_TAYLOR[3], ATAN_TAYLOR[2]);
76+
double c2 = fputil::multiply_add(x2, ATAN_TAYLOR[5], ATAN_TAYLOR[4]);
77+
double p = fputil::polyeval(x4, c0, c1, c2);
78+
double r = fputil::multiply_add(x_d, p, const_term);
79+
return static_cast<float>(r);
80+
}
81+
82+
// Range reduction steps:
83+
// 1) atan(x) = sign(x) * atan(|x|)
84+
// 2) If |x| > 1, atan(|x|) = pi/2 - atan(1/|x|)
85+
// 3) For 1/16 < x <= 1, we find k such that: |x - k/16| <= 1/32.
86+
// 4) Then we use polynomial approximation:
87+
// atan(x) ~ atan((k/16) + (x - (k/16)) * Q(x - k/16)
88+
// = P(x - k/16)
89+
double x_d = 0, const_term = 0, final_sign = 0;
90+
int idx = 0;
91+
92+
if (x_abs > 0x3f80'0000U) {
93+
// |x| > 1, we need to invert x, so we will perform range reduction in
94+
// double precision.
95+
x_d = 1.0 / static_cast<double>(x_bits.get_val());
96+
double k_d = fputil::nearest_integer(x_d * 0x1.0p4);
97+
x_d = fputil::multiply_add(k_d, -0x1.0p-4, x_d);
98+
idx = static_cast<int>(k_d);
99+
final_sign = FINAL_SIGN[sign.is_pos()];
100+
// Adjust constant term of the polynomial by +- pi/2.
101+
const_term = fputil::multiply_add(final_sign, ATAN_COEFFS[idx][0],
102+
SIGNED_PI_OVER_2[sign.is_neg()]);
103+
} else {
104+
// Exceptional value:
105+
if (LIBC_UNLIKELY(x_abs == 0x3d8d'6b23U)) { // |x| = 0x1.1ad646p-4
106+
return sign.is_pos() ? fputil::round_result_slightly_down(0x1.1a6386p-4f)
107+
: fputil::round_result_slightly_up(-0x1.1a6386p-4f);
108+
}
109+
// Perform range reduction in single precision.
110+
float x_f = x_bits.get_val();
111+
float k_f = fputil::nearest_integer(x_f * 0x1.0p4f);
112+
x_f = fputil::multiply_add(k_f, -0x1.0p-4f, x_f);
113+
x_d = static_cast<double>(x_f);
114+
idx = static_cast<int>(k_f);
115+
final_sign = FINAL_SIGN[sign.is_neg()];
116+
const_term = final_sign * ATAN_COEFFS[idx][0];
117+
}
118+
119+
double p = atan_eval(x_d, idx);
120+
double r = fputil::multiply_add(final_sign * x_d, p, const_term);
121+
122+
return static_cast<float>(r);
123+
}
124+
125+
} // namespace math
126+
127+
} // namespace LIBC_NAMESPACE_DECL
128+
129+
#endif // LLVM_LIBC_SRC___SUPPORT_MATH_ATANF_H

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

Lines changed: 1 addition & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -4027,14 +4027,7 @@ add_entrypoint_object(
40274027
HDRS
40284028
../atanf.h
40294029
DEPENDS
4030-
libc.src.__support.FPUtil.except_value_utils
4031-
libc.src.__support.FPUtil.fp_bits
4032-
libc.src.__support.FPUtil.multiply_add
4033-
libc.src.__support.FPUtil.nearest_integer
4034-
libc.src.__support.FPUtil.polyeval
4035-
libc.src.__support.FPUtil.rounding_mode
4036-
libc.src.__support.macros.optimization
4037-
libc.src.__support.math.inv_trigf_utils
4030+
libc.src.__support.math.atanf
40384031
)
40394032

40404033
add_entrypoint_object(

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

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

99
#include "src/math/atanf.h"
10-
#include "src/__support/FPUtil/FPBits.h"
11-
#include "src/__support/FPUtil/PolyEval.h"
12-
#include "src/__support/FPUtil/except_value_utils.h"
13-
#include "src/__support/FPUtil/multiply_add.h"
14-
#include "src/__support/FPUtil/nearest_integer.h"
15-
#include "src/__support/FPUtil/rounding_mode.h"
16-
#include "src/__support/macros/config.h"
17-
#include "src/__support/macros/optimization.h" // LIBC_UNLIKELY
18-
#include "src/__support/math/inv_trigf_utils.h"
10+
#include "src/__support/math/atanf.h"
1911

2012
namespace LIBC_NAMESPACE_DECL {
2113

22-
LLVM_LIBC_FUNCTION(float, atanf, (float x)) {
23-
using namespace inv_trigf_utils_internal;
24-
using FPBits = typename fputil::FPBits<float>;
25-
26-
constexpr double FINAL_SIGN[2] = {1.0, -1.0};
27-
constexpr double SIGNED_PI_OVER_2[2] = {0x1.921fb54442d18p0,
28-
-0x1.921fb54442d18p0};
29-
30-
FPBits x_bits(x);
31-
Sign sign = x_bits.sign();
32-
x_bits.set_sign(Sign::POS);
33-
uint32_t x_abs = x_bits.uintval();
34-
35-
// x is inf or nan, |x| < 2^-4 or |x|= > 16.
36-
if (LIBC_UNLIKELY(x_abs <= 0x3d80'0000U || x_abs >= 0x4180'0000U)) {
37-
double x_d = static_cast<double>(x);
38-
double const_term = 0.0;
39-
if (LIBC_UNLIKELY(x_abs >= 0x4180'0000)) {
40-
// atan(+-Inf) = +-pi/2.
41-
if (x_bits.is_inf()) {
42-
volatile double sign_pi_over_2 = SIGNED_PI_OVER_2[sign.is_neg()];
43-
return static_cast<float>(sign_pi_over_2);
44-
}
45-
if (x_bits.is_nan())
46-
return x;
47-
// x >= 16
48-
x_d = -1.0 / x_d;
49-
const_term = SIGNED_PI_OVER_2[sign.is_neg()];
50-
}
51-
// 0 <= x < 1/16;
52-
if (LIBC_UNLIKELY(x_bits.is_zero()))
53-
return x;
54-
// x <= 2^-12;
55-
if (LIBC_UNLIKELY(x_abs < 0x3980'0000)) {
56-
#if defined(LIBC_TARGET_CPU_HAS_FMA_FLOAT)
57-
return fputil::multiply_add(x, -0x1.0p-25f, x);
58-
#else
59-
double x_d = static_cast<double>(x);
60-
return static_cast<float>(fputil::multiply_add(x_d, -0x1.0p-25, x_d));
61-
#endif // LIBC_TARGET_CPU_HAS_FMA_FLOAT
62-
}
63-
// Use Taylor polynomial:
64-
// atan(x) ~ x * (1 - x^2 / 3 + x^4 / 5 - x^6 / 7 + x^8 / 9 - x^10 / 11).
65-
constexpr double ATAN_TAYLOR[6] = {
66-
0x1.0000000000000p+0, -0x1.5555555555555p-2, 0x1.999999999999ap-3,
67-
-0x1.2492492492492p-3, 0x1.c71c71c71c71cp-4, -0x1.745d1745d1746p-4,
68-
};
69-
double x2 = x_d * x_d;
70-
double x4 = x2 * x2;
71-
double c0 = fputil::multiply_add(x2, ATAN_TAYLOR[1], ATAN_TAYLOR[0]);
72-
double c1 = fputil::multiply_add(x2, ATAN_TAYLOR[3], ATAN_TAYLOR[2]);
73-
double c2 = fputil::multiply_add(x2, ATAN_TAYLOR[5], ATAN_TAYLOR[4]);
74-
double p = fputil::polyeval(x4, c0, c1, c2);
75-
double r = fputil::multiply_add(x_d, p, const_term);
76-
return static_cast<float>(r);
77-
}
78-
79-
// Range reduction steps:
80-
// 1) atan(x) = sign(x) * atan(|x|)
81-
// 2) If |x| > 1, atan(|x|) = pi/2 - atan(1/|x|)
82-
// 3) For 1/16 < x <= 1, we find k such that: |x - k/16| <= 1/32.
83-
// 4) Then we use polynomial approximation:
84-
// atan(x) ~ atan((k/16) + (x - (k/16)) * Q(x - k/16)
85-
// = P(x - k/16)
86-
double x_d, const_term, final_sign;
87-
int idx;
88-
89-
if (x_abs > 0x3f80'0000U) {
90-
// |x| > 1, we need to invert x, so we will perform range reduction in
91-
// double precision.
92-
x_d = 1.0 / static_cast<double>(x_bits.get_val());
93-
double k_d = fputil::nearest_integer(x_d * 0x1.0p4);
94-
x_d = fputil::multiply_add(k_d, -0x1.0p-4, x_d);
95-
idx = static_cast<int>(k_d);
96-
final_sign = FINAL_SIGN[sign.is_pos()];
97-
// Adjust constant term of the polynomial by +- pi/2.
98-
const_term = fputil::multiply_add(final_sign, ATAN_COEFFS[idx][0],
99-
SIGNED_PI_OVER_2[sign.is_neg()]);
100-
} else {
101-
// Exceptional value:
102-
if (LIBC_UNLIKELY(x_abs == 0x3d8d'6b23U)) { // |x| = 0x1.1ad646p-4
103-
return sign.is_pos() ? fputil::round_result_slightly_down(0x1.1a6386p-4f)
104-
: fputil::round_result_slightly_up(-0x1.1a6386p-4f);
105-
}
106-
// Perform range reduction in single precision.
107-
float x_f = x_bits.get_val();
108-
float k_f = fputil::nearest_integer(x_f * 0x1.0p4f);
109-
x_f = fputil::multiply_add(k_f, -0x1.0p-4f, x_f);
110-
x_d = static_cast<double>(x_f);
111-
idx = static_cast<int>(k_f);
112-
final_sign = FINAL_SIGN[sign.is_neg()];
113-
const_term = final_sign * ATAN_COEFFS[idx][0];
114-
}
115-
116-
double p = atan_eval(x_d, idx);
117-
double r = fputil::multiply_add(final_sign * x_d, p, const_term);
118-
119-
return static_cast<float>(r);
120-
}
14+
LLVM_LIBC_FUNCTION(float, atanf, (float x)) { return math::atanf(x); }
12115

12216
} // namespace LIBC_NAMESPACE_DECL

‎libc/test/shared/CMakeLists.txt‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -19,6 +19,7 @@ add_fp_unittest(
1919
libc.src.__support.math.asinhf
2020
libc.src.__support.math.asinhf16
2121
libc.src.__support.math.atan
22+
libc.src.__support.math.atanf
2223
libc.src.__support.math.erff
2324
libc.src.__support.math.exp
2425
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
@@ -44,6 +44,7 @@ TEST(LlvmLibcSharedMathTest, AllFloat) {
4444
EXPECT_FP_EQ(0x0p+0f, LIBC_NAMESPACE::shared::acoshf(1.0f));
4545
EXPECT_FP_EQ(0x0p+0f, LIBC_NAMESPACE::shared::asinf(0.0f));
4646
EXPECT_FP_EQ(0x0p+0f, LIBC_NAMESPACE::shared::asinhf(0.0f));
47+
EXPECT_FP_EQ(0x0p+0f, LIBC_NAMESPACE::shared::atanf(0.0f));
4748
EXPECT_FP_EQ(0x0p+0f, LIBC_NAMESPACE::shared::erff(0.0f));
4849
EXPECT_FP_EQ(0x1p+0f, LIBC_NAMESPACE::shared::exp10f(0.0f));
4950
EXPECT_FP_EQ(0x1p+0f, LIBC_NAMESPACE::shared::expf(0.0f));

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

Lines changed: 15 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -2274,6 +2274,20 @@ libc_support_library(
22742274
],
22752275
)
22762276

2277+
libc_support_library(
2278+
name = "__support_math_atanf",
2279+
hdrs = ["src/__support/math/atanf.h"],
2280+
deps = [
2281+
":__support_fputil_except_value_utils",
2282+
":__support_fputil_multiply_add",
2283+
":__support_fputil_nearest_integer",
2284+
":__support_fputil_polyeval",
2285+
":__support_fputil_rounding_mode",
2286+
":__support_macros_optimization",
2287+
":__support_math_inv_trigf_utils",
2288+
],
2289+
)
2290+
22772291
libc_support_library(
22782292
name = "__support_math_asinf",
22792293
hdrs = ["src/__support/math/asinf.h"],
@@ -2886,13 +2900,7 @@ libc_math_function(
28862900
libc_math_function(
28872901
name = "atanf",
28882902
additional_deps = [
2889-
":__support_fputil_fma",
2890-
":__support_fputil_multiply_add",
2891-
":__support_fputil_nearest_integer",
2892-
":__support_fputil_polyeval",
2893-
":__support_fputil_rounding_mode",
2894-
":__support_macros_optimization",
2895-
":__support_math_inv_trigf_utils",
2903+
":__support_math_atanf"
28962904
],
28972905
)
28982906

0 commit comments

Comments
 (0)