[libc-commits] [libc] [libc][math][c23] Add asinf16() function (PR #124212)
via libc-commits
libc-commits at lists.llvm.org
Wed Feb 5 14:33:52 PST 2025
https://github.com/wldfngrs updated https://github.com/llvm/llvm-project/pull/124212
>From 9c6454fb3614d61a196de03a1eb7a3b06a9fa87e Mon Sep 17 00:00:00 2001
From: wldfngrs <wldfngrs at gmail.com>
Date: Fri, 24 Jan 2025 00:54:08 +0100
Subject: [PATCH 1/8] Add asinf16() function
---
libc/config/linux/x86_64/entrypoints.txt | 1 +
libc/docs/headers/math/index.rst | 2 +-
libc/include/math.yaml | 7 ++
libc/src/math/CMakeLists.txt | 1 +
libc/src/math/asinf16.h | 21 ++++
libc/src/math/generic/CMakeLists.txt | 19 ++++
libc/src/math/generic/asinf16.cpp | 120 ++++++++++++++++++++++
libc/test/src/math/CMakeLists.txt | 11 ++
libc/test/src/math/asinf16_test.cpp | 40 ++++++++
libc/test/src/math/smoke/CMakeLists.txt | 11 ++
libc/test/src/math/smoke/asinf16_test.cpp | 40 ++++++++
11 files changed, 272 insertions(+), 1 deletion(-)
create mode 100644 libc/src/math/asinf16.h
create mode 100644 libc/src/math/generic/asinf16.cpp
create mode 100644 libc/test/src/math/asinf16_test.cpp
create mode 100644 libc/test/src/math/smoke/asinf16_test.cpp
diff --git a/libc/config/linux/x86_64/entrypoints.txt b/libc/config/linux/x86_64/entrypoints.txt
index 366e4d34294d15..bdf010876910be 100644
--- a/libc/config/linux/x86_64/entrypoints.txt
+++ b/libc/config/linux/x86_64/entrypoints.txt
@@ -642,6 +642,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..0a0bdc0e28c3d1
--- /dev/null
+++ b/libc/src/math/generic/asinf16.cpp
@@ -0,0 +1,120 @@
+//===-- 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/cast.h"
+#include "src/__support/FPUtil/multiply_add.h"
+#include "src/__support/FPUtil/PolyEval.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, 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));
+}
+
+}
diff --git a/libc/test/src/math/CMakeLists.txt b/libc/test/src/math/CMakeLists.txt
index bbcdf2363c1e2b..7af20dcaab4bd7 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..43b3aeccd66364
--- /dev/null
+++ b/libc/test/src/math/asinf16_test.cpp
@@ -0,0 +1,40 @@
+//===-- 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 e4501eb75fa48a..43a580ac6b8942 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);
+}
>From d8fa4dddfd88f0bbac407410d77a0edd426a26fc Mon Sep 17 00:00:00 2001
From: wldfngrs <wldfngrs at gmail.com>
Date: Fri, 24 Jan 2025 00:55:47 +0100
Subject: [PATCH 2/8] clang-format
---
libc/src/math/generic/asinf16.cpp | 42 ++++++++++++++++-------------
libc/test/src/math/asinf16_test.cpp | 6 +++--
2 files changed, 28 insertions(+), 20 deletions(-)
diff --git a/libc/src/math/generic/asinf16.cpp b/libc/src/math/generic/asinf16.cpp
index 0a0bdc0e28c3d1..d340a7244cf34d 100644
--- a/libc/src/math/generic/asinf16.cpp
+++ b/libc/src/math/generic/asinf16.cpp
@@ -11,9 +11,9 @@
#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/PolyEval.h"
#include "src/__support/FPUtil/sqrt.h"
#include "src/__support/macros/optimization.h"
@@ -32,17 +32,18 @@ LLVM_LIBC_FUNCTION(float16, asinf16, (float16 x)) {
float xf = x;
float xsq = xf * xf;
-
+
// |x| <= 0x1p-1, |x| <= 0.5
- if (x_abs <= 0x3800) {
+ 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, then,
+ // 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
@@ -57,10 +58,12 @@ LLVM_LIBC_FUNCTION(float16, asinf16, (float16 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);
+ 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
@@ -72,7 +75,7 @@ LLVM_LIBC_FUNCTION(float16, asinf16, (float16 x)) {
return FPBits::quiet_nan().get_val();
}
- // When |x| > 0.5, assume that 0.5 < |x| <= 1,
+ // 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)
@@ -95,26 +98,29 @@ LLVM_LIBC_FUNCTION(float16, asinf16, (float16 x)) {
// 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.
-
+ //
+ // 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));
-}
+ 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/asinf16_test.cpp b/libc/test/src/math/asinf16_test.cpp
index 43b3aeccd66364..9593cad16ac772 100644
--- a/libc/test/src/math/asinf16_test.cpp
+++ b/libc/test/src/math/asinf16_test.cpp
@@ -27,7 +27,8 @@ 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);
+ EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Asin, x,
+ LIBC_NAMESPACE::asinf16(x), 0.5);
}
}
@@ -35,6 +36,7 @@ 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);
+ EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Asin, x,
+ LIBC_NAMESPACE::asinf16(x), 0.5);
}
}
>From c6002ffd3a652d8ba3cb1f04a0e6c1ed819f8921 Mon Sep 17 00:00:00 2001
From: wldfngrs <wldfngrs at gmail.com>
Date: Mon, 27 Jan 2025 17:23:48 +0100
Subject: [PATCH 3/8] nit
---
libc/src/math/generic/asinf16.cpp | 40 ++++++++++-------------
libc/test/src/math/smoke/asinf16_test.cpp | 3 --
2 files changed, 18 insertions(+), 25 deletions(-)
diff --git a/libc/src/math/generic/asinf16.cpp b/libc/src/math/generic/asinf16.cpp
index d340a7244cf34d..df1d7ec1cd1ed2 100644
--- a/libc/src/math/generic/asinf16.cpp
+++ b/libc/src/math/generic/asinf16.cpp
@@ -35,24 +35,23 @@ LLVM_LIBC_FUNCTION(float16, asinf16, (float16 x)) {
// |x| <= 0x1p-1, |x| <= 0.5
if (x_abs <= 0x3800) {
+ // asinf16(+/-0) = +/-0
+ if (LIBC_UNLIKELY(x_abs == 0))
+ return x;
+
+ // Exhaustive tests show that,
+ // for |x| <= 0x1.878p-9, 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
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)) {
+ (xbits.is_neg() && rounding == FE_DOWNWARD))
return fputil::cast<float16>(fputil::multiply_add(xf, 0x1.0p-11f, xf));
- }
return x;
}
@@ -66,11 +65,8 @@ LLVM_LIBC_FUNCTION(float16, asinf16, (float16 x)) {
// |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);
- }
+ fputil::set_errno_if_required(EDOM);
+ fputil::raise_except_if_required(FE_INVALID);
return FPBits::quiet_nan().get_val();
}
@@ -91,13 +87,13 @@ LLVM_LIBC_FUNCTION(float16, asinf16, (float16 x)) {
// 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
+ // asin(sqrt(u)) = z/2
+ // 2 * asin(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))
+ // asin(x) = pi/2 - 2 * asin(sqrt(u))
//
// WHY?
// 11: Recall [7], u = (1 - x)/2
@@ -109,7 +105,7 @@ LLVM_LIBC_FUNCTION(float16, asinf16, (float16 x)) {
// 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 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);
diff --git a/libc/test/src/math/smoke/asinf16_test.cpp b/libc/test/src/math/smoke/asinf16_test.cpp
index 303a3a302c2d5b..76ce36f673a8af 100644
--- a/libc/test/src/math/smoke/asinf16_test.cpp
+++ b/libc/test/src/math/smoke/asinf16_test.cpp
@@ -17,9 +17,6 @@ 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);
>From 299a753074882b5d0689802e1cf6d043e14afa85 Mon Sep 17 00:00:00 2001
From: wldfngrs <wldfngrs at gmail.com>
Date: Fri, 31 Jan 2025 00:19:34 +0100
Subject: [PATCH 4/8] fix behaviour for 'signaling NaN'/'quiet nan' inputs
---
libc/src/math/generic/asinf16.cpp | 11 ++++++++---
libc/test/src/math/smoke/asinf16_test.cpp | 2 ++
2 files changed, 10 insertions(+), 3 deletions(-)
diff --git a/libc/src/math/generic/asinf16.cpp b/libc/src/math/generic/asinf16.cpp
index df1d7ec1cd1ed2..56e5b0905d957c 100644
--- a/libc/src/math/generic/asinf16.cpp
+++ b/libc/src/math/generic/asinf16.cpp
@@ -65,8 +65,13 @@ LLVM_LIBC_FUNCTION(float16, asinf16, (float16 x)) {
// |x| > 1, asinf16(x) = NaN
if (LIBC_UNLIKELY(x_abs > 0x3c00)) {
- fputil::set_errno_if_required(EDOM);
- fputil::raise_except_if_required(FE_INVALID);
+ // |x| <= +/-inf
+ if (LIBC_UNLIKELY(x_abs <= 0x7c00)) {
+ if (!xbits.is_signaling_nan())
+ fputil::set_errno_if_required(EDOM);
+
+ fputil::raise_except_if_required(FE_INVALID);
+ }
return FPBits::quiet_nan().get_val();
}
@@ -91,7 +96,7 @@ LLVM_LIBC_FUNCTION(float16, asinf16, (float16 x)) {
// 2 * asin(sqrt(u)) = z
// 9: Recall [3], z = pi/2 - y. Therefore:
// y = pi/2 - z
- // y = pi/2 - 2 * arcsin(sqrt(u))
+ // y = pi/2 - 2 * asin(sqrt(u))
// 10: Recall [1], y = asin(x). Therefore:
// asin(x) = pi/2 - 2 * asin(sqrt(u))
//
diff --git a/libc/test/src/math/smoke/asinf16_test.cpp b/libc/test/src/math/smoke/asinf16_test.cpp
index 76ce36f673a8af..4b4836cc579581 100644
--- a/libc/test/src/math/smoke/asinf16_test.cpp
+++ b/libc/test/src/math/smoke/asinf16_test.cpp
@@ -16,6 +16,8 @@ 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);
>From dc173e35db3a16407c053cb444d2d4174c569c3f Mon Sep 17 00:00:00 2001
From: wldfngrs <wldfngrs at gmail.com>
Date: Mon, 3 Feb 2025 15:31:51 +0100
Subject: [PATCH 5/8] convert signaling NaN (sNaN) to quiet NaN (qNaN)
---
libc/src/math/generic/asinf16.cpp | 14 +++++++++++---
1 file changed, 11 insertions(+), 3 deletions(-)
diff --git a/libc/src/math/generic/asinf16.cpp b/libc/src/math/generic/asinf16.cpp
index 56e5b0905d957c..4c42624c2a6747 100644
--- a/libc/src/math/generic/asinf16.cpp
+++ b/libc/src/math/generic/asinf16.cpp
@@ -35,6 +35,16 @@ LLVM_LIBC_FUNCTION(float16, asinf16, (float16 x)) {
// |x| <= 0x1p-1, |x| <= 0.5
if (x_abs <= 0x3800) {
+ //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;
+ }
+
// asinf16(+/-0) = +/-0
if (LIBC_UNLIKELY(x_abs == 0))
return x;
@@ -67,9 +77,7 @@ LLVM_LIBC_FUNCTION(float16, asinf16, (float16 x)) {
if (LIBC_UNLIKELY(x_abs > 0x3c00)) {
// |x| <= +/-inf
if (LIBC_UNLIKELY(x_abs <= 0x7c00)) {
- if (!xbits.is_signaling_nan())
- fputil::set_errno_if_required(EDOM);
-
+ fputil::set_errno_if_required(EDOM);
fputil::raise_except_if_required(FE_INVALID);
}
>From 96ab6365395de2e7cf0954f30af2be24a5cc973a Mon Sep 17 00:00:00 2001
From: wldfngrs <wldfngrs at gmail.com>
Date: Mon, 3 Feb 2025 15:40:31 +0100
Subject: [PATCH 6/8] clang-format
---
libc/src/math/generic/asinf16.cpp | 6 +++---
1 file changed, 3 insertions(+), 3 deletions(-)
diff --git a/libc/src/math/generic/asinf16.cpp b/libc/src/math/generic/asinf16.cpp
index 4c42624c2a6747..f359f375e35963 100644
--- a/libc/src/math/generic/asinf16.cpp
+++ b/libc/src/math/generic/asinf16.cpp
@@ -35,16 +35,16 @@ LLVM_LIBC_FUNCTION(float16, asinf16, (float16 x)) {
// |x| <= 0x1p-1, |x| <= 0.5
if (x_abs <= 0x3800) {
- //asinf16(NaN) = NaN
+ // 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 FPBits::quiet_nan().get_val();
}
return x;
}
-
+
// asinf16(+/-0) = +/-0
if (LIBC_UNLIKELY(x_abs == 0))
return x;
>From 5a5136cbffb67936778ba12891c969a4ba735c10 Mon Sep 17 00:00:00 2001
From: wldfngrs <wldfngrs at gmail.com>
Date: Wed, 5 Feb 2025 19:05:32 +0100
Subject: [PATCH 7/8] minor changes
---
libc/src/math/generic/asinf16.cpp | 28 +++++++++++------------
libc/test/src/math/smoke/asinf16_test.cpp | 3 +++
2 files changed, 16 insertions(+), 15 deletions(-)
diff --git a/libc/src/math/generic/asinf16.cpp b/libc/src/math/generic/asinf16.cpp
index f359f375e35963..bbcd133863e3e4 100644
--- a/libc/src/math/generic/asinf16.cpp
+++ b/libc/src/math/generic/asinf16.cpp
@@ -31,10 +31,8 @@ LLVM_LIBC_FUNCTION(float16, asinf16, (float16 x)) {
uint16_t x_abs = x_u & 0x7fff;
float xf = x;
- float xsq = xf * xf;
-
- // |x| <= 0x1p-1, |x| <= 0.5
- if (x_abs <= 0x3800) {
+ // |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()) {
@@ -44,7 +42,18 @@ LLVM_LIBC_FUNCTION(float16, asinf16, (float16 x)) {
return x;
}
+
+ // 1 < |x| <= +/-inf
+ fputil::raise_except_if_required(FE_INVALID);
+ fputil::set_errno_if_required(EDOM);
+
+ return FPBits::quiet_nan().get_val();
+ }
+
+ float xsq = xf * xf;
+ // |x| <= 0x1p-1, |x| <= 0.5
+ if (x_abs <= 0x3800) {
// asinf16(+/-0) = +/-0
if (LIBC_UNLIKELY(x_abs == 0))
return x;
@@ -73,17 +82,6 @@ LLVM_LIBC_FUNCTION(float16, asinf16, (float16 x)) {
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:
diff --git a/libc/test/src/math/smoke/asinf16_test.cpp b/libc/test/src/math/smoke/asinf16_test.cpp
index 4b4836cc579581..9f675b08319c00 100644
--- a/libc/test/src/math/smoke/asinf16_test.cpp
+++ b/libc/test/src/math/smoke/asinf16_test.cpp
@@ -19,6 +19,9 @@ TEST_F(LlvmLibcAsinf16Test, SpecialNumbers) {
EXPECT_FP_EQ(aNaN, LIBC_NAMESPACE::asinf16(aNaN));
EXPECT_MATH_ERRNO(0);
+ EXPECT_FP_EQ_WITH_EXCEPTION(aNaN, LIBC_NAMESPACE::asinf16(sNaN), FE_INVALID);
+ EXPECT_MATH_ERRNO(0);
+
EXPECT_FP_EQ(zero, LIBC_NAMESPACE::asinf16(zero));
EXPECT_MATH_ERRNO(0);
>From cdf231a12979d68dc1933cc347ca42fe095fc9aa Mon Sep 17 00:00:00 2001
From: wldfngrs <wldfngrs at gmail.com>
Date: Wed, 5 Feb 2025 23:33:14 +0100
Subject: [PATCH 8/8] clang formatt
---
libc/src/math/generic/asinf16.cpp | 6 +++---
1 file changed, 3 insertions(+), 3 deletions(-)
diff --git a/libc/src/math/generic/asinf16.cpp b/libc/src/math/generic/asinf16.cpp
index bbcd133863e3e4..518c384a615304 100644
--- a/libc/src/math/generic/asinf16.cpp
+++ b/libc/src/math/generic/asinf16.cpp
@@ -42,14 +42,14 @@ LLVM_LIBC_FUNCTION(float16, asinf16, (float16 x)) {
return x;
}
-
+
// 1 < |x| <= +/-inf
fputil::raise_except_if_required(FE_INVALID);
fputil::set_errno_if_required(EDOM);
-
+
return FPBits::quiet_nan().get_val();
}
-
+
float xsq = xf * xf;
// |x| <= 0x1p-1, |x| <= 0.5
More information about the libc-commits
mailing list