[libc-commits] [libc] e2f065c - [libc][math] Implement asinf function correctly rounded for all rounding modes.
Tue Ly via libc-commits
libc-commits at lists.llvm.org
Wed Sep 7 16:28:03 PDT 2022
Author: Tue Ly
Date: 2022-09-07T19:27:47-04:00
New Revision: e2f065c2a3538c85ad7f6964b696dcdc1c3f95e4
URL: https://github.com/llvm/llvm-project/commit/e2f065c2a3538c85ad7f6964b696dcdc1c3f95e4
DIFF: https://github.com/llvm/llvm-project/commit/e2f065c2a3538c85ad7f6964b696dcdc1c3f95e4.diff
LOG: [libc][math] Implement asinf function correctly rounded for all rounding modes.
Implement asinf function correctly rounded for all rounding modes.
For `|x| <= 0.5`, we approximate `asin(x)` by
```
asin(x) = x * P(x^2)
```
where `P(X^2) = Q(X)` is a degree-20 minimax even polynomial approximating
`asin(x)/x` on `[0, 0.5]` generated by Sollya with:
```
> Q = fpminimax(asin(x)/x, [|0, 2, 4, 6, 8, 10, 12, 14, 16, 18, 20|],
[|1, D...|], [0, 0.5]);
```
When `|x| > 0.5`, we perform range reduction as follow:
Assume further that `0.5 < x <= 1`, and let:
```
y = asin(x)
```
We will use the double angle formula:
```
cos(2X) = 1 - 2 sin^2(X)
```
and the complement angle identity:
```
x = sin(y) = cos(pi/2 - y)
= 1 - 2 sin^2 (pi/4 - y/2)
```
So:
```
sin(pi/4 - y/2) = sqrt( (1 - x)/2 )
```
And hence:
```
pi/4 - y/2 = asin( sqrt( (1 - x)/2 ) )
```
Equivalently:
```
asin(x) = y = pi/2 - 2 * asin( sqrt( (1 - x)/2 ) )
```
Let `u = (1 - x)/2`, then
```
asin(x) = pi/2 - 2 * asin(u)
```
Moreover, since `0.5 < x <= 1`,
```
0 <= u < 1/4, and 0 <= sqrt(u) < 0.5.
```
And hence we can reuse the same polynomial approximation of `asin(x)` when
`|x| <= 0.5`:
```
asin(x) = pi/2 - 2 * u * P(u^2).
```
Performance benchmark using `perf` tool from the CORE-MATH project on Ryzen 1700:
```
$ CORE_MATH_PERF_MODE="rdtsc" ./perf.sh asinf
CORE-MATH reciprocal throughput : 23.418
System LIBC reciprocal throughput : 27.310
LIBC reciprocal throughput : 22.741
$ CORE_MATH_PERF_MODE="rdtsc" ./perf.sh asinf --latency
GNU libc version: 2.35
GNU libc release: stable
CORE-MATH latency : 58.884
System LIBC latency : 62.055
LIBC latency : 62.037
```
Reviewed By: orex, zimmermann6
Differential Revision: https://reviews.llvm.org/D133400
Added:
libc/src/math/asinf.h
libc/src/math/generic/asinf.cpp
libc/test/src/math/asinf_test.cpp
libc/test/src/math/exhaustive/asinf_test.cpp
Modified:
libc/config/darwin/arm/entrypoints.txt
libc/config/linux/aarch64/entrypoints.txt
libc/config/linux/x86_64/entrypoints.txt
libc/config/windows/entrypoints.txt
libc/docs/math.rst
libc/spec/stdc.td
libc/src/math/CMakeLists.txt
libc/src/math/generic/CMakeLists.txt
libc/test/src/math/CMakeLists.txt
libc/test/src/math/exhaustive/CMakeLists.txt
libc/utils/MPFRWrapper/MPFRUtils.cpp
libc/utils/MPFRWrapper/MPFRUtils.h
Removed:
################################################################################
diff --git a/libc/config/darwin/arm/entrypoints.txt b/libc/config/darwin/arm/entrypoints.txt
index 0b5e21984a3f9..433c245e640b5 100644
--- a/libc/config/darwin/arm/entrypoints.txt
+++ b/libc/config/darwin/arm/entrypoints.txt
@@ -105,6 +105,7 @@ set(TARGET_LIBM_ENTRYPOINTS
libc.src.fenv.feupdateenv
# math.h entrypoints
+ libc.src.math.asinf
libc.src.math.atanf
libc.src.math.atanhf
libc.src.math.copysign
diff --git a/libc/config/linux/aarch64/entrypoints.txt b/libc/config/linux/aarch64/entrypoints.txt
index 0e8e8df8326d4..abc6d77397de3 100644
--- a/libc/config/linux/aarch64/entrypoints.txt
+++ b/libc/config/linux/aarch64/entrypoints.txt
@@ -149,6 +149,7 @@ set(TARGET_LIBM_ENTRYPOINTS
libc.src.fenv.feupdateenv
# math.h entrypoints
+ libc.src.math.asinf
libc.src.math.atanf
libc.src.math.atanhf
libc.src.math.copysign
diff --git a/libc/config/linux/x86_64/entrypoints.txt b/libc/config/linux/x86_64/entrypoints.txt
index 80e178aeaa1b5..60ec5f0748509 100644
--- a/libc/config/linux/x86_64/entrypoints.txt
+++ b/libc/config/linux/x86_64/entrypoints.txt
@@ -149,6 +149,7 @@ set(TARGET_LIBM_ENTRYPOINTS
libc.src.fenv.feupdateenv
# math.h entrypoints
+ libc.src.math.asinf
libc.src.math.atanf
libc.src.math.atanhf
libc.src.math.copysign
diff --git a/libc/config/windows/entrypoints.txt b/libc/config/windows/entrypoints.txt
index bf13a5c082c08..a6791813cab2d 100644
--- a/libc/config/windows/entrypoints.txt
+++ b/libc/config/windows/entrypoints.txt
@@ -106,6 +106,7 @@ set(TARGET_LIBM_ENTRYPOINTS
libc.src.fenv.feupdateenv
# math.h entrypoints
+ libc.src.math.asinf
libc.src.math.atanf
libc.src.math.atanhf
libc.src.math.copysign
diff --git a/libc/docs/math.rst b/libc/docs/math.rst
index 9356f023ed91b..0cbd53329776b 100644
--- a/libc/docs/math.rst
+++ b/libc/docs/math.rst
@@ -118,7 +118,7 @@ Higher Math Functions
============== ================ =============== ======================
acos
acosh
-asin
+asin |check|
asinh
atan |check|
atan2
@@ -154,6 +154,7 @@ Accuracy of Higher Math Functions
============== ================ =============== ======================
<Func> <Func_f> (float) <Func> (double) <Func_l> (long double)
============== ================ =============== ======================
+asin |check|
atan |check|
atanh |check|
cos |check| large
@@ -203,6 +204,8 @@ Performance
| +-----------+-------------------+-----------+-------------------+ +------------+-------------------------+--------------+---------------+
| | LLVM libc | Reference (glibc) | LLVM libc | Reference (glibc) | | CPU | OS | Compiler | Special flags |
+==============+===========+===================+===========+===================+=====================================+============+=========================+==============+===============+
+| asinf | 23 | 27 | 62 | 62 | :math:`[-1, 1]` | Ryzen 1700 | Ubuntu 22.04 LTS x86_64 | Clang 14.0.0 | FMA |
++--------------+-----------+-------------------+-----------+-------------------+-------------------------------------+------------+-------------------------+--------------+---------------+
| atanf | 27 | 29 | 79 | 68 | :math:`[-10, 10]` | Ryzen 1700 | Ubuntu 22.04 LTS x86_64 | Clang 14.0.0 | FMA |
+--------------+-----------+-------------------+-----------+-------------------+-------------------------------------+------------+-------------------------+--------------+---------------+
| atanhf | 20 | 66 | 71 | 133 | :math:`[-1, 1]` | Ryzen 1700 | Ubuntu 22.04 LTS x86_64 | Clang 14.0.0 | FMA |
diff --git a/libc/spec/stdc.td b/libc/spec/stdc.td
index d7daa1d26b05c..fdb6beb1cdf02 100644
--- a/libc/spec/stdc.td
+++ b/libc/spec/stdc.td
@@ -480,6 +480,7 @@ def StdC : StandardSpec<"stdc"> {
FunctionSpec<"sinhf", RetValSpec<FloatType>, [ArgSpec<FloatType>]>,
FunctionSpec<"tanhf", RetValSpec<FloatType>, [ArgSpec<FloatType>]>,
+ FunctionSpec<"asinf", RetValSpec<FloatType>, [ArgSpec<FloatType>]>,
FunctionSpec<"atanf", RetValSpec<FloatType>, [ArgSpec<FloatType>]>,
FunctionSpec<"atanhf", RetValSpec<FloatType>, [ArgSpec<FloatType>]>,
diff --git a/libc/src/math/CMakeLists.txt b/libc/src/math/CMakeLists.txt
index eeb78d5aa45c1..313ff6a226ec1 100644
--- a/libc/src/math/CMakeLists.txt
+++ b/libc/src/math/CMakeLists.txt
@@ -64,7 +64,9 @@ add_entrypoint_object(
-O3
)
+add_math_entrypoint_object(asinf)
add_math_entrypoint_object(atanf)
+
add_math_entrypoint_object(atanhf)
add_math_entrypoint_object(ceil)
diff --git a/libc/src/math/asinf.h b/libc/src/math/asinf.h
new file mode 100644
index 0000000000000..71821277ef761
--- /dev/null
+++ b/libc/src/math/asinf.h
@@ -0,0 +1,18 @@
+//===-- Implementation header for asinf -------------------------*- 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_ASINF_H
+#define LLVM_LIBC_SRC_MATH_ASINF_H
+
+namespace __llvm_libc {
+
+float asinf(float x);
+
+} // namespace __llvm_libc
+
+#endif // LLVM_LIBC_SRC_MATH_ASINF_H
diff --git a/libc/src/math/generic/CMakeLists.txt b/libc/src/math/generic/CMakeLists.txt
index 3436c504cf0b2..a946ae4437952 100644
--- a/libc/src/math/generic/CMakeLists.txt
+++ b/libc/src/math/generic/CMakeLists.txt
@@ -1288,6 +1288,22 @@ add_object_library(
libc.include.math
)
+add_entrypoint_object(
+ asinf
+ SRCS
+ asinf.cpp
+ HDRS
+ ../asinf.h
+ DEPENDS
+ libc.src.__support.FPUtil.except_value_utils
+ libc.src.__support.FPUtil.fp_bits
+ libc.src.__support.FPUtil.multiply_add
+ libc.src.__support.FPUtil.polyeval
+ libc.src.__support.FPUtil.sqrt
+ COMPILE_OPTIONS
+ -O3
+)
+
add_entrypoint_object(
atanf
SRCS
diff --git a/libc/src/math/generic/asinf.cpp b/libc/src/math/generic/asinf.cpp
new file mode 100644
index 0000000000000..262af0fcbd3bc
--- /dev/null
+++ b/libc/src/math/generic/asinf.cpp
@@ -0,0 +1,169 @@
+//===-- Single-precision asin 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/asinf.h"
+#include "src/__support/FPUtil/FEnvImpl.h"
+#include "src/__support/FPUtil/FPBits.h"
+#include "src/__support/FPUtil/PolyEval.h"
+#include "src/__support/FPUtil/except_value_utils.h"
+#include "src/__support/FPUtil/multiply_add.h"
+#include "src/__support/FPUtil/sqrt.h"
+
+#include <errno.h>
+
+namespace __llvm_libc {
+
+// PI / 2
+constexpr double M_MATH_PI_2 = 0x1.921fb54442d18p+0;
+
+static constexpr size_t N_EXCEPTS = 2;
+
+// Exceptional values when |x| <= 0.5
+static constexpr fputil::ExceptValues<float, N_EXCEPTS> ASINF_EXCEPTS_LO = {{
+ // (inputs, RZ output, RU offset, RD offset, RN offset)
+ // x = 0x1.137f0cp-5, asinf(x) = 0x1.138c58p-5 (RZ)
+ {0x3d09bf86, 0x3d09c62c, 1, 0, 1},
+ // x = 0x1.cbf43cp-4, asinf(x) = 0x1.cced1cp-4 (RZ)
+ {0x3de5fa1e, 0x3de6768e, 1, 0, 0},
+}};
+
+// Exceptional values when 0.5 < |x| <= 1
+static constexpr fputil::ExceptValues<float, N_EXCEPTS> ASINF_EXCEPTS_HI = {{
+ // (inputs, RZ output, RU offset, RD offset, RN offset)
+ // x = 0x1.107434p-1, asinf(x) = 0x1.1f4b64p-1 (RZ)
+ {0x3f083a1a, 0x3f0fa5b2, 1, 0, 0},
+ // x = 0x1.ee836cp-1, asinf(x) = 0x1.4f0654p0 (RZ)
+ {0x3f7741b6, 0x3fa7832a, 1, 0, 0},
+}};
+
+// > Q = fpminimax(asin(x)/x, [|0, 2, 4, 6, 8, 10, 12, 14, 16, 18, 20|],
+// [|1, D...|], [0, 0.5]);
+static constexpr double COEFFS[10] = {
+ 0x1.5555555540fa1p-3, 0x1.333333512edc2p-4, 0x1.6db6cc1541b31p-5,
+ 0x1.f1caff324770ep-6, 0x1.6e43899f5f4f4p-6, 0x1.1f847cf652577p-6,
+ 0x1.9b60f47f87146p-7, 0x1.259e2634c494fp-6, -0x1.df946fa875ddp-8,
+ 0x1.02311ecf99c28p-5};
+
+LLVM_LIBC_FUNCTION(float, asinf, (float x)) {
+ using FPBits = typename fputil::FPBits<float>;
+ FPBits xbits(x);
+ uint32_t x_uint = xbits.uintval();
+ uint32_t x_abs = xbits.uintval() & 0x7fff'ffffU;
+ constexpr double SIGN[2] = {1.0, -1.0};
+ uint32_t x_sign = x_uint >> 31;
+
+ // |x| <= 0.5-ish
+ if (x_abs < 0x3f04'471dU) {
+ // |x| < 0x1.d12edp-12
+ if (unlikely(x_abs < 0x39e8'9768U)) {
+ // When |x| < 2^-12, the relative error of the approximation asin(x) ~ x
+ // is:
+ // |asin(x) - x| / |asin(x)| < |x^3| / (6|x|)
+ // = x^2 / 6
+ // < 2^-25
+ // < epsilon(1)/2.
+ // So the correctly rounded values of asin(x) are:
+ // = x + sign(x)*eps(x) if rounding mode = FE_TOWARDZERO,
+ // or (rounding mode = FE_UPWARD and x is
+ // negative),
+ // = x otherwise.
+ // To simplify the rounding decision and make it more efficient, we use
+ // fma(x, 2^-25, x) instead.
+ // An exhaustive test shows that this formula work correctly for all
+ // rounding modes up to |x| < 0x1.d12edp-12.
+ // Note: to use the formula x + 2^-25*x to decide the correct rounding, we
+ // do need fma(x, 2^-25, x) to prevent underflow caused by 2^-25*x when
+ // |x| < 2^-125. For targets without FMA instructions, we simply use
+ // double for intermediate results as it is more efficient than using an
+ // emulated version of FMA.
+#if defined(LIBC_TARGET_HAS_FMA)
+ return fputil::multiply_add(x, 0x1.0p-25f, x);
+#else
+ double xd = static_cast<double>(x);
+ return static_cast<float>(fputil::multiply_add(xd, 0x1.0p-25, xd));
+#endif // LIBC_TARGET_HAS_FMA
+ }
+
+ // Check for exceptional values
+ if (auto r = ASINF_EXCEPTS_LO.lookup_odd(x_abs, x_sign);
+ unlikely(r.has_value()))
+ return r.value();
+
+ // For |x| <= 0.5, we approximate asinf(x) by:
+ // asin(x) = x * P(x^2)
+ // Where P(X^2) = Q(X) is a degree-20 minimax even polynomial approximating
+ // asin(x)/x on [0, 0.5] generated by Sollya with:
+ // > Q = fpminimax(asin(x)/x, [|0, 2, 4, 6, 8, 10, 12, 14, 16, 18, 20|],
+ // [|1, D...|], [0, 0.5]);
+ // An exhaustive test shows that this approximation works well up to a
+ // little more than 0.5.
+ double xd = static_cast<double>(x);
+ double xsq = xd * xd;
+ double x4 = xsq * xsq;
+ double r1 = fputil::polyeval(x4, COEFFS[0], COEFFS[2], COEFFS[4], COEFFS[6],
+ COEFFS[8]);
+ double r2 = fputil::polyeval(x4, COEFFS[1], COEFFS[3], COEFFS[5], COEFFS[7],
+ COEFFS[9]);
+ double r3 = fputil::multiply_add(xsq, r2, r1);
+ double x3 = xd * xsq;
+ return fputil::multiply_add(x3, r3, xd);
+ }
+
+ // |x| > 1, return NaNs.
+ if (unlikely(x_abs > 0x3f80'0000U)) {
+ if (x_abs <= 0x7f80'0000U) {
+ errno = EDOM;
+ fputil::set_except(FE_INVALID);
+ }
+ return x +
+ FPBits::build_nan(1 << (fputil::MantissaWidth<float>::VALUE - 1));
+ }
+
+ // Check for exceptional values
+ if (auto r = ASINF_EXCEPTS_HI.lookup_odd(x_abs, x_sign);
+ unlikely(r.has_value()))
+ return r.value();
+
+ // When |x| > 0.5, we perform range reduction as follow:
+ // Assume further that 0.5 < x <= 1, and let:
+ // y = asin(x)
+ // We will use the double angle formula:
+ // cos(2y) = 1 - 2 sin^2(y)
+ // and the complement angle identity:
+ // x = sin(y) = cos(pi/2 - y)
+ // = 1 - 2 sin^2 (pi/4 - y/2)
+ // So:
+ // sin(pi/4 - y/2) = sqrt( (1 - x)/2 )
+ // And hence:
+ // pi/4 - y/2 = asin( sqrt( (1 - x)/2 ) )
+ // Equivalently:
+ // asin(x) = y = pi/2 - 2 * asin( sqrt( (1 - x)/2 ) )
+ // Let u = (1 - x)/2, then
+ // asin(x) = pi/2 - 2 * asin(u)
+ // Moreover, since 0.5 < x <= 1,
+ // 0 <= u < 1/4, and 0 <= sqrt(u) < 0.5,
+ // And hence we can reuse the same polynomial approximation of asin(x) when
+ // |x| <= 0.5:
+ // asin(x) = pi/2 - 2 * u * P(u^2),
+
+ xbits.set_sign(false);
+ double xd = static_cast<double>(xbits.get_val());
+ double u = fputil::multiply_add(-0.5, xd, 0.5);
+ double cv = -2 * fputil::sqrt(u);
+
+ double usq = u * u;
+ double r1 = fputil::polyeval(usq, COEFFS[0], COEFFS[2], COEFFS[4], COEFFS[6],
+ COEFFS[8]);
+ double r2 = fputil::polyeval(usq, COEFFS[1], COEFFS[3], COEFFS[5], COEFFS[7],
+ COEFFS[9]);
+ double r3 = fputil::multiply_add(u, r2, r1);
+ double r = fputil::multiply_add(cv * u, r3, M_MATH_PI_2 + cv);
+ return SIGN[x_sign] * r;
+}
+
+} // namespace __llvm_libc
diff --git a/libc/test/src/math/CMakeLists.txt b/libc/test/src/math/CMakeLists.txt
index e69d855c15fa5..b5669057c0509 100644
--- a/libc/test/src/math/CMakeLists.txt
+++ b/libc/test/src/math/CMakeLists.txt
@@ -1437,6 +1437,20 @@ add_fp_unittest(
libc.src.__support.FPUtil.fp_bits
)
+add_fp_unittest(
+ asinf_test
+ NEED_MPFR
+ SUITE
+ libc_math_unittests
+ SRCS
+ asinf_test.cpp
+ DEPENDS
+ libc.include.errno
+ libc.src.errno.errno
+ libc.src.math.asinf
+ libc.src.__support.FPUtil.fp_bits
+)
+
add_fp_unittest(
atanf_test
NEED_MPFR
diff --git a/libc/test/src/math/asinf_test.cpp b/libc/test/src/math/asinf_test.cpp
new file mode 100644
index 0000000000000..93443fc517bce
--- /dev/null
+++ b/libc/test/src/math/asinf_test.cpp
@@ -0,0 +1,80 @@
+//===-- Unittests for asinf
+//------------------------------------------------===//
+//
+// 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/FPUtil/FPBits.h"
+#include "src/math/asinf.h"
+#include "utils/MPFRWrapper/MPFRUtils.h"
+#include "utils/UnitTest/FPMatcher.h"
+#include "utils/UnitTest/Test.h"
+#include <math.h>
+
+#include <errno.h>
+#include <stdint.h>
+
+using FPBits = __llvm_libc::fputil::FPBits<float>;
+
+namespace mpfr = __llvm_libc::testing::mpfr;
+
+DECLARE_SPECIAL_CONSTANTS(float)
+
+TEST(LlvmLibcAsinfTest, SpecialNumbers) {
+ errno = 0;
+
+ EXPECT_FP_EQ(aNaN, __llvm_libc::asinf(aNaN));
+ EXPECT_MATH_ERRNO(0);
+
+ EXPECT_FP_EQ(0.0f, __llvm_libc::asinf(0.0f));
+ EXPECT_MATH_ERRNO(0);
+
+ EXPECT_FP_EQ(-0.0f, __llvm_libc::asinf(-0.0f));
+ EXPECT_MATH_ERRNO(0);
+
+ EXPECT_FP_EQ(aNaN, __llvm_libc::asinf(inf));
+ EXPECT_MATH_ERRNO(EDOM);
+
+ EXPECT_FP_EQ(aNaN, __llvm_libc::asinf(neg_inf));
+ EXPECT_MATH_ERRNO(EDOM);
+}
+
+TEST(LlvmLibcAsinfTest, InFloatRange) {
+ constexpr uint32_t COUNT = 1000000;
+ constexpr uint32_t STEP = UINT32_MAX / COUNT;
+ for (uint32_t i = 0, v = 0; i <= COUNT; ++i, v += STEP) {
+ float x = float(FPBits(v));
+ if (isnan(x) || isinf(x))
+ continue;
+ ASSERT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Asin, x,
+ __llvm_libc::asinf(x), 0.5);
+ }
+}
+
+TEST(LlvmLibcAsinfTest, SpecificBitPatterns) {
+ constexpr int N = 11;
+ constexpr uint32_t INPUTS[N] = {
+ 0x3f000000, // x = 0.5f
+ 0x3f3504f3, // x = sqrt(2)/2, FE_DOWNWARD
+ 0x3f3504f4, // x = sqrt(2)/2, FE_UPWARD
+ 0x3f5db3d7, // x = sqrt(3)/2, FE_DOWNWARD
+ 0x3f5db3d8, // x = sqrt(3)/2, FE_UPWARD
+ 0x3f800000, // x = 1.0f
+ 0x40000000, // x = 2.0f
+ 0x3d09bf86, // x = 0x1.137f0cp-5f
+ 0x3de5fa1e, // x = 0x1.cbf43cp-4f
+ 0x3f083a1a, // x = 0x1.107434p-1f
+ 0x3f7741b6, // x = 0x1.ee836cp-1f
+ };
+
+ for (int i = 0; i < N; ++i) {
+ float x = float(FPBits(INPUTS[i]));
+ EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Asin, x,
+ __llvm_libc::asinf(x), 0.5);
+ EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Asin, -x,
+ __llvm_libc::asinf(-x), 0.5);
+ }
+}
diff --git a/libc/test/src/math/exhaustive/CMakeLists.txt b/libc/test/src/math/exhaustive/CMakeLists.txt
index 5a45e3febe441..1051e87c40192 100644
--- a/libc/test/src/math/exhaustive/CMakeLists.txt
+++ b/libc/test/src/math/exhaustive/CMakeLists.txt
@@ -325,3 +325,19 @@ add_fp_unittest(
-lpthread
)
+add_fp_unittest(
+ asinf_test
+ NO_RUN_POSTBUILD
+ NEED_MPFR
+ SUITE
+ libc_math_exhaustive_tests
+ SRCS
+ asinf_test.cpp
+ DEPENDS
+ .exhaustive_test
+ libc.include.math
+ libc.src.math.asinf
+ libc.src.__support.FPUtil.fp_bits
+ LINK_LIBRARIES
+ -lpthread
+)
diff --git a/libc/test/src/math/exhaustive/asinf_test.cpp b/libc/test/src/math/exhaustive/asinf_test.cpp
new file mode 100644
index 0000000000000..fa16acbe27eb2
--- /dev/null
+++ b/libc/test/src/math/exhaustive/asinf_test.cpp
@@ -0,0 +1,78 @@
+//===-- Exhaustive test for asinf -----------------------------------------===//
+//
+// 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 "exhaustive_test.h"
+#include "src/__support/FPUtil/FPBits.h"
+#include "src/math/asinf.h"
+#include "utils/MPFRWrapper/MPFRUtils.h"
+
+#include <thread>
+
+using FPBits = __llvm_libc::fputil::FPBits<float>;
+
+namespace mpfr = __llvm_libc::testing::mpfr;
+
+struct LlvmLibcAsinfExhaustiveTest : public LlvmLibcExhaustiveTest<uint32_t> {
+ bool check(uint32_t start, uint32_t stop,
+ mpfr::RoundingMode rounding) override {
+ mpfr::ForceRoundingMode r(rounding);
+ uint32_t bits = start;
+ bool result = true;
+ do {
+ FPBits xbits(bits);
+ float x = float(xbits);
+ result &= EXPECT_MPFR_MATCH(mpfr::Operation::Asin, x,
+ __llvm_libc::asinf(x), 0.5, rounding);
+ // if (!result) break;
+ } while (bits++ < stop);
+ return result;
+ }
+};
+
+static const int NUM_THREADS = std::thread::hardware_concurrency();
+
+// Range: [0, Inf];
+static const uint32_t POS_START = 0x0000'0000U;
+static const uint32_t POS_STOP = 0x7f80'0000U;
+/
+
+ TEST_F(LlvmLibcAsinfExhaustiveTest, PostiveRangeRoundNearestTieToEven) {
+ test_full_range(POS_START, POS_STOP, mpfr::RoundingMode::Nearest);
+}
+
+TEST_F(LlvmLibcAsinfExhaustiveTest, PostiveRangeRoundUp) {
+ test_full_range(POS_START, POS_STOP, mpfr::RoundingMode::Upward);
+}
+
+TEST_F(LlvmLibcAsinfExhaustiveTest, PostiveRangeRoundDown) {
+ test_full_range(POS_START, POS_STOP, mpfr::RoundingMode::Downward);
+}
+
+TEST_F(LlvmLibcAsinfExhaustiveTest, PostiveRangeRoundTowardZero) {
+ test_full_range(POS_START, POS_STOP, mpfr::RoundingMode::TowardZero);
+}
+
+// Range: [-Inf, 0];
+static const uint32_t NEG_START = 0x8000'0000U;
+static const uint32_t NEG_STOP = 0xff80'0000U;
+
+TEST_F(LlvmLibcAsinfExhaustiveTest, NegativeRangeRoundNearestTieToEven) {
+ test_full_range(NEG_START, NEG_STOP, mpfr::RoundingMode::Nearest);
+}
+
+TEST_F(LlvmLibcAsinfExhaustiveTest, NegativeRangeRoundUp) {
+ test_full_range(NEG_START, NEG_STOP, mpfr::RoundingMode::Upward);
+}
+
+TEST_F(LlvmLibcAsinfExhaustiveTest, NegativeRangeRoundDown) {
+ test_full_range(NEG_START, NEG_STOP, mpfr::RoundingMode::Downward);
+}
+
+TEST_F(LlvmLibcAsinfExhaustiveTest, NegativeRangeRoundTowardZero) {
+ test_full_range(NEG_START, NEG_STOP, mpfr::RoundingMode::TowardZero);
+}
diff --git a/libc/utils/MPFRWrapper/MPFRUtils.cpp b/libc/utils/MPFRWrapper/MPFRUtils.cpp
index 22d056b699526..910233e87dda7 100644
--- a/libc/utils/MPFRWrapper/MPFRUtils.cpp
+++ b/libc/utils/MPFRWrapper/MPFRUtils.cpp
@@ -176,12 +176,20 @@ class MPFRNumber {
return *this;
}
+ bool is_nan() const { return mpfr_nan_p(value); }
+
MPFRNumber abs() const {
MPFRNumber result(*this);
mpfr_abs(result.value, value, mpfr_rounding);
return result;
}
+ MPFRNumber asin() const {
+ MPFRNumber result(*this);
+ mpfr_asin(result.value, value, mpfr_rounding);
+ return result;
+ }
+
MPFRNumber atan() const {
MPFRNumber result(*this);
mpfr_atan(result.value, value, mpfr_rounding);
@@ -436,6 +444,12 @@ class MPFRNumber {
if (thisAsT == input)
return T(0.0);
+ if (is_nan()) {
+ if (fputil::FPBits<T>(input).is_nan())
+ return T(0.0);
+ return T(fputil::FPBits<T>::inf());
+ }
+
int thisExponent = fputil::FPBits<T>(thisAsT).get_exponent();
int inputExponent = fputil::FPBits<T>(input).get_exponent();
// Adjust the exponents for denormal numbers.
@@ -512,6 +526,8 @@ unary_operation(Operation op, InputType input, unsigned int precision,
switch (op) {
case Operation::Abs:
return mpfrInput.abs();
+ case Operation::Asin:
+ return mpfrInput.asin();
case Operation::Atan:
return mpfrInput.atan();
case Operation::Atanh:
diff --git a/libc/utils/MPFRWrapper/MPFRUtils.h b/libc/utils/MPFRWrapper/MPFRUtils.h
index f098c29138877..2f9f6f03d6b30 100644
--- a/libc/utils/MPFRWrapper/MPFRUtils.h
+++ b/libc/utils/MPFRWrapper/MPFRUtils.h
@@ -25,6 +25,7 @@ enum class Operation : int {
// and output floating point numbers are of the same kind.
BeginUnaryOperationsSingleOutput,
Abs,
+ Asin,
Atan,
Atanh,
Ceil,
More information about the libc-commits
mailing list