diff --git a/libc/config/linux/x86_64/entrypoints.txt b/libc/config/linux/x86_64/entrypoints.txt index e8f59c9d5cd73..91bfaa70ad748 100644 --- a/libc/config/linux/x86_64/entrypoints.txt +++ b/libc/config/linux/x86_64/entrypoints.txt @@ -705,6 +705,7 @@ if(LIBC_TYPES_HAS_FLOAT16) libc.src.math.asinhf16 libc.src.math.atanf16 libc.src.math.atanhf16 + libc.src.math.atanpif16 libc.src.math.canonicalizef16 libc.src.math.ceilf16 libc.src.math.copysignf16 diff --git a/libc/docs/headers/math/index.rst b/libc/docs/headers/math/index.rst index 9679c4a6c807f..d5fb397342d4f 100644 --- a/libc/docs/headers/math/index.rst +++ b/libc/docs/headers/math/index.rst @@ -269,7 +269,7 @@ Higher Math Functions +-----------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+----------------------------+ | atanh | |check| | | | |check| | | 7.12.5.3 | F.10.2.3 | +-----------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+----------------------------+ -| atanpi | | | | | | 7.12.4.10 | F.10.1.10 | +| atanpi | | | | |check| | | 7.12.4.10 | F.10.1.10 | +-----------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+----------------------------+ | cbrt | |check| | |check| | | | | 7.12.7.1 | F.10.4.1 | +-----------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+----------------------------+ diff --git a/libc/include/math.yaml b/libc/include/math.yaml index 007be235f4380..539e42c0bbecb 100644 --- a/libc/include/math.yaml +++ b/libc/include/math.yaml @@ -134,6 +134,13 @@ functions: arguments: - type: float - name: atanhf16 + standatds: + - stdc + return_type: _Float16 + arguments: + - type: _Float16 + guard: LIBC_TYPES_HAS_FLOAT16 + - name: atanpif16 standards: - stdc return_type: _Float16 diff --git a/libc/src/math/CMakeLists.txt b/libc/src/math/CMakeLists.txt index 455ad3456573a..8ed6668ce71c7 100644 --- a/libc/src/math/CMakeLists.txt +++ b/libc/src/math/CMakeLists.txt @@ -71,6 +71,8 @@ add_math_entrypoint_object(atanh) add_math_entrypoint_object(atanhf) add_math_entrypoint_object(atanhf16) +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..8f2391a6e47e5 --- /dev/null +++ b/libc/src/math/atanpif16.h @@ -0,0 +1,21 @@ +//===-- 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 a001d99e032d6..b1fb6aa562288 100644 --- a/libc/src/math/generic/CMakeLists.txt +++ b/libc/src/math/generic/CMakeLists.txt @@ -3951,6 +3951,25 @@ add_entrypoint_object( libc.src.__support.macros.properties.types ) +add_entrypoint_object( + atanpif16 + SRCS + atanpif16.cpp + HDRS + ../atanpif16.h + DEPENDS + libc.hdr.errno_macros + libc.hdr.fenv_macros + libc.src.__support.FPUtil.cast + libc.src.__support.FPUtil.fenv_impl + libc.src.__support.FPUtil.fp_bits + libc.src.__support.FPUtil.multiply_add + libc.src.__support.FPUtil.polyeval + libc.src.__support.FPUtil.sqrt + libc.src.__support.macros.optimization + libc.src.__support.macros.properties.types +) + add_entrypoint_object( asinf SRCS diff --git a/libc/src/math/generic/atanpif16.cpp b/libc/src/math/generic/atanpif16.cpp new file mode 100644 index 0000000000000..4db8be09bfd27 --- /dev/null +++ b/libc/src/math/generic/atanpif16.cpp @@ -0,0 +1,158 @@ +//===-- 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" +#include "src/__support/macros/optimization.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; + + FPBits xbits(x); + bool is_neg = xbits.is_neg(); + + auto signed_result = [is_neg](double r) -> float16 { + return fputil::cast(is_neg ? -r : r); + }; + + if (LIBC_UNLIKELY(xbits.is_inf_or_nan())) { + if (xbits.is_nan()) { + if (xbits.is_signaling_nan()) { + fputil::raise_except_if_required(FE_INVALID); + return FPBits::quiet_nan().get_val(); + } + return x; + } + // atanpi(±∞) = ±0.5 + return signed_result(0.5); + } + + if (LIBC_UNLIKELY(xbits.is_zero())) + return x; + + double x_abs = fputil::cast(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 x_sq = x * x; + return x * fputil::polyeval(x_sq, 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(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 x_sq_recip = x_recip * x_recip; + double sqrt_term = fputil::sqrt(1.0 + x_sq_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 43cde0d68873e..da355e392ed7c 100644 --- a/libc/test/src/math/CMakeLists.txt +++ b/libc/test/src/math/CMakeLists.txt @@ -2169,6 +2169,17 @@ add_fp_unittest( libc.src.math.atanhf16 ) +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..38771f0147d60 --- /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; + +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); + } +} diff --git a/libc/test/src/math/smoke/CMakeLists.txt b/libc/test/src/math/smoke/CMakeLists.txt index ec4c09c5b70f7..8101030bb883e 100644 --- a/libc/test/src/math/smoke/CMakeLists.txt +++ b/libc/test/src/math/smoke/CMakeLists.txt @@ -3965,6 +3965,18 @@ add_fp_unittest( libc.src.__support.FPUtil.cast ) +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..5accb5d34193a --- /dev/null +++ b/libc/test/src/math/smoke/atanpif16_test.cpp @@ -0,0 +1,64 @@ +//===-- 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/__support/libc_errno.h" +#include "src/math/atanpif16.h" +#include "test/UnitTest/FPMatcher.h" + +using LIBC_NAMESPACE::cpp::array; +using LlvmLibcAtanpif16Test = LIBC_NAMESPACE::testing::FPTest; + +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(aNaN, LIBC_NAMESPACE::atanpif16(aNaN)); + + // 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 array 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 < TEST_VALS.size() - 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 07339a06fff80..8a164b8dc25ad 100644 --- a/libc/utils/MPFRWrapper/MPCommon.cpp +++ b/libc/utils/MPFRWrapper/MPCommon.cpp @@ -123,6 +123,21 @@ MPFRNumber MPFRNumber::atanh() const { return result; } +MPFRNumber MPFRNumber::atanpi() const { + MPFRNumber result(*this); +#if MPFR_VERSION >= MPFR_VERSION_NUM(4, 2, 0) + mpfr_atanpi(result.value, value, mpfr_rounding); + return result; +#else + MPFRNumber value_atan(0.0, 1280); + mpfr_atan(value_atan.value, value, MPFR_RNDN); + MPFRNumber value_pi(0.0, 1280); + mpfr_const_pi(value_pi.value, MPFR_RNDN); + mpfr_div(result.value, value_atan.value, value_pi.value, mpfr_rounding); + return result; +#endif +} + 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 8bcc69c247a34..5932b8fff6d60 100644 --- a/libc/utils/MPFRWrapper/MPCommon.h +++ b/libc/utils/MPFRWrapper/MPCommon.h @@ -192,6 +192,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 8853f96ef8f92..2b0dd6e883c31 100644 --- a/libc/utils/MPFRWrapper/MPFRUtils.cpp +++ b/libc/utils/MPFRWrapper/MPFRUtils.cpp @@ -42,6 +42,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 45468c6cb19ae..338253a8cbf8b 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,