[libc][math][c23] implement C23 math function asinpif16 #146226
+344
−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) ChangesThe function is implemented using the following Taylor series that's generated using python-sympy, and it is very accurate for x $$ Geogebra graph
Closes #132210 Full diff: https://github.com/llvm/llvm-project/pull/146226.diff 10 Files Affected:
diff --git a/libc/config/linux/aarch64/entrypoints.txt b/libc/config/linux/aarch64/entrypoints.txt
index be5f5a66016b5..240d48f5bcdb8 100644
--- a/libc/config/linux/aarch64/entrypoints.txt
+++ b/libc/config/linux/aarch64/entrypoints.txt
@@ -651,6 +651,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 73dfeae1a2c94..c9e5b0096099f 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.asinpif16
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 947bd4b60b391..aaacc3ca92832 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 | | | | 7.12.4.3 | F.10.1.3 |
+-----------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+----------------------------+
diff --git a/libc/include/math.yaml b/libc/include/math.yaml
index fef829422244d..b49c17b871a93 100644
--- a/libc/include/math.yaml
+++ b/libc/include/math.yaml
@@ -71,7 +71,14 @@ functions:
- stdc
return_type: double
arguments:
- - type: double
+ - 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 d177ff79141c0..e7ccb3726c54c 100644
--- a/libc/src/math/CMakeLists.txt
+++ b/libc/src/math/CMakeLists.txt
@@ -56,6 +56,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)
diff --git a/libc/src/math/asinpif16.h b/libc/src/math/asinpif16.h
new file mode 100644
index 0000000000000..67ccb4ff4ac3d
--- /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_ASINFPI16_H
+#define LLVM_LIBC_SRC_MATH_ASINFPI16_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 adbed5b2de48c..76accf72058d3 100644
--- a/libc/src/math/generic/CMakeLists.txt
+++ b/libc/src/math/generic/CMakeLists.txt
@@ -4017,6 +4017,22 @@ 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.fenv_impl
+ libc.src.__support.FPUtil.fp_bits
+ libc.src.__support.FPUtil.multiply_add
+ 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..34b0b7e7779fc
--- /dev/null
+++ b/libc/src/math/generic/asinpif16.cpp
@@ -0,0 +1,177 @@
+//===-- Half-precision asinf16(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 = 3;
+static constexpr float16 ONE_OVER_SIX = 0.166748046875f16;
+
+static constexpr fputil::ExceptValues<float16, N_ASINFPI_EXCEPTS>
+ 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.0) = 0.5
+ {(fputil::FPBits<float16>(1.0f16)).uintval(),
+ (fputil::FPBits<float16>(ONE_OVER_TWO)).uintval(), 0, 0, 0},
+
+ // x = 0.5, asinfpi(0.5) = 1/6
+ {(fputil::FPBits<float16>(0.5f16)).uintval(),
+ (fputil::FPBits<float16>(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<float16>;
+
+ FPBits xbits(x);
+ uint16_t x_uint = xbits.uintval();
+ bool is_neg = static_cast<bool>(x_uint >> 15);
+ float16 x_abs = is_neg ? -x : x;
+
+ auto __signed_result = [is_neg](float16 r) -> float16 {
+ 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
+
+ // zero
+ if (LIBC_UNLIKELY(x_abs == 0.0f16)) {
+ return 0.0f16;
+ }
+
+ // +/-1.0
+ // if x is +/-1.0, return +/-0.5
+ if (LIBC_UNLIKELY(x_abs == 1.0f16)) {
+ return fputil::cast<float16>(__signed_result(ONE_OVER_TWO));
+ }
+
+ // 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 float16 POLY_COEFFS[10] = {
+ 0.318309886183791f16, // x^1
+ 0.0530516476972984f16, // x^3
+ 0.0238732414637843f16, // x^5
+ 0.0142102627760621f16, // x^7
+ 0.00967087327815336f16, // x^9
+ 0.00712127941391293f16, // x^11
+ 0.00552355646848375f16, // x^13
+ 0.00444514782463692f16, // x^15
+ 0.00367705242846804f16, // x^17
+ 0.00310721681820837f16 // x^19
+ };
+
+ // polynomial evaluation using horner's method
+ // work only for |x| in [0, 0.5]
+ auto __asinpi_polyeval = [](float16 xsq) -> float16 {
+ 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;
+ float16 result = x * __asinpi_polyeval(xsq);
+ return fputil::cast<float16>(__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)
+
+ float16 u = fputil::multiply_add(-ONE_OVER_TWO, x_abs, ONE_OVER_TWO);
+ float16 u_sqrt = fputil::sqrt<float16>(u);
+ float16 asinpi_sqrt_u = u_sqrt * __asinpi_polyeval(u);
+ float16 result = fputil::multiply_add(-2.0f16, asinpi_sqrt_u, ONE_OVER_TWO);
+
+ return fputil::cast<float16>(__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..bd62b76817ceb 100644
--- a/libc/test/src/math/CMakeLists.txt
+++ b/libc/test/src/math/CMakeLists.txt
@@ -2245,6 +2245,19 @@ 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.fabs
+ libc.src.math.asinpif16
+ libc.src.__support.FPUtil.fp_bits
+)
+
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..6d8f76568902e
--- /dev/null
+++ b/libc/test/src/math/asinpif16_test.cpp
@@ -0,0 +1,110 @@
+//===-- 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/errno/libc_errno.h"
+#include "src/math/asinpif16.h"
+#include "src/math/fabs.h"
+#include "test/UnitTest/FPMatcher.h"
+
+#include <errno.h>
+#include <stdint.h>
+
+using LlvmLibcAsinpif16Test = LIBC_NAMESPACE::testing::FPTest<float16>;
+
+TEST_F(LlvmLibcAsinpif16Test, SpecialNumbers) {
+ using FPBits = LIBC_NAMESPACE::fputil::FPBits<float16>;
+
+ // 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<float16>;
+
+ // 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, LIBC_NAMESPACE::fabs(neg_result));
+ }
+ }
+}
+
+TEST_F(LlvmLibcAsinpif16Test, RangeValidation) {
+ // Test that output is always in [-0.5, 0.5] for valid inputs
+ constexpr int num_tests = 1000;
+
+ for (int i = 0; i <= num_tests; ++i) {
+ float16 t = -1.0f16 + (2.0f16 * static_cast<float16>(i)) /
+ static_cast<float16>(num_tests);
+
+ if (LIBC_NAMESPACE::fabs(t) <= 1.0) {
+ float16 result = LIBC_NAMESPACE::asinpif16(t);
+
+ // should be in [-0.5, 0.5]
+ EXPECT_TRUE(result >= -0.5f16);
+ EXPECT_TRUE(result <= 0.5f16);
+ }
+ }
+}
|
|
✅ With the latest revision this PR passed the C/C++ code formatter. |
overmighty
requested changes
Jun 28, 2025
There was a problem hiding this comment.
Please split this test into a smoke test and an MPFR test that checks the entire input range for this function. See tests for existing float16 math functions.
overmighty
reviewed
Jun 28, 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.
The function is implemented using the following Taylor series that's generated using python-sympy, and it is very accurate for |x|$$\in [0, 0.5]$$ and has been verified using Geogebra. The range reduction is used for the rest range (0.5, 1].
Geogebra graph
Closes #132210