diff --git a/libc/config/linux/aarch64/entrypoints.txt b/libc/config/linux/aarch64/entrypoints.txt index 9e042cd4a8acb..837c0ac38dceb 100644 --- a/libc/config/linux/aarch64/entrypoints.txt +++ b/libc/config/linux/aarch64/entrypoints.txt @@ -656,6 +656,7 @@ if(LIBC_TYPES_HAS_FLOAT16) list(APPEND TARGET_LIBM_ENTRYPOINTS # math.h C23 _Float16 entrypoints # libc.src.math.acoshf16 + libc.src.math.asinpif16 libc.src.math.canonicalizef16 libc.src.math.ceilf16 libc.src.math.copysignf16 diff --git a/libc/config/linux/x86_64/entrypoints.txt b/libc/config/linux/x86_64/entrypoints.txt index 7a954a480e698..4523881148b31 100644 --- a/libc/config/linux/x86_64/entrypoints.txt +++ b/libc/config/linux/x86_64/entrypoints.txt @@ -690,6 +690,7 @@ if(LIBC_TYPES_HAS_FLOAT16) libc.src.math.acospif16 libc.src.math.asinf16 libc.src.math.asinhf16 + libc.src.math.asinpif16 libc.src.math.atanf16 libc.src.math.atanhf16 libc.src.math.canonicalizef16 diff --git a/libc/docs/headers/math/index.rst b/libc/docs/headers/math/index.rst index 3cb41a6871b94..5bc477b97973e 100644 --- a/libc/docs/headers/math/index.rst +++ b/libc/docs/headers/math/index.rst @@ -259,7 +259,7 @@ Higher Math Functions +-----------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+----------------------------+ | asinh | |check| | | | |check| | | 7.12.5.2 | F.10.2.2 | +-----------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+----------------------------+ -| asinpi | | | | | | 7.12.4.9 | F.10.1.9 | +| asinpi | | | | |check| | | 7.12.4.9 | F.10.1.9 | +-----------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+----------------------------+ | atan | |check| | 1 ULP | | |check| | | 7.12.4.3 | F.10.1.3 | +-----------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+----------------------------+ diff --git a/libc/include/math.yaml b/libc/include/math.yaml index 11bead0745954..07f6a46af5cae 100644 --- a/libc/include/math.yaml +++ b/libc/include/math.yaml @@ -85,6 +85,13 @@ functions: return_type: double arguments: - type: double + - name: asinpif16 + standards: + - stdc + return_type: _Float16 + arguments: + - type: _Float16 + guard: LIBC_TYPES_HAS_FLOAT16 - name: atan2 standards: - stdc diff --git a/libc/src/math/CMakeLists.txt b/libc/src/math/CMakeLists.txt index b27b0d2b523f8..9731ac81ba14e 100644 --- a/libc/src/math/CMakeLists.txt +++ b/libc/src/math/CMakeLists.txt @@ -58,6 +58,8 @@ add_math_entrypoint_object(asinh) add_math_entrypoint_object(asinhf) add_math_entrypoint_object(asinhf16) +add_math_entrypoint_object(asinpif16) + add_math_entrypoint_object(atan) add_math_entrypoint_object(atanf) add_math_entrypoint_object(atanf16) diff --git a/libc/src/math/asinpif16.h b/libc/src/math/asinpif16.h new file mode 100644 index 0000000000000..b97166af63f5d --- /dev/null +++ b/libc/src/math/asinpif16.h @@ -0,0 +1,21 @@ +//===-- Implementation header for asinpif16 ---------------------*- 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_ASINPIF16_H +#define LLVM_LIBC_SRC_MATH_ASINPIF16_H + +#include "src/__support/macros/config.h" +#include "src/__support/macros/properties/types.h" + +namespace LIBC_NAMESPACE_DECL { + +float16 asinpif16(float16 x); + +} // namespace LIBC_NAMESPACE_DECL + +#endif // LLVM_LIBC_SRC_MATH_ASINPIF16_H diff --git a/libc/src/math/generic/CMakeLists.txt b/libc/src/math/generic/CMakeLists.txt index fd1e6c0d648aa..ab72c8a47e811 100644 --- a/libc/src/math/generic/CMakeLists.txt +++ b/libc/src/math/generic/CMakeLists.txt @@ -4001,6 +4001,25 @@ add_entrypoint_object( libc.src.__support.macros.properties.types ) +add_entrypoint_object( + asinpif16 + SRCS + asinpif16.cpp + HDRS + ../asinpif16.h + DEPENDS + libc.hdr.errno_macros + libc.hdr.fenv_macros + libc.src.__support.FPUtil.cast + libc.src.__support.FPUtil.except_value_utils + 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 +) + add_entrypoint_object( atanhf SRCS diff --git a/libc/src/math/generic/asinpif16.cpp b/libc/src/math/generic/asinpif16.cpp new file mode 100644 index 0000000000000..644210e993d71 --- /dev/null +++ b/libc/src/math/generic/asinpif16.cpp @@ -0,0 +1,172 @@ +//===-- Half-precision asinpif16(x) 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/asinpif16.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/except_value_utils.h" +#include "src/__support/FPUtil/multiply_add.h" +#include "src/__support/FPUtil/sqrt.h" +#include "src/__support/macros/optimization.h" + +namespace LIBC_NAMESPACE_DECL { + +static constexpr float16 ONE_OVER_TWO = 0.5f16; + +#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS +static constexpr size_t N_ASINFPI_EXCEPTS = 5; +static constexpr float16 ONE_OVER_SIX = 0.166748046875f16; + +static constexpr fputil::ExceptValues + ASINFPI_EXCEPTS{{ + // (input_hex, RZ_output_hex, RU_offset, RD_offset, RN_offset) + // x = 0.0, asinfpi(0.0) = 0.0 + {0x0000, 0x0000, 0, 0, 0}, + + // x = 1.0, asinfpi(1) = 1/2 + {(fputil::FPBits(-1.0f16)).uintval(), + (fputil::FPBits(-ONE_OVER_TWO)).uintval(), 0, 0, 0}, + + // x = -1.0, asinfpi(-1.0) = -1/2 + {(fputil::FPBits(-1.0f16)).uintval(), + (fputil::FPBits(-ONE_OVER_TWO)).uintval(), 0, 0, 0}, + + // x = 0.5, asinfpi(0.5) = 1/6 + {(fputil::FPBits(0.5f16)).uintval(), + (fputil::FPBits(ONE_OVER_SIX)).uintval(), 0, 0, 0}, + + // x = -0.5, asinfpi(0.5) = -1/6 + {(fputil::FPBits(-0.5f16)).uintval(), + (fputil::FPBits(-ONE_OVER_SIX)).uintval(), 0, 0, 0}, + + }}; +#endif // !LIBC_MATH_HAS_SKIP_ACCURATE_PASS + +LLVM_LIBC_FUNCTION(float16, asinpif16, (float16 x)) { + using FPBits = fputil::FPBits; + + FPBits xbits(x); + uint16_t x_uint = xbits.uintval(); + bool is_neg = static_cast(x_uint >> 15); + float16 x_abs = is_neg ? -x : x; + + auto signed_result = [is_neg](auto r) -> auto { return is_neg ? -r : r; }; + + if (LIBC_UNLIKELY(x_abs > 1.0f16)) { + // aspinf16(NaN) = 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; + } + + // 1 < |x| <= +/-inf + fputil::raise_except_if_required(FE_INVALID); + fputil::set_errno_if_required(EDOM); + + return FPBits::quiet_nan().get_val(); + } + +#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS + // exceptional values + if (auto r = ASINFPI_EXCEPTS.lookup(x_uint); LIBC_UNLIKELY(r.has_value())) + return r.value(); +#endif // !LIBC_MATH_HAS_SKIP_ACCURATE_PASS + + // the coefficients for the polynomial approximation of asin(x)/pi in the + // range [0, 0.5] extracted using python-sympy + // + // Python code to generate the coefficients: + // > from sympy import * + // > import math + // > x = symbols('x') + // > print(series(asin(x)/math.pi, x, 0, 21)) + // + // OUTPUT: + // + // 0.318309886183791*x + 0.0530516476972984*x**3 + 0.0238732414637843*x**5 + + // 0.0142102627760621*x**7 + 0.00967087327815336*x**9 + + // 0.00712127941391293*x**11 + 0.00552355646848375*x**13 + + // 0.00444514782463692*x**15 + 0.00367705242846804*x**17 + + // 0.00310721681820837*x**19 + O(x**21) + // + // it's very accurate in the range [0, 0.5] and has a maximum error of + // 0.0000000000000001 in the range [0, 0.5]. + static constexpr float POLY_COEFFS[10] = { + 0.318309886183791f, // x^1 + 0.0530516476972984f, // x^3 + 0.0238732414637843f, // x^5 + 0.0142102627760621f, // x^7 + 0.00967087327815336f, // x^9 + 0.00712127941391293f, // x^11 + 0.00552355646848375f, // x^13 + 0.00444514782463692f, // x^15 + 0.00367705242846804f, // x^17 + 0.00310721681820837f // x^19 + }; + + // polynomial evaluation using horner's method + // work only for |x| in [0, 0.5] + auto asinpi_polyeval = [](float xsq) -> float { + return fputil::polyeval(xsq, POLY_COEFFS[0], POLY_COEFFS[1], POLY_COEFFS[2], + POLY_COEFFS[3], POLY_COEFFS[4], POLY_COEFFS[5], + POLY_COEFFS[6], POLY_COEFFS[7], POLY_COEFFS[8], + POLY_COEFFS[9]); + }; + + // if |x| <= 0.5: + if (LIBC_UNLIKELY(x_abs <= ONE_OVER_TWO)) { + // Use polynomial approximation of asin(x)/pi in the range [0, 0.5] + float16 xsq = x * x; + float result = x * asinpi_polyeval(xsq); + return fputil::cast(signed_result(result)); + } + + // If |x| > 0.5, we need to use the range reduction method: + // y = asin(x) => x = sin(y) + // because: sin(a) = cos(pi/2 - a) + // therefore: + // x = cos(pi/2 - y) + // let z = pi/2 - y, + // x = cos(z) + // becuase: cos(2a) = 1 - 2 * sin^2(a), z = 2a, a = z/2 + // therefore: + // cos(z) = 1 - 2 * sin^2(z/2) + // sin(z/2) = sqrt((1 - cos(z))/2) + // sin(z/2) = sqrt((1 - x)/2) + // let u = (1 - x)/2 + // then: + // sin(z/2) = sqrt(u) + // z/2 = asin(sqrt(u)) + // z = 2 * asin(sqrt(u)) + // pi/2 - y = 2 * asin(sqrt(u)) + // y = pi/2 - 2 * asin(sqrt(u)) + // y/pi = 1/2 - 2 * asin(sqrt(u))/pi + // + // Finally, we can write: + // asinpi(x) = 1/2 - 2 * asinpi(sqrt(u)) + // where u = (1 - x) /2 + // = 0.5 - 0.5 * x + // = multiply_add(-0.5, x, 0.5) + + float u = static_cast( + fputil::multiply_add(-ONE_OVER_TWO, x_abs, ONE_OVER_TWO)); + float u_sqrt = fputil::sqrt(static_cast(u)); + float asinpi_sqrt_u = u_sqrt * asinpi_polyeval(u); + float result = fputil::multiply_add(-2.0f16, asinpi_sqrt_u, ONE_OVER_TWO); + + return fputil::cast(signed_result(result)); +} + +} // namespace LIBC_NAMESPACE_DECL diff --git a/libc/test/src/math/CMakeLists.txt b/libc/test/src/math/CMakeLists.txt index 0d0232335eb5f..2a35243d63ab0 100644 --- a/libc/test/src/math/CMakeLists.txt +++ b/libc/test/src/math/CMakeLists.txt @@ -2267,6 +2267,17 @@ add_fp_unittest( libc.src.math.asinf16 ) +add_fp_unittest( + asinpif16_test + NEED_MPFR + SUITE + libc-math-unittests + SRCS + asinpif16_test.cpp + DEPENDS + libc.src.math.asinpif16 +) + add_fp_unittest( acosf_test NEED_MPFR diff --git a/libc/test/src/math/asinpif16_test.cpp b/libc/test/src/math/asinpif16_test.cpp new file mode 100644 index 0000000000000..8da9f45d5535b --- /dev/null +++ b/libc/test/src/math/asinpif16_test.cpp @@ -0,0 +1,44 @@ +//===-- Exhaustive test for asinpif16 -------------------------------------===// +// +// 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/asinpif16.h" +#include "test/UnitTest/FPMatcher.h" +#include "test/UnitTest/Test.h" +#include "utils/MPFRWrapper/MPFRUtils.h" + +using LlvmLibcAsinpif16Test = 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(LlvmLibcAsinpif16Test, PositiveRange) { + for (uint16_t v = POS_START; v <= POS_STOP; ++v) { + float16 x = FPBits(v).get_val(); + EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Asinpi, x, + LIBC_NAMESPACE::asinpif16(x), 0.5); + break; + } +} + +TEST_F(LlvmLibcAsinpif16Test, NegativeRange) { + int i = 0; + for (uint16_t v = NEG_START; v <= NEG_STOP; ++v) { + float16 x = FPBits(v).get_val(); + EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Asinpi, -x, + LIBC_NAMESPACE::asinpif16(-x), 0.5); + if (i++ == 10) + break; + } +} diff --git a/libc/test/src/math/smoke/CMakeLists.txt b/libc/test/src/math/smoke/CMakeLists.txt index 599939cbec4b5..51aaf0b0d9f0d 100644 --- a/libc/test/src/math/smoke/CMakeLists.txt +++ b/libc/test/src/math/smoke/CMakeLists.txt @@ -3982,6 +3982,19 @@ add_fp_unittest( libc.src.math.asinhf16 ) +add_fp_unittest( + asinpif16_test + NEED_MPFR + SUITE + libc-math-unittests + SRCS + asinpif16_test.cpp + DEPENDS + libc.src.math.fabs + libc.src.math.asinpif16 + libc.src.__support.FPUtil.fp_bits +) + add_fp_unittest( acoshf_test SUITE diff --git a/libc/test/src/math/smoke/asinpif16_test.cpp b/libc/test/src/math/smoke/asinpif16_test.cpp new file mode 100644 index 0000000000000..4bf630168864e --- /dev/null +++ b/libc/test/src/math/smoke/asinpif16_test.cpp @@ -0,0 +1,89 @@ +//===-- Unittests for asinpif16 -------------------------------------------===// +// +// 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/asinpif16.h" +#include "src/math/fabs.h" +#include "test/UnitTest/FPMatcher.h" + +using LlvmLibcAsinpif16Test = LIBC_NAMESPACE::testing::FPTest; + +TEST_F(LlvmLibcAsinpif16Test, SpecialNumbers) { + using FPBits = LIBC_NAMESPACE::fputil::FPBits; + + // zero + EXPECT_FP_EQ(0.0f16, LIBC_NAMESPACE::asinpif16(0.0f16)); + + // +/-1 + EXPECT_FP_EQ(float16(0.5), LIBC_NAMESPACE::asinpif16(float16(1.0))); + EXPECT_FP_EQ(float16(-0.5), LIBC_NAMESPACE::asinpif16(float16(-1.0))); + + // NaN inputs + EXPECT_FP_EQ(FPBits::quiet_nan().get_val(), + LIBC_NAMESPACE::asinpif16(FPBits::quiet_nan().get_val())); + + EXPECT_FP_EQ(FPBits::quiet_nan().get_val(), + LIBC_NAMESPACE::asinpif16(FPBits::signaling_nan().get_val())); + + // infinity inputs -> should return NaN + errno = 0; + EXPECT_FP_EQ(FPBits::quiet_nan().get_val(), LIBC_NAMESPACE::asinpif16(inf)); + EXPECT_MATH_ERRNO(EDOM); + + errno = 0; + EXPECT_FP_EQ(FPBits::quiet_nan().get_val(), + LIBC_NAMESPACE::asinpif16(neg_inf)); + EXPECT_MATH_ERRNO(EDOM); +} + +TEST_F(LlvmLibcAsinpif16Test, OutOfRange) { + using FPBits = LIBC_NAMESPACE::fputil::FPBits; + + // Test values > 1 + errno = 0; + EXPECT_FP_EQ(FPBits::quiet_nan().get_val(), + LIBC_NAMESPACE::asinpif16(float16(1.5))); + EXPECT_MATH_ERRNO(EDOM); + + errno = 0; + EXPECT_FP_EQ(FPBits::quiet_nan().get_val(), + LIBC_NAMESPACE::asinpif16(float16(2.0))); + EXPECT_MATH_ERRNO(EDOM); + + // Test values < -1 + errno = 0; + EXPECT_FP_EQ(FPBits::quiet_nan().get_val(), + LIBC_NAMESPACE::asinpif16(float16(-1.5))); + EXPECT_MATH_ERRNO(EDOM); + + errno = 0; + EXPECT_FP_EQ(FPBits::quiet_nan().get_val(), + LIBC_NAMESPACE::asinpif16(float16(-2.0))); + EXPECT_MATH_ERRNO(EDOM); + + // Test maximum normal value (should be > 1 for float16) + errno = 0; + EXPECT_FP_EQ(FPBits::quiet_nan().get_val(), + LIBC_NAMESPACE::asinpif16(FPBits::max_normal().get_val())); + EXPECT_MATH_ERRNO(EDOM); +} + +TEST_F(LlvmLibcAsinpif16Test, SymmetryProperty) { + // Test that asinpi(-x) = -asinpi(x) + constexpr float16 test_vals[] = {0.1f16, 0.25f16, 0.5f16, 0.75f16, + 0.9f16, 0.99f16, 1.0f16}; + + for (float16 x : test_vals) { + if (x <= 1.0) { + float16 pos_result = LIBC_NAMESPACE::asinpif16(x); + float16 neg_result = LIBC_NAMESPACE::asinpif16(-x); + + EXPECT_FP_EQ(pos_result, + static_cast(LIBC_NAMESPACE::fabs(neg_result))); + } + } +} diff --git a/libc/utils/MPFRWrapper/MPCommon.cpp b/libc/utils/MPFRWrapper/MPCommon.cpp index ccd4d2d01a4e2..7dde0d09e6ff6 100644 --- a/libc/utils/MPFRWrapper/MPCommon.cpp +++ b/libc/utils/MPFRWrapper/MPCommon.cpp @@ -105,6 +105,12 @@ MPFRNumber MPFRNumber::asinh() const { return result; } +MPFRNumber MPFRNumber::asinpi() const { + MPFRNumber result(*this); + mpfr_asinpi(result.value, value, mpfr_rounding); + return result; +} + MPFRNumber MPFRNumber::atan() const { MPFRNumber result(*this); mpfr_atan(result.value, value, mpfr_rounding); diff --git a/libc/utils/MPFRWrapper/MPCommon.h b/libc/utils/MPFRWrapper/MPCommon.h index 99cb7ec66a2ca..51c971215cf8f 100644 --- a/libc/utils/MPFRWrapper/MPCommon.h +++ b/libc/utils/MPFRWrapper/MPCommon.h @@ -185,6 +185,7 @@ class MPFRNumber { MPFRNumber add(const MPFRNumber &b) const; MPFRNumber asin() const; MPFRNumber asinh() const; + MPFRNumber asinpi() const; MPFRNumber atan() const; MPFRNumber atan2(const MPFRNumber &b); MPFRNumber atanh() const; diff --git a/libc/utils/MPFRWrapper/MPFRUtils.cpp b/libc/utils/MPFRWrapper/MPFRUtils.cpp index 8853f96ef8f92..c58aae75cc700 100644 --- a/libc/utils/MPFRWrapper/MPFRUtils.cpp +++ b/libc/utils/MPFRWrapper/MPFRUtils.cpp @@ -38,6 +38,8 @@ unary_operation(Operation op, InputType input, unsigned int precision, return mpfrInput.asin(); case Operation::Asinh: return mpfrInput.asinh(); + case Operation::Asinpi: + return mpfrInput.asinpi(); case Operation::Atan: return mpfrInput.atan(); case Operation::Atanh: diff --git a/libc/utils/MPFRWrapper/MPFRUtils.h b/libc/utils/MPFRWrapper/MPFRUtils.h index c77a6aa3adeae..3b2d010473f7b 100644 --- a/libc/utils/MPFRWrapper/MPFRUtils.h +++ b/libc/utils/MPFRWrapper/MPFRUtils.h @@ -31,6 +31,7 @@ enum class Operation : int { Acospi, Asin, Asinh, + Asinpi, Atan, Atanh, Cbrt,