[libc][math][c23] Implement C23 math function atanpif16 #150400
+355
−1
Conversation
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
|
@llvm/pr-subscribers-libc Author: Mohamed Emad (hulxv) ChangesThis PR implements Mathematical ImplementationThe implementation uses a 15th-degree Taylor polynomial expansion of $$ To ensure accuracy across all real inputs, the domain is divided into three cases with appropriate transformations: Case 1: Case 2: Case 3: Closes #132212 Full diff: https://github.com/llvm/llvm-project/pull/150400.diff 14 Files Affected:
diff --git a/libc/config/linux/x86_64/entrypoints.txt b/libc/config/linux/x86_64/entrypoints.txt
index 73dfeae1a2c94..1585f803229b6 100644
--- a/libc/config/linux/x86_64/entrypoints.txt
+++ b/libc/config/linux/x86_64/entrypoints.txt
@@ -664,6 +664,7 @@ if(LIBC_TYPES_HAS_FLOAT16)
libc.src.math.acoshf16
libc.src.math.asinf16
libc.src.math.asinhf16
+ libc.src.math.atanpif16
libc.src.math.canonicalizef16
libc.src.math.ceilf16
libc.src.math.copysignf16
diff --git a/libc/include/math.yaml b/libc/include/math.yaml
index fef829422244d..06a67fdb44a9c 100644
--- a/libc/include/math.yaml
+++ b/libc/include/math.yaml
@@ -113,6 +113,12 @@ functions:
return_type: float
arguments:
- type: float
+ - name: atanpif16
+ standards:
+ - stdc
+ return_type: _Float16
+ arguments:
+ - type: _Float16
- name: canonicalize
standards:
- stdc
diff --git a/libc/src/math/CMakeLists.txt b/libc/src/math/CMakeLists.txt
index d177ff79141c0..033aa2396b86b 100644
--- a/libc/src/math/CMakeLists.txt
+++ b/libc/src/math/CMakeLists.txt
@@ -67,6 +67,8 @@ add_math_entrypoint_object(atan2f128)
add_math_entrypoint_object(atanh)
add_math_entrypoint_object(atanhf)
+add_math_entrypoint_object(atanpif16)
+
add_math_entrypoint_object(canonicalize)
add_math_entrypoint_object(canonicalizef)
add_math_entrypoint_object(canonicalizel)
diff --git a/libc/src/math/atanpif16.h b/libc/src/math/atanpif16.h
new file mode 100644
index 0000000000000..14dc150014578
--- /dev/null
+++ b/libc/src/math/atanpif16.h
@@ -0,0 +1,22 @@
+//===-- Implementation header for atanpif16 ---------------------*- C++ -*-===//
+//
+// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
+// See https://llvm.org/LICENSE.txt for license information.
+// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception.
+//
+//===----------------------------------------------------------------------===//
+
+#ifndef LLVM_LIBC_SRC_MATH_ATANPIF16_H
+#define LLVM_LIBC_SRC_MATH_ATANPIF16_H
+
+#include "src/__support/macros/config.h"
+#include "src/__support/macros/properties/types.h"
+
+namespace LIBC_NAMESPACE_DECL {
+
+float16 atanpif16(float16 x);
+
+} // namespace LIBC_NAMESPACE_DECL
+
+#endif // LLVM_LIBC_SRC_MATH_ASINF16_H
+
diff --git a/libc/src/math/generic/CMakeLists.txt b/libc/src/math/generic/CMakeLists.txt
index adbed5b2de48c..421d8163ea792 100644
--- a/libc/src/math/generic/CMakeLists.txt
+++ b/libc/src/math/generic/CMakeLists.txt
@@ -4029,6 +4029,23 @@ add_entrypoint_object(
libc.src.__support.macros.optimization
)
+add_entrypoint_object(
+ atanpif16
+ SRCS
+ atanpif16.cpp
+ HDRS
+ ../atanpif16.h
+ DEPENDS
+ libc.hdr.errno_macros
+ libc.hdr.fenv_macros
+ libc.src.__support.FPUtil.fenv_impl
+ libc.src.__support.FPUtil.fp_bits
+ libc.src.__support.FPUtil.polyeval
+ libc.src.__support.FPUtil.cast
+ libc.src.__support.FPUtil.multiply_add
+ libc.src.__support.FPUtil.sqrt
+)
+
add_object_library(
inv_trigf_utils
HDRS
diff --git a/libc/src/math/generic/atanpif16.cpp b/libc/src/math/generic/atanpif16.cpp
new file mode 100644
index 0000000000000..985cf3738614e
--- /dev/null
+++ b/libc/src/math/generic/atanpif16.cpp
@@ -0,0 +1,155 @@
+//===-- Half-precision atanpi function ------------------------------------===//
+//
+// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
+// See https://llvm.org/LICENSE.txt for license information.
+// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
+//
+//===----------------------------------------------------------------------===//
+
+#include "src/math/atanpif16.h"
+#include "hdr/errno_macros.h"
+#include "hdr/fenv_macros.h"
+#include "src/__support/FPUtil/FEnvImpl.h"
+#include "src/__support/FPUtil/FPBits.h"
+#include "src/__support/FPUtil/PolyEval.h"
+#include "src/__support/FPUtil/cast.h"
+#include "src/__support/FPUtil/multiply_add.h"
+#include "src/__support/FPUtil/sqrt.h"
+
+namespace LIBC_NAMESPACE_DECL {
+
+// Using Python's SymPy library, we can obtain the polynomial approximation of
+// arctan(x)/pi. The steps are as follows:
+// >>> from sympy import *
+// >>> import math
+// >>> x = symbols('x')
+// >>> print(series(atan(x)/math.pi, x, 0, 17))
+//
+// Output:
+// 0.318309886183791*x - 0.106103295394597*x**3 + 0.0636619772367581*x**5 -
+// 0.0454728408833987*x**7 + 0.0353677651315323*x**9 - 0.0289372623803446*x**11
+// + 0.0244853758602916*x**13 - 0.0212206590789194*x**15 + O(x**17)
+//
+// We will assign this 19-degree Taylor polynomial as g(x). This polynomial
+// approximation is accurate for arctan(x)/pi when |x| is in the range [0, 0.5].
+//
+//
+// To compute arctan(x) for all real x, we divide the domain into the following
+// cases:
+//
+// * Case 1: |x| <= 0.5
+// In this range, the direct polynomial approximation is used:
+// arctan(x)/pi = sign(x) * g(|x|)
+// or equivalently, arctan(x) = sign(x) * pi * g(|x|).
+//
+// * Case 2: 0.5 < |x| <= 1
+// We use the double-angle identity for the tangent function, specifically:
+// arctan(x) = 2 * arctan(x / (1 + sqrt(1 + x^2))).
+// Applying this, we have:
+// arctan(x)/pi = sign(x) * 2 * arctan(x')/pi,
+// where x' = |x| / (1 + sqrt(1 + x^2)).
+// Thus, arctan(x)/pi = sign(x) * 2 * g(x')
+//
+// When |x| is in (0.5, 1], the value of x' will always fall within the
+// interval [0.207, 0.414], which is within the accurate range of g(x).
+//
+// * Case 3: |x| > 1
+// For values of |x| greater than 1, we use the reciprocal transformation
+// identity:
+// arctan(x) = pi/2 - arctan(1/x) for x > 0.
+// For any x (real number), this generalizes to:
+// arctan(x)/pi = sign(x) * (1/2 - arctan(1/|x|)/pi).
+// Then, using g(x) for arctan(1/|x|)/pi:
+// arctan(x)/pi = sign(x) * (1/2 - g(1/|x|)).
+//
+// Note that if 1/|x| still falls outside the
+// g(x)'s primary range of accuracy (i.e., if 0.5 < 1/|x| <= 1), the rule
+// from Case 2 must be applied recursively to 1/|x|.
+
+LLVM_LIBC_FUNCTION(float16, atanpif16, (float16 x)) {
+ using FPBits = fputil::FPBits<float16>;
+
+ FPBits xbits(x);
+ bool is_neg = xbits.is_neg();
+
+ auto signed_result = [is_neg](double r) -> float16 {
+ return fputil::cast<float16>(is_neg ? -r : r);
+ };
+
+ if (LIBC_UNLIKELY(xbits.is_inf_or_nan())) {
+ if (xbits.is_nan()) {
+ return x;
+ }
+ // atanpi(±∞) = ±0.5
+ return signed_result(0.5);
+ }
+
+ if (LIBC_UNLIKELY(xbits.is_zero())) {
+ return x;
+ }
+
+ double x_abs = fputil::cast<double>(xbits.abs().get_val());
+
+ if (LIBC_UNLIKELY(x_abs == 1.0)) {
+ return signed_result(0.25);
+ }
+
+ // polynomial coefficients for atan(x)/pi taylor series
+ // generated using sympy: series(atan(x)/pi, x, 0, 17)
+ constexpr double POLY_COEFFS[] = {
+ 0x1.45f306dc9c889p-2, // x^1: 1/pi
+ -0x1.b2995e7b7b60bp-4, // x^3: -1/(3*pi)
+ 0x1.04c26be3b06ccp-4, // x^5: 1/(5*pi)
+ -0x1.7483758e69c08p-5, // x^7: -1/(7*pi)
+ 0x1.21bb945252403p-5, // x^9: 1/(9*pi)
+ -0x1.da1bace3cc68ep-6, // x^11: -1/(11*pi)
+ 0x1.912b1c2336cf2p-6, // x^13: 1/(13*pi)
+ -0x1.5bade52f95e7p-6, // x^15: -1/(15*pi)
+ };
+
+ // evaluate atan(x)/pi using polynomial approximation, valid for |x| <= 0.5
+ constexpr auto atanpi_eval = [](double x) -> double {
+ double xx = x * x;
+ return x * fputil::polyeval(xx, POLY_COEFFS[0], POLY_COEFFS[1],
+ POLY_COEFFS[2], POLY_COEFFS[3], POLY_COEFFS[4],
+ POLY_COEFFS[5], POLY_COEFFS[6], POLY_COEFFS[7]);
+ };
+
+ // Case 1: |x| <= 0.5 - Direct polynomial evaluation
+ if (LIBC_LIKELY(x_abs <= 0.5)) {
+ double result = atanpi_eval(x_abs);
+ return signed_result(result);
+ }
+
+ // case 2: 0.5 < |x| <= 1 - use double-angle reduction
+ // atan(x) = 2 * atan(x / (1 + sqrt(1 + x^2)))
+ // so atanpi(x) = 2 * atanpi(x') where x' = x / (1 + sqrt(1 + x^2))
+ if (x_abs <= 1.0) {
+ double x2 = x_abs * x_abs;
+ double sqrt_term = fputil::sqrt<double>(1.0 + x2);
+ double x_prime = x_abs / (1.0 + sqrt_term);
+ double result = 2.0 * atanpi_eval(x_prime);
+ return signed_result(result);
+ }
+
+ // case 3: |x| > 1 - use reciprocal transformation
+ // atan(x) = pi/2 - atan(1/x) for x > 0
+ // so atanpi(x) = 1/2 - atanpi(1/x)
+ double x_recip = 1.0 / x_abs;
+ double result;
+
+ // if 1/|x| > 0.5, we need to apply Case 2 transformation to 1/|x|
+ if (x_recip > 0.5) {
+ double xx_recip = x_recip * x_recip;
+ double sqrt_term = fputil::sqrt<double>(1.0 + xx_recip);
+ double x_prime = x_recip / (1.0 + sqrt_term);
+ result = fputil::multiply_add(-2.0, atanpi_eval(x_prime), 0.5);
+ } else {
+ // direct evaluation since 1/|x| <= 0.5
+ result = 0.5 - atanpi_eval(x_recip);
+ }
+
+ return signed_result(result);
+}
+
+} // namespace LIBC_NAMESPACE_DECL
diff --git a/libc/test/src/math/CMakeLists.txt b/libc/test/src/math/CMakeLists.txt
index 7ee8b86135557..b17c7bc50fa42 100644
--- a/libc/test/src/math/CMakeLists.txt
+++ b/libc/test/src/math/CMakeLists.txt
@@ -2143,6 +2143,17 @@ add_fp_unittest(
libc.src.__support.FPUtil.fp_bits
)
+add_fp_unittest(
+ atanpif16_test
+ NEED_MPFR
+ SUITE
+ libc-math-unittests
+ SRCS
+ atanpif16_test.cpp
+ DEPENDS
+ libc.src.math.atanpif16
+)
+
add_fp_unittest(
fmul_test
NEED_MPFR
diff --git a/libc/test/src/math/atanpif16_test.cpp b/libc/test/src/math/atanpif16_test.cpp
new file mode 100644
index 0000000000000..2c24b5048a2fb
--- /dev/null
+++ b/libc/test/src/math/atanpif16_test.cpp
@@ -0,0 +1,40 @@
+//===-- Exhaustive test for atanpif16 -------------------------------------===//
+//
+// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
+// See https://llvm.org/LICENSE.txt for license information.
+// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
+//
+//===----------------------------------------------------------------------===//
+
+#include "src/math/atanpif16.h"
+#include "test/UnitTest/FPMatcher.h"
+#include "test/UnitTest/Test.h"
+#include "utils/MPFRWrapper/MPFRUtils.h"
+
+using LlvmLibcAtanpif16Test = LIBC_NAMESPACE::testing::FPTest<float16>;
+
+namespace mpfr = LIBC_NAMESPACE::testing::mpfr;
+
+// Range: [0, Inf]
+static constexpr uint16_t POS_START = 0x0000U;
+static constexpr uint16_t POS_STOP = 0x7c00U;
+
+// Range: [-Inf, 0]
+static constexpr uint16_t NEG_START = 0x8000U;
+static constexpr uint16_t NEG_STOP = 0xfc00U;
+
+TEST_F(LlvmLibcAtanpif16Test, PositiveRange) {
+ for (uint16_t v = POS_START; v <= POS_STOP; ++v) {
+ float16 x = FPBits(v).get_val();
+ EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Atanpi, x,
+ LIBC_NAMESPACE::atanpif16(x), 0.5);
+ }
+}
+
+TEST_F(LlvmLibcAtanpif16Test, NegativeRange) {
+ for (uint16_t v = NEG_START; v <= NEG_STOP; ++v) {
+ float16 x = FPBits(v).get_val();
+ EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Atanpi, x,
+ LIBC_NAMESPACE::atanpif16(x), 0.5);
+ }
+}
\ No newline at end of file
diff --git a/libc/test/src/math/smoke/CMakeLists.txt b/libc/test/src/math/smoke/CMakeLists.txt
index 223d1933bca38..2461478e2ffb5 100644
--- a/libc/test/src/math/smoke/CMakeLists.txt
+++ b/libc/test/src/math/smoke/CMakeLists.txt
@@ -3947,6 +3947,18 @@ add_fp_unittest(
libc.src.__support.FPUtil.fp_bits
)
+add_fp_unittest(
+ atanpif16_test
+ NEED_MPFR
+ SUITE
+ libc-math-unittests
+ SRCS
+ atanpif16_test.cpp
+ DEPENDS
+ libc.src.math.atanpif16
+ libc.src.errno.errno
+)
+
add_fp_unittest(
asinhf_test
SUITE
diff --git a/libc/test/src/math/smoke/atanpif16_test.cpp b/libc/test/src/math/smoke/atanpif16_test.cpp
new file mode 100644
index 0000000000000..e0f3ce5f52392
--- /dev/null
+++ b/libc/test/src/math/smoke/atanpif16_test.cpp
@@ -0,0 +1,66 @@
+//===-- Unittests for atanpif16 -------------------------------------------===//
+//
+// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
+// See https://llvm.org/LICENSE.txt for license information.
+// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
+//
+//===----------------------------------------------------------------------===//
+
+#include "src/errno/libc_errno.h"
+#include "src/math/atanpif16.h"
+#include "test/UnitTest/FPMatcher.h"
+
+using LlvmLibcAtanpif16Test = LIBC_NAMESPACE::testing::FPTest<float16>;
+
+TEST_F(LlvmLibcAtanpif16Test, SpecialNumbers) {
+ // zero
+ EXPECT_FP_EQ(zero, LIBC_NAMESPACE::atanpif16(zero));
+ EXPECT_FP_EQ(neg_zero, LIBC_NAMESPACE::atanpif16(neg_zero));
+
+ // NaN inputs
+ EXPECT_FP_EQ(FPBits::quiet_nan().get_val(),
+ LIBC_NAMESPACE::atanpif16(FPBits::quiet_nan().get_val()));
+
+ EXPECT_FP_EQ(FPBits::quiet_nan().get_val(),
+ LIBC_NAMESPACE::atanpif16(FPBits::signaling_nan().get_val()));
+
+ // infinity inputs -> should return +/-0.5
+ EXPECT_FP_EQ(0.5f16, LIBC_NAMESPACE::atanpif16(inf));
+ EXPECT_FP_EQ(-0.5f16, LIBC_NAMESPACE::atanpif16(neg_inf));
+}
+
+
+TEST_F(LlvmLibcAtanpif16Test, SymmetryProperty) {
+ // Test that atanpi(-x) = -atanpi(x)
+ constexpr float16 TEST_VALS[] = {0.1f16, 0.25f16, 0.5f16, 0.75f16,
+ 1.0f16, 1.5f16, 2.0f16, 5.0f16,
+ 10.0f16, 50.0f16, 100.0f16, 1000.0f16};
+
+ for (float16 x : TEST_VALS) {
+ FPBits neg_x_bits(x);
+ neg_x_bits.set_sign(Sign::NEG);
+ float16 neg_x = neg_x_bits.get_val();
+
+ float16 pos_result = LIBC_NAMESPACE::atanpif16(x);
+ float16 neg_result = LIBC_NAMESPACE::atanpif16(neg_x);
+
+ EXPECT_FP_EQ(pos_result, FPBits(neg_result).abs().get_val());
+ }
+}
+
+TEST_F(LlvmLibcAtanpif16Test, MonotonicityProperty) {
+ // Test that atanpi is monotonically increasing
+ constexpr float16 TEST_VALS[] = {-1000.0f16, -100.0f16, -10.0f16, -2.0f16,
+ -1.0f16, -0.5f16, -0.1f16, 0.0f16,
+ 0.1f16, 0.5f16, 1.0f16, 2.0f16,
+ 10.0f16, 100.0f16, 1000.0f16};
+
+ for (size_t i = 0; i < sizeof(TEST_VALS) / sizeof(TEST_VALS[0]) - 1; ++i) {
+ float16 x1 = TEST_VALS[i];
+ float16 x2 = TEST_VALS[i + 1];
+ float16 result1 = LIBC_NAMESPACE::atanpif16(x1);
+ float16 result2 = LIBC_NAMESPACE::atanpif16(x2);
+
+ EXPECT_TRUE(result1 < result2);
+ }
+}
diff --git a/libc/utils/MPFRWrapper/MPCommon.cpp b/libc/utils/MPFRWrapper/MPCommon.cpp
index ac8cde2f97221..5bb95df02fcf0 100644
--- a/libc/utils/MPFRWrapper/MPCommon.cpp
+++ b/libc/utils/MPFRWrapper/MPCommon.cpp
@@ -106,6 +106,12 @@ MPFRNumber MPFRNumber::atanh() const {
return result;
}
+MPFRNumber MPFRNumber::atanpi() const {
+ MPFRNumber result(*this);
+ mpfr_atanpi(result.value, value, mpfr_rounding);
+ return result;
+}
+
MPFRNumber MPFRNumber::cbrt() const {
MPFRNumber result(*this);
mpfr_cbrt(result.value, value, mpfr_rounding);
diff --git a/libc/utils/MPFRWrapper/MPCommon.h b/libc/utils/MPFRWrapper/MPCommon.h
index eaa512e30bc86..6888fcfad1329 100644
--- a/libc/utils/MPFRWrapper/MPCommon.h
+++ b/libc/utils/MPFRWrapper/MPCommon.h
@@ -187,6 +187,7 @@ class MPFRNumber {
MPFRNumber atan() const;
MPFRNumber atan2(const MPFRNumber &b);
MPFRNumber atanh() const;
+ MPFRNumber atanpi() const;
MPFRNumber cbrt() const;
MPFRNumber ceil() const;
MPFRNumber cos() const;
diff --git a/libc/utils/MPFRWrapper/MPFRUtils.cpp b/libc/utils/MPFRWrapper/MPFRUtils.cpp
index 4a033dbc2049b..cb65f080192b8 100644
--- a/libc/utils/MPFRWrapper/MPFRUtils.cpp
+++ b/libc/utils/MPFRWrapper/MPFRUtils.cpp
@@ -40,6 +40,8 @@ unary_operation(Operation op, InputType input, unsigned int precision,
return mpfrInput.atan();
case Operation::Atanh:
return mpfrInput.atanh();
+ case Operation::Atanpi:
+ return mpfrInput.atanpi();
case Operation::Cbrt:
return mpfrInput.cbrt();
case Operation::Ceil:
diff --git a/libc/utils/MPFRWrapper/MPFRUtils.h b/libc/utils/MPFRWrapper/MPFRUtils.h
index bc65f87c6b5ab..a650a09428f18 100644
--- a/libc/utils/MPFRWrapper/MPFRUtils.h
+++ b/libc/utils/MPFRWrapper/MPFRUtils.h
@@ -32,6 +32,7 @@ enum class Operation : int {
Asinh,
Atan,
Atanh,
+ Atanpi,
Cbrt,
Ceil,
Cos,
|
|
✅ With the latest revision this PR passed the C/C++ code formatter. |
hulxv
force-pushed
the
libc/math/teagt/impl-atanpif16
branch
from
56aefa1 to
1e3a59f
Compare
overmighty
requested changes
Jul 27, 2025
lntue
reviewed
Aug 15, 2025
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Add this suggestion to a batch that can be applied as a single commit.
This suggestion is invalid because no changes were made to the code.
Suggestions cannot be applied while the pull request is closed.
Suggestions cannot be applied while viewing a subset of changes.
Only one suggestion per line can be applied in a batch.
Add this suggestion to a batch that can be applied as a single commit.
Applying suggestions on deleted lines is not supported.
You must change the existing code in this line in order to create a valid suggestion.
Outdated suggestions cannot be applied.
This suggestion has been applied or marked resolved.
Suggestions cannot be applied from pending reviews.
Suggestions cannot be applied on multi-line comments.
Suggestions cannot be applied while the pull request is queued to merge.
Suggestion cannot be applied right now. Please check back later.
This PR implements$\frac{\arctan(x)}{\pi}$ for half-precision floating-point numbers using polynomial approximation with domain reduction.
atanpif16(x)which computesMathematical Implementation
The implementation uses a 15th-degree Taylor polynomial expansion of$\frac{\arctan(x)}{\pi}$ that's computed using $|x| \in [0, 0.5)$ :
python-sympyand it's accurate inTo ensure accuracy across all real inputs, the domain is divided into three cases with appropriate transformations:
Case 1:$|x| \leq 0.5$
Direct polynomial evaluation:
Case 2:$0.5 < |x| \leq 1$
Double-angle reduction using:
Case 3:$|x| > 1$
Reciprocal transformation using
Closes #132212