[libc][math][c23] Add asinf16() function #124212
Merged
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
Contributor
wldfngrs
commented
Jan 23, 2025
wldfngrs
commented
- Implementation of asin() function for 16-bit floating point inputs.
- Exhaustive tests across the 16-bit input range
|
@llvm/pr-subscribers-libc Author: wldfngrs (wldfngrs) Changes
Full diff: https://github.com./llvm/llvm-project/pull/124212.diff 11 Files Affected:
diff --git a/libc/config/linux/x86_64/entrypoints.txt b/libc/config/linux/x86_64/entrypoints.txt
index 0c1ae9561a7e69..f784432cb25f9e 100644
--- a/libc/config/linux/x86_64/entrypoints.txt
+++ b/libc/config/linux/x86_64/entrypoints.txt
@@ -641,6 +641,7 @@ endif()
if(LIBC_TYPES_HAS_FLOAT16)
list(APPEND TARGET_LIBM_ENTRYPOINTS
# math.h C23 _Float16 entrypoints
+ libc.src.math.asinf16
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 8548e4a5773bc0..3e45e3e618abbb 100644
--- a/libc/docs/headers/math/index.rst
+++ b/libc/docs/headers/math/index.rst
@@ -256,7 +256,7 @@ Higher Math Functions
+-----------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+----------------------------+
| acospi | | | | | | 7.12.4.8 | F.10.1.8 |
+-----------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+----------------------------+
-| asin | |check| | | | | | 7.12.4.2 | F.10.1.2 |
+| asin | |check| | | | |check| | | 7.12.4.2 | F.10.1.2 |
+-----------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+----------------------------+
| asinh | |check| | | | | | 7.12.5.2 | F.10.2.2 |
+-----------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+----------------------------+
diff --git a/libc/include/math.yaml b/libc/include/math.yaml
index 3a660a59d3605f..b98bc55f6cc537 100644
--- a/libc/include/math.yaml
+++ b/libc/include/math.yaml
@@ -32,6 +32,13 @@ functions:
return_type: float
arguments:
- type: float
+ - name: asinf16
+ standards:
+ - stdc
+ return_type: _Float16
+ arguments:
+ - type: _Float16
+ guard: LIBC_TYPES_HAS_FLOAT16
- name: asinhf
standards:
- stdc
diff --git a/libc/src/math/CMakeLists.txt b/libc/src/math/CMakeLists.txt
index fe5ebd793b40af..87004bf6ddd0a3 100644
--- a/libc/src/math/CMakeLists.txt
+++ b/libc/src/math/CMakeLists.txt
@@ -49,6 +49,7 @@ add_math_entrypoint_object(asin)
add_math_entrypoint_object(asinf)
add_math_entrypoint_object(asinh)
add_math_entrypoint_object(asinhf)
+add_math_entrypoint_object(asinf16)
add_math_entrypoint_object(atan)
add_math_entrypoint_object(atanf)
diff --git a/libc/src/math/asinf16.h b/libc/src/math/asinf16.h
new file mode 100644
index 00000000000000..f16647ec2a6f96
--- /dev/null
+++ b/libc/src/math/asinf16.h
@@ -0,0 +1,21 @@
+//===-- Implementation header for asinf16 -----------------------*- 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_ASINF16_H
+#define LLVM_LIBC_SRC_MATH_ASINF16_H
+
+#include "src/__support/macros/config.h"
+#include "src/__support/macros/properties/types.h"
+
+namespace LIBC_NAMESPACE_DECL {
+
+float16 asinf16(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 0e57051807b332..cf77987ff3add9 100644
--- a/libc/src/math/generic/CMakeLists.txt
+++ b/libc/src/math/generic/CMakeLists.txt
@@ -4598,6 +4598,25 @@ add_entrypoint_object(
${libc_opt_high_flag}
)
+add_entrypoint_object(
+ asinf16
+ SRCS
+ asinf16.cpp
+ HDRS
+ ../asinf16.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(
acosf
SRCS
diff --git a/libc/src/math/generic/asinf16.cpp b/libc/src/math/generic/asinf16.cpp
new file mode 100644
index 00000000000000..d340a7244cf34d
--- /dev/null
+++ b/libc/src/math/generic/asinf16.cpp
@@ -0,0 +1,126 @@
+//===-- 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/asinf16.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 {
+
+// Generated by Sollya using the following command:
+// > round(pi/2, D, RN);
+static constexpr float PI_2 = 0x1.921fb54442d18p0f;
+
+LLVM_LIBC_FUNCTION(float16, asinf16, (float16 x)) {
+ using FPBits = fputil::FPBits<float16>;
+ FPBits xbits(x);
+
+ uint16_t x_u = xbits.uintval();
+ uint16_t x_abs = x_u & 0x7fff;
+ float xf = x;
+
+ float xsq = xf * xf;
+
+ // |x| <= 0x1p-1, |x| <= 0.5
+ if (x_abs <= 0x3800) {
+ if (LIBC_UNLIKELY(x_abs <= 0x1a1e)) {
+ // asinf16(+/-0) = +/-0
+ if (LIBC_UNLIKELY(x_abs == 0U)) {
+ return x;
+ }
+
+ // Exhaustive tests show that, when:
+ // x > 0, and rounding upward, or
+ // x < 0, and rounding downward, then,
+ // asin(x) = x * 2^-11 + x
+ // else, in other rounding modes,
+ // asin(x) = x
+ int rounding = fputil::quick_get_round();
+
+ if ((xbits.is_pos() && rounding == FE_UPWARD) ||
+ (xbits.is_neg() && rounding == FE_DOWNWARD)) {
+ return fputil::cast<float16>(fputil::multiply_add(xf, 0x1.0p-11f, xf));
+ }
+ return x;
+ }
+
+ // Degree-6 minimax odd polynomial of asin(x) generated by Sollya with:
+ // > P = fpminimax(asin(x)/x, [|0, 2, 4, 6, 8|], [|SG...|], [0, 0.5]);
+ float result =
+ fputil::polyeval(xsq, 0x1.000002p0f, 0x1.554c2ap-3f, 0x1.3541ccp-4f,
+ 0x1.43b2d6p-5f, 0x1.a0d73ep-5f);
+ return fputil::cast<float16>(xf * result);
+ }
+
+ // |x| > 1, asinf16(x) = NaN
+ if (LIBC_UNLIKELY(x_abs > 0x3c00)) {
+ // |x| <= +/-inf
+ if (LIBC_UNLIKELY(x_abs <= 0x7c00)) {
+ fputil::set_errno_if_required(EDOM);
+ fputil::raise_except_if_required(FE_INVALID);
+ }
+
+ return FPBits::quiet_nan().get_val();
+ }
+
+ // When |x| > 0.5, assume that 0.5 < |x| <= 1,
+ //
+ // Step-by-step range-reduction proof:
+ // 1: Let y = asin(x), such that, x = sin(y)
+ // 2: From complimentary angle identity:
+ // x = sin(y) = cos(pi/2 - y)
+ // 3: Let z = pi/2 - y, such that x = cos(z)
+ // 4: From double angle formula; cos(2A) = 1 - sin^2(A):
+ // z = 2A, z/2 = A
+ // cos(z) = 1 - 2 * sin^2(z/2)
+ // 5: Make sin(z/2) subject of the formula:
+ // sin(z/2) = sqrt((1 - cos(z))/2)
+ // 6: Recall [3]; x = cos(z). Therefore:
+ // sin(z/2) = sqrt((1 - x)/2)
+ // 7: Let u = (1 - x)/2
+ // 8: Therefore:
+ // arcsin(sqrt(u)) = z/2
+ // 2 * arcsin(sqrt(u)) = z
+ // 9: Recall [3], z = pi/2 - y. Therefore:
+ // y = pi/2 - z
+ // y = pi/2 - 2 * arcsin(sqrt(u))
+ // 10: Recall [1], y = asin(x). Therefore:
+ // arcsin(x) = pi/2 - 2 * arcsin(sqrt(u))
+ //
+ // WHY?
+ // 11: Recall [7], u = (1 - x)/2
+ // 12: Since 0.5 < x <= 1, therefore:
+ // 0 <= u <= 0.25 and 0 <= sqrt(u) <= 0.5
+ //
+ // Hence, we can reuse the same [0, 0.5] domain polynomial approximation for
+ // Step [10] as `sqrt(u)` is in range.
+
+ // 0x1p-1 < |x| <= 0x1p0, 0.5 < |x| <= 1.0
+ float xf_abs = (xf < 0 ? -xf : xf);
+ float sign = ((xbits.uintval() >> 15) == 1 ? -1.0 : 1.0);
+ float u = fputil::multiply_add(-0.5f, xf_abs, 0.5f);
+ float u_sqrt = fputil::sqrt<float>(u);
+
+ // Degree-6 minimax odd polynomial of asin(x) generated by Sollya with:
+ // > P = fpminimax(asin(x)/x, [|0, 2, 4, 6, 8|], [|SG...|], [0, 0.5]);
+ float asin_sqrt_u =
+ u_sqrt * fputil::polyeval(u, 0x1.000002p0f, 0x1.554c2ap-3f,
+ 0x1.3541ccp-4f, 0x1.43b2d6p-5f, 0x1.a0d73ep-5f);
+
+ return fputil::cast<float16>(sign *
+ fputil::multiply_add(-2.0f, asin_sqrt_u, PI_2));
+}
+
+} // namespace LIBC_NAMESPACE_DECL
diff --git a/libc/test/src/math/CMakeLists.txt b/libc/test/src/math/CMakeLists.txt
index ae8518ee4b4cc1..06c505037426af 100644
--- a/libc/test/src/math/CMakeLists.txt
+++ b/libc/test/src/math/CMakeLists.txt
@@ -2186,6 +2186,17 @@ add_fp_unittest(
libc.src.__support.FPUtil.fp_bits
)
+add_fp_unittest(
+ asinf16_test
+ NEED_MPFR
+ SUITE
+ libc-math-unittests
+ SRCS
+ asinf16_test.cpp
+ DEPENDS
+ libc.src.math.asinf16
+)
+
add_fp_unittest(
acosf_test
NEED_MPFR
diff --git a/libc/test/src/math/asinf16_test.cpp b/libc/test/src/math/asinf16_test.cpp
new file mode 100644
index 00000000000000..9593cad16ac772
--- /dev/null
+++ b/libc/test/src/math/asinf16_test.cpp
@@ -0,0 +1,42 @@
+//===-- Exhaustive test for asinf16 ---------------------------------------===//
+//
+// 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/asinf16.h"
+#include "test/UnitTest/FPMatcher.h"
+#include "test/UnitTest/Test.h"
+#include "utils/MPFRWrapper/MPFRUtils.h"
+
+using LlvmLibcAsinf16Test = 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(LlvmLibcAsinf16Test, PositiveRange) {
+ for (uint16_t v = POS_START; v <= POS_STOP; ++v) {
+ float16 x = FPBits(v).get_val();
+
+ EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Asin, x,
+ LIBC_NAMESPACE::asinf16(x), 0.5);
+ }
+}
+
+TEST_F(LlvmLibcAsinf16Test, NegativeRange) {
+ for (uint16_t v = NEG_START; v <= NEG_STOP; ++v) {
+ float16 x = FPBits(v).get_val();
+
+ EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Asin, x,
+ LIBC_NAMESPACE::asinf16(x), 0.5);
+ }
+}
diff --git a/libc/test/src/math/smoke/CMakeLists.txt b/libc/test/src/math/smoke/CMakeLists.txt
index e23e7f41222d4a..7af12e133ed236 100644
--- a/libc/test/src/math/smoke/CMakeLists.txt
+++ b/libc/test/src/math/smoke/CMakeLists.txt
@@ -3954,6 +3954,17 @@ add_fp_unittest(
libc.src.__support.FPUtil.fp_bits
)
+add_fp_unittest(
+ asinf16_test
+ SUITE
+ libc-math-smoke-tests
+ SRCS
+ asinf16_test.cpp
+ DEPENDS
+ libc.src.errno.errno
+ libc.src.math.asinf16
+)
+
add_fp_unittest(
acosf_test
SUITE
diff --git a/libc/test/src/math/smoke/asinf16_test.cpp b/libc/test/src/math/smoke/asinf16_test.cpp
new file mode 100644
index 00000000000000..303a3a302c2d5b
--- /dev/null
+++ b/libc/test/src/math/smoke/asinf16_test.cpp
@@ -0,0 +1,40 @@
+//===-- Unittests for asinf16 ---------------------------------------------===//
+//
+//
+// 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/asinf16.h"
+#include "test/UnitTest/FPMatcher.h"
+#include "test/UnitTest/Test.h"
+
+using LlvmLibcAsinf16Test = LIBC_NAMESPACE::testing::FPTest<float16>;
+
+TEST_F(LlvmLibcAsinf16Test, SpecialNumbers) {
+ LIBC_NAMESPACE::libc_errno = 0;
+
+ EXPECT_FP_EQ(aNaN, LIBC_NAMESPACE::asinf16(aNaN));
+ EXPECT_MATH_ERRNO(0);
+
+ EXPECT_FP_EQ(zero, LIBC_NAMESPACE::asinf16(zero));
+ EXPECT_MATH_ERRNO(0);
+
+ EXPECT_FP_EQ(neg_zero, LIBC_NAMESPACE::asinf16(neg_zero));
+ EXPECT_MATH_ERRNO(0);
+
+ EXPECT_FP_EQ(aNaN, LIBC_NAMESPACE::asinf16(inf));
+ EXPECT_MATH_ERRNO(EDOM);
+
+ EXPECT_FP_EQ(aNaN, LIBC_NAMESPACE::asinf16(neg_inf));
+ EXPECT_MATH_ERRNO(EDOM);
+
+ EXPECT_FP_EQ(aNaN, LIBC_NAMESPACE::asinf16(2.0f));
+ EXPECT_MATH_ERRNO(EDOM);
+
+ EXPECT_FP_EQ(aNaN, LIBC_NAMESPACE::asinf16(-2.0f));
+ EXPECT_MATH_ERRNO(EDOM);
+}
|
|
@overmighty @lntue please review |
overmighty
requested changes
Jan 27, 2025
libc/src/math/generic/asinf16.cpp
Outdated
Comment on lines
111
to
112
| float xf_abs = (xf < 0 ? -xf : xf); | ||
| float sign = ((xbits.uintval() >> 15) == 1 ? -1.0 : 1.0); |
There was a problem hiding this comment.
Nit: unnecessary parentheses.
wldfngrs
force-pushed
the
add_asinf16_function
branch
from
5a1414f to
c6002ff
Compare
overmighty
reviewed
Jan 29, 2025
wldfngrs
force-pushed
the
add_asinf16_function
branch
from
de9fdff to
299a753
Compare
overmighty
requested changes
Feb 3, 2025
|
✅ With the latest revision this PR passed the C/C++ code formatter. |
overmighty
requested changes
Feb 4, 2025
libc/src/math/generic/asinf16.cpp
Outdated
Comment on lines
34
to
39
| float xsq = xf * xf; | ||
|
|
||
| // |x| <= 0x1p-1, |x| <= 0.5 | ||
| if (x_abs <= 0x3800) { | ||
| // asinf16(NaN) = NaN | ||
| if (xbits.is_nan()) { |
There was a problem hiding this comment.
This still computes xf * xf before checking that x is not NaN, and xbits.is_nan() is always false here since it implies x_abs > 0x7c00.
Suggested change
| float xsq = xf * xf; | |
| // |x| <= 0x1p-1, |x| <= 0.5 | |
| if (x_abs <= 0x3800) { | |
| // asinf16(NaN) = NaN | |
| if (xbits.is_nan()) { | |
| // |x| <= 0x1p-1, |x| <= 0.5, or x is NaN. | |
| if (LIBC_UNLIKELY(x_abs <= 0x3800 || x_abs > 0x7c00)) { | |
| // asinf16(NaN) = NaN | |
| if (xbits.is_nan()) { | |
| // ... | |
| } | |
| // ... | |
| } | |
| // ... | |
| float xsq = xf * xf; |
There was a problem hiding this comment.
What about directly checking xbits.is_nan(), like so?
...
float xf = x;
if (xbits.is_nan()) {
if (xbits.is_signaling_nan()) {
fputil::raise_except_if_required(FE_INVALID);
return FPBits::quiet_nan().get_val();
}
return x;
}
float xsq = xf * xf;
...
I mean, that'd be equivalent to checking if x_abs > 0x7c00, right?
Also, why include the x_abs <= 0x3800 in the NaN condition? If it falls within the domain of the function, it couldn't be NaN, right?
There was a problem hiding this comment.
What about directly checking
xbits.is_nan(), like so?
It might generate more instructions than x_abs > 0x7c00, although I haven't verified it in this case, and I haven't benchmarked it.
Also, why include the
x_abs <= 0x3800in the NaN condition? If it falls within the domain of the function, it couldn't be NaN, right?
The point is to group all unlikely cases under a single if (LIBC_UNLIKELY(...)) statement so that the compiler generates better code for the likely case: https://godbolt.org/z/bjzM5fv1z.
There was a problem hiding this comment.
It might generate more instructions than x_abs > 0x7c00, although I haven't verified it in this case, and I haven't benchmarked it.
How can I quickly verify the instructions generated?
There was a problem hiding this comment.
llvm-objdump -Cdr path/to/*.o
There was a problem hiding this comment.
// |x| > 0x1p0, |x| > 1, or x is NaN.
if (LIBC_UNLIKELY(x_abs > 0x3c00)) {
// asinf16(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();
}
It might generate more instructions than x_abs > 0x7c00, although I haven't verified it in this case, and I haven't benchmarked it
xbits.is_nan() and x_abs > 0x7c00 generate the same instructions. Also, what do you think of having the final check for inf values in the same conditional scope as the check for NaN-ness like above?
There was a problem hiding this comment.
Also, what do you think of having the final check for
infvalues in the same conditional scope as the check for NaN-ness like above?
It's what I would do. I would even move the x_abs <= 0x3800 case under the same if (LIBC_UNLIKELY(...)) statement, but I guess it's fine either way.
overmighty
approved these changes
Feb 10, 2025
overmighty
reviewed
Feb 10, 2025
Icohedron
pushed a commit
to Icohedron/llvm-project
that referenced
this pull request
Feb 11, 2025
Co-authored-by: OverMighty <[email protected]>
52 tasks
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.