[libc-commits] [libc] 78cc822 - [libc][math] Implement double precision acos correctly rounded for all rounding modes. (#138308)
via libc-commits
libc-commits at lists.llvm.org
Thu May 8 20:23:12 PDT 2025
Author: lntue
Date: 2025-05-08T23:23:09-04:00
New Revision: 78cc822aa6f5af0eda55089d22ba915b6d8e0216
URL: https://github.com/llvm/llvm-project/commit/78cc822aa6f5af0eda55089d22ba915b6d8e0216
DIFF: https://github.com/llvm/llvm-project/commit/78cc822aa6f5af0eda55089d22ba915b6d8e0216.diff
LOG: [libc][math] Implement double precision acos correctly rounded for all rounding modes. (#138308)
We reduce computation of `acos` to `asin` as follow:
When `|x| < 0.5`:
```math
acos(x) = \frac{\pi}{2} - asin(x).
```
For `0.5 <= |x| < 1`, let
```math
u = \frac{1 - \left| x \right|}{2},
```
then
```math
acos(x) = \begin{cases}
2 \cdot asin \left( \sqrt{u} \right) &, 0.5 \leq x < 1 \\
\pi - 2 \cdot asin \left( \sqrt{u} \right) &, -1 < x \leq 0.5
\end{cases}
```
Added:
libc/src/math/generic/acos.cpp
libc/test/src/math/acos_test.cpp
libc/test/src/math/smoke/acos_test.cpp
Modified:
libc/config/darwin/arm/entrypoints.txt
libc/config/linux/aarch64/entrypoints.txt
libc/config/linux/arm/entrypoints.txt
libc/config/linux/riscv/entrypoints.txt
libc/config/linux/x86_64/entrypoints.txt
libc/config/windows/entrypoints.txt
libc/docs/headers/math/index.rst
libc/include/math.yaml
libc/src/math/generic/CMakeLists.txt
libc/src/math/generic/asin.cpp
libc/src/math/generic/asin_utils.h
libc/test/src/math/CMakeLists.txt
libc/test/src/math/smoke/CMakeLists.txt
Removed:
################################################################################
diff --git a/libc/config/darwin/arm/entrypoints.txt b/libc/config/darwin/arm/entrypoints.txt
index 70c888aec064c..308fc49d681d7 100644
--- a/libc/config/darwin/arm/entrypoints.txt
+++ b/libc/config/darwin/arm/entrypoints.txt
@@ -135,6 +135,7 @@ set(TARGET_LIBM_ENTRYPOINTS
libc.src.fenv.feupdateenv
# math.h entrypoints
+ libc.src.math.acos
libc.src.math.acosf
libc.src.math.acoshf
libc.src.math.asin
diff --git a/libc/config/linux/aarch64/entrypoints.txt b/libc/config/linux/aarch64/entrypoints.txt
index 66d7576ffb882..520046f768b5d 100644
--- a/libc/config/linux/aarch64/entrypoints.txt
+++ b/libc/config/linux/aarch64/entrypoints.txt
@@ -410,6 +410,7 @@ set(TARGET_LIBM_ENTRYPOINTS
libc.src.fenv.feupdateenv
# math.h entrypoints
+ libc.src.math.acos
libc.src.math.acosf
libc.src.math.acoshf
libc.src.math.asin
diff --git a/libc/config/linux/arm/entrypoints.txt b/libc/config/linux/arm/entrypoints.txt
index b5e2f59d25a54..7432a7e912e81 100644
--- a/libc/config/linux/arm/entrypoints.txt
+++ b/libc/config/linux/arm/entrypoints.txt
@@ -242,6 +242,7 @@ set(TARGET_LIBM_ENTRYPOINTS
libc.src.fenv.feupdateenv
# math.h entrypoints
+ libc.src.math.acos
libc.src.math.acosf
libc.src.math.acoshf
libc.src.math.asin
diff --git a/libc/config/linux/riscv/entrypoints.txt b/libc/config/linux/riscv/entrypoints.txt
index 30d9d00dfefc9..c2a31b9f5c964 100644
--- a/libc/config/linux/riscv/entrypoints.txt
+++ b/libc/config/linux/riscv/entrypoints.txt
@@ -416,6 +416,7 @@ set(TARGET_LIBM_ENTRYPOINTS
libc.src.fenv.feupdateenv
# math.h entrypoints
+ libc.src.math.acos
libc.src.math.acosf
libc.src.math.acoshf
libc.src.math.asin
diff --git a/libc/config/linux/x86_64/entrypoints.txt b/libc/config/linux/x86_64/entrypoints.txt
index 2ac016f065b2e..9f447dd0d35d2 100644
--- a/libc/config/linux/x86_64/entrypoints.txt
+++ b/libc/config/linux/x86_64/entrypoints.txt
@@ -415,6 +415,7 @@ set(TARGET_LIBM_ENTRYPOINTS
libc.src.fenv.feupdateenv
# math.h entrypoints
+ libc.src.math.acos
libc.src.math.acosf
libc.src.math.acoshf
libc.src.math.asin
diff --git a/libc/config/windows/entrypoints.txt b/libc/config/windows/entrypoints.txt
index 37fa888d6498a..09021a08cf731 100644
--- a/libc/config/windows/entrypoints.txt
+++ b/libc/config/windows/entrypoints.txt
@@ -127,6 +127,7 @@ set(TARGET_LIBM_ENTRYPOINTS
libc.src.fenv.feupdateenv
# math.h entrypoints
+ libc.src.math.acos
libc.src.math.acosf
libc.src.math.acoshf
libc.src.math.asin
diff --git a/libc/docs/headers/math/index.rst b/libc/docs/headers/math/index.rst
index 9d4d6698b6122..6b0365f481a4c 100644
--- a/libc/docs/headers/math/index.rst
+++ b/libc/docs/headers/math/index.rst
@@ -249,7 +249,7 @@ Higher Math Functions
+-----------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+----------------------------+
| <Func> | <Func_f> (float) | <Func> (double) | <Func_l> (long double) | <Func_f16> (float16) | <Func_f128> (float128) | C23 Definition Section | C23 Error Handling Section |
+===========+==================+=================+========================+======================+========================+========================+============================+
-| acos | |check| | | | |check| | | 7.12.4.1 | F.10.1.1 |
+| acos | |check| | |check| | | |check| | | 7.12.4.1 | F.10.1.1 |
+-----------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+----------------------------+
| acosh | |check| | | | |check| | | 7.12.5.1 | F.10.2.1 |
+-----------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+----------------------------+
diff --git a/libc/include/math.yaml b/libc/include/math.yaml
index f8c57f7d48ab2..0a69653d156af 100644
--- a/libc/include/math.yaml
+++ b/libc/include/math.yaml
@@ -8,6 +8,12 @@ types:
enums: []
objects: []
functions:
+ - name: acos
+ standards:
+ - stdc
+ return_type: double
+ arguments:
+ - type: double
- name: acosf
standards:
- stdc
diff --git a/libc/src/math/generic/CMakeLists.txt b/libc/src/math/generic/CMakeLists.txt
index b08486ac61366..e5d1f9075bcb1 100644
--- a/libc/src/math/generic/CMakeLists.txt
+++ b/libc/src/math/generic/CMakeLists.txt
@@ -4117,6 +4117,7 @@ add_entrypoint_object(
HDRS
../asin.h
DEPENDS
+ .asin_utils
libc.src.__support.FPUtil.double_double
libc.src.__support.FPUtil.dyadic_float
libc.src.__support.FPUtil.fenv_impl
@@ -4164,6 +4165,26 @@ add_entrypoint_object(
libc.src.__support.macros.properties.types
)
+add_entrypoint_object(
+ acos
+ SRCS
+ acos.cpp
+ HDRS
+ ../acos.h
+ DEPENDS
+ .asin_utils
+ libc.src.__support.FPUtil.double_double
+ libc.src.__support.FPUtil.dyadic_float
+ 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
+ libc.src.__support.macros.properties.cpu_features
+)
+
add_entrypoint_object(
acospif16
SRCS
diff --git a/libc/src/math/generic/acos.cpp b/libc/src/math/generic/acos.cpp
new file mode 100644
index 0000000000000..c14721faef3ce
--- /dev/null
+++ b/libc/src/math/generic/acos.cpp
@@ -0,0 +1,278 @@
+//===-- Double-precision acos 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/acos.h"
+#include "asin_utils.h"
+#include "src/__support/FPUtil/FEnvImpl.h"
+#include "src/__support/FPUtil/FPBits.h"
+#include "src/__support/FPUtil/PolyEval.h"
+#include "src/__support/FPUtil/double_double.h"
+#include "src/__support/FPUtil/dyadic_float.h"
+#include "src/__support/FPUtil/multiply_add.h"
+#include "src/__support/FPUtil/sqrt.h"
+#include "src/__support/macros/config.h"
+#include "src/__support/macros/optimization.h" // LIBC_UNLIKELY
+#include "src/__support/macros/properties/cpu_features.h" // LIBC_TARGET_CPU_HAS_FMA
+
+namespace LIBC_NAMESPACE_DECL {
+
+using DoubleDouble = fputil::DoubleDouble;
+using Float128 = fputil::DyadicFloat<128>;
+
+LLVM_LIBC_FUNCTION(double, acos, (double x)) {
+ using FPBits = fputil::FPBits<double>;
+
+ FPBits xbits(x);
+ int x_exp = xbits.get_biased_exponent();
+
+ // |x| < 0.5.
+ if (x_exp < FPBits::EXP_BIAS - 1) {
+ // |x| < 2^-55.
+ if (LIBC_UNLIKELY(x_exp < FPBits::EXP_BIAS - 55)) {
+ // When |x| < 2^-55, acos(x) = pi/2
+#if defined(LIBC_MATH_HAS_SKIP_ACCURATE_PASS)
+ return PI_OVER_TWO.hi;
+#else
+ // Force the evaluation and prevent constant propagation so that it
+ // is rounded correctly for FE_UPWARD rounding mode.
+ return (xbits.abs().get_val() + 0x1.0p-160) + PI_OVER_TWO.hi;
+#endif // LIBC_MATH_HAS_SKIP_ACCURATE_PASS
+ }
+
+#ifdef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
+ // acos(x) = pi/2 - asin(x)
+ // = pi/2 - x * P(x^2)
+ double p = asin_eval(x * x);
+ return PI_OVER_TWO.hi + fputil::multiply_add(-x, p, PI_OVER_TWO.lo);
+#else
+ unsigned idx;
+ DoubleDouble x_sq = fputil::exact_mult(x, x);
+ double err = xbits.abs().get_val() * 0x1.0p-51;
+ // Polynomial approximation:
+ // p ~ asin(x)/x
+ DoubleDouble p = asin_eval(x_sq, idx, err);
+ // asin(x) ~ x * p
+ DoubleDouble r0 = fputil::exact_mult(x, p.hi);
+ // acos(x) = pi/2 - asin(x)
+ // ~ pi/2 - x * p
+ // = pi/2 - x * (p.hi + p.lo)
+ double r_hi = fputil::multiply_add(-x, p.hi, PI_OVER_TWO.hi);
+ // Use Dekker's 2SUM algorithm to compute the lower part.
+ double r_lo = ((PI_OVER_TWO.hi - r_hi) - r0.hi) - r0.lo;
+ r_lo = fputil::multiply_add(-x, p.lo, r_lo + PI_OVER_TWO.lo);
+
+ // Ziv's accuracy test.
+
+ double r_upper = r_hi + (r_lo + err);
+ double r_lower = r_hi + (r_lo - err);
+
+ if (LIBC_LIKELY(r_upper == r_lower))
+ return r_upper;
+
+ // Ziv's accuracy test failed, perform 128-bit calculation.
+
+ // Recalculate mod 1/64.
+ idx = static_cast<unsigned>(fputil::nearest_integer(x_sq.hi * 0x1.0p6));
+
+ // Get x^2 - idx/64 exactly. When FMA is available, double-double
+ // multiplication will be correct for all rounding modes. Otherwise we use
+ // Float128 directly.
+ Float128 x_f128(x);
+
+#ifdef LIBC_TARGET_CPU_HAS_FMA_DOUBLE
+ // u = x^2 - idx/64
+ Float128 u_hi(
+ fputil::multiply_add(static_cast<double>(idx), -0x1.0p-6, x_sq.hi));
+ Float128 u = fputil::quick_add(u_hi, Float128(x_sq.lo));
+#else
+ Float128 x_sq_f128 = fputil::quick_mul(x_f128, x_f128);
+ Float128 u = fputil::quick_add(
+ x_sq_f128, Float128(static_cast<double>(idx) * (-0x1.0p-6)));
+#endif // LIBC_TARGET_CPU_HAS_FMA_DOUBLE
+
+ Float128 p_f128 = asin_eval(u, idx);
+ // Flip the sign of x_f128 to perform subtraction.
+ x_f128.sign = x_f128.sign.negate();
+ Float128 r =
+ fputil::quick_add(PI_OVER_TWO_F128, fputil::quick_mul(x_f128, p_f128));
+
+ return static_cast<double>(r);
+#endif // LIBC_MATH_HAS_SKIP_ACCURATE_PASS
+ }
+ // |x| >= 0.5
+
+ double x_abs = xbits.abs().get_val();
+
+ // Maintaining the sign:
+ constexpr double SIGN[2] = {1.0, -1.0};
+ double x_sign = SIGN[xbits.is_neg()];
+ // |x| >= 1
+ if (LIBC_UNLIKELY(x_exp >= FPBits::EXP_BIAS)) {
+ // x = +-1, asin(x) = +- pi/2
+ if (x_abs == 1.0) {
+ // x = 1, acos(x) = 0,
+ // x = -1, acos(x) = pi
+ return x == 1.0 ? 0.0 : fputil::multiply_add(-x_sign, PI.hi, PI.lo);
+ }
+ // |x| > 1, return NaN.
+ if (xbits.is_quiet_nan())
+ return x;
+
+ // Set domain error for non-NaN input.
+ if (!xbits.is_nan())
+ fputil::set_errno_if_required(EDOM);
+
+ fputil::raise_except_if_required(FE_INVALID);
+ return FPBits::quiet_nan().get_val();
+ }
+
+ // When |x| >= 0.5, we perform range reduction as follow:
+ //
+ // When 0.5 <= x < 1, let:
+ // y = acos(x)
+ // We will use the double angle formula:
+ // cos(2y) = 1 - 2 sin^2(y)
+ // and the complement angle identity:
+ // x = cos(y) = 1 - 2 sin^2 (y/2)
+ // So:
+ // sin(y/2) = sqrt( (1 - x)/2 )
+ // And hence:
+ // y/2 = asin( sqrt( (1 - x)/2 ) )
+ // Equivalently:
+ // acos(x) = y = 2 * asin( sqrt( (1 - x)/2 ) )
+ // Let u = (1 - x)/2, then:
+ // acos(x) = 2 * asin( sqrt(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:
+ // acos(x) ~ 2 * sqrt(u) * P(u).
+ //
+ // When -1 < x <= -0.5, we reduce to the previous case using the formula:
+ // acos(x) = pi - acos(-x)
+ // = pi - 2 * asin ( sqrt( (1 + x)/2 ) )
+ // ~ pi - 2 * sqrt(u) * P(u),
+ // where u = (1 - |x|)/2.
+
+ // u = (1 - |x|)/2
+ double u = fputil::multiply_add(x_abs, -0.5, 0.5);
+ // v_hi + v_lo ~ sqrt(u).
+ // Let:
+ // h = u - v_hi^2 = (sqrt(u) - v_hi) * (sqrt(u) + v_hi)
+ // Then:
+ // sqrt(u) = v_hi + h / (sqrt(u) + v_hi)
+ // ~ v_hi + h / (2 * v_hi)
+ // So we can use:
+ // v_lo = h / (2 * v_hi).
+ double v_hi = fputil::sqrt<double>(u);
+
+#ifdef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
+ constexpr DoubleDouble CONST_TERM[2] = {{0.0, 0.0}, PI};
+ DoubleDouble const_term = CONST_TERM[xbits.is_neg()];
+
+ double p = asin_eval(u);
+ double scale = x_sign * 2.0 * v_hi;
+ double r = const_term.hi + fputil::multiply_add(scale, p, const_term.lo);
+ return r;
+#else
+
+#ifdef LIBC_TARGET_CPU_HAS_FMA_DOUBLE
+ double h = fputil::multiply_add(v_hi, -v_hi, u);
+#else
+ DoubleDouble v_hi_sq = fputil::exact_mult(v_hi, v_hi);
+ double h = (u - v_hi_sq.hi) - v_hi_sq.lo;
+#endif // LIBC_TARGET_CPU_HAS_FMA_DOUBLE
+
+ // Scale v_lo and v_hi by 2 from the formula:
+ // vh = v_hi * 2
+ // vl = 2*v_lo = h / v_hi.
+ double vh = v_hi * 2.0;
+ double vl = h / v_hi;
+
+ // Polynomial approximation:
+ // p ~ asin(sqrt(u))/sqrt(u)
+ unsigned idx;
+ double err = vh * 0x1.0p-51;
+
+ DoubleDouble p = asin_eval(DoubleDouble{0.0, u}, idx, err);
+
+ // Perform computations in double-double arithmetic:
+ // asin(x) = pi/2 - (v_hi + v_lo) * (ASIN_COEFFS[idx][0] + p)
+ DoubleDouble r0 = fputil::quick_mult(DoubleDouble{vl, vh}, p);
+
+ double r_hi, r_lo;
+ if (xbits.is_pos()) {
+ r_hi = r0.hi;
+ r_lo = r0.lo;
+ } else {
+ DoubleDouble r = fputil::exact_add(PI.hi, -r0.hi);
+ r_hi = r.hi;
+ r_lo = (PI.lo - r0.lo) + r.lo;
+ }
+
+ // Ziv's accuracy test.
+
+ double r_upper = r_hi + (r_lo + err);
+ double r_lower = r_hi + (r_lo - err);
+
+ if (LIBC_LIKELY(r_upper == r_lower))
+ return r_upper;
+
+ // Ziv's accuracy test failed, we redo the computations in Float128.
+ // Recalculate mod 1/64.
+ idx = static_cast<unsigned>(fputil::nearest_integer(u * 0x1.0p6));
+
+ // After the first step of Newton-Raphson approximating v = sqrt(u), we have
+ // that:
+ // sqrt(u) = v_hi + h / (sqrt(u) + v_hi)
+ // v_lo = h / (2 * v_hi)
+ // With error:
+ // sqrt(u) - (v_hi + v_lo) = h * ( 1/(sqrt(u) + v_hi) - 1/(2*v_hi) )
+ // = -h^2 / (2*v * (sqrt(u) + v)^2).
+ // Since:
+ // (sqrt(u) + v_hi)^2 ~ (2sqrt(u))^2 = 4u,
+ // we can add another correction term to (v_hi + v_lo) that is:
+ // v_ll = -h^2 / (2*v_hi * 4u)
+ // = -v_lo * (h / 4u)
+ // = -vl * (h / 8u),
+ // making the errors:
+ // sqrt(u) - (v_hi + v_lo + v_ll) = O(h^3)
+ // well beyond 128-bit precision needed.
+
+ // Get the rounding error of vl = 2 * v_lo ~ h / vh
+ // Get full product of vh * vl
+#ifdef LIBC_TARGET_CPU_HAS_FMA_DOUBLE
+ double vl_lo = fputil::multiply_add(-v_hi, vl, h) / v_hi;
+#else
+ DoubleDouble vh_vl = fputil::exact_mult(v_hi, vl);
+ double vl_lo = ((h - vh_vl.hi) - vh_vl.lo) / v_hi;
+#endif // LIBC_TARGET_CPU_HAS_FMA_DOUBLE
+ // vll = 2*v_ll = -vl * (h / (4u)).
+ double t = h * (-0.25) / u;
+ double vll = fputil::multiply_add(vl, t, vl_lo);
+ // m_v = -(v_hi + v_lo + v_ll).
+ Float128 m_v = fputil::quick_add(
+ Float128(vh), fputil::quick_add(Float128(vl), Float128(vll)));
+ m_v.sign = xbits.sign();
+
+ // Perform computations in Float128:
+ // acos(x) = (v_hi + v_lo + vll) * P(u) , when 0.5 <= x < 1,
+ // = pi - (v_hi + v_lo + vll) * P(u) , when -1 < x <= -0.5.
+ Float128 y_f128(fputil::multiply_add(static_cast<double>(idx), -0x1.0p-6, u));
+
+ Float128 p_f128 = asin_eval(y_f128, idx);
+ Float128 r_f128 = fputil::quick_mul(m_v, p_f128);
+
+ if (xbits.is_neg())
+ r_f128 = fputil::quick_add(PI_F128, r_f128);
+
+ return static_cast<double>(r_f128);
+#endif // LIBC_MATH_HAS_SKIP_ACCURATE_PASS
+}
+
+} // namespace LIBC_NAMESPACE_DECL
diff --git a/libc/src/math/generic/asin.cpp b/libc/src/math/generic/asin.cpp
index ee85cb1988bb2..ad77683d1f880 100644
--- a/libc/src/math/generic/asin.cpp
+++ b/libc/src/math/generic/asin.cpp
@@ -74,7 +74,7 @@ LLVM_LIBC_FUNCTION(double, asin, (double x)) {
#else
unsigned idx;
DoubleDouble x_sq = fputil::exact_mult(x, x);
- double err = x * 0x1.0p-51;
+ double err = xbits.abs().get_val() * 0x1.0p-51;
// Polynomial approximation:
// p ~ asin(x)/x
@@ -135,12 +135,14 @@ LLVM_LIBC_FUNCTION(double, asin, (double x)) {
x_sign * PI_OVER_TWO.lo);
}
// |x| > 1, return NaN.
- if (xbits.is_finite()) {
+ if (xbits.is_quiet_nan())
+ return x;
+
+ // Set domain error for non-NaN input.
+ if (!xbits.is_nan())
fputil::set_errno_if_required(EDOM);
- fputil::raise_except_if_required(FE_INVALID);
- } else if (xbits.is_signaling_nan()) {
- fputil::raise_except_if_required(FE_INVALID);
- }
+
+ fputil::raise_except_if_required(FE_INVALID);
return FPBits::quiet_nan().get_val();
}
diff --git a/libc/src/math/generic/asin_utils.h b/libc/src/math/generic/asin_utils.h
index 0bdf25367f820..44913d573de2c 100644
--- a/libc/src/math/generic/asin_utils.h
+++ b/libc/src/math/generic/asin_utils.h
@@ -25,6 +25,8 @@ namespace {
using DoubleDouble = fputil::DoubleDouble;
using Float128 = fputil::DyadicFloat<128>;
+constexpr DoubleDouble PI = {0x1.1a62633145c07p-53, 0x1.921fb54442d18p1};
+
constexpr DoubleDouble PI_OVER_TWO = {0x1.1a62633145c07p-54,
0x1.921fb54442d18p0};
@@ -172,7 +174,8 @@ LIBC_INLINE DoubleDouble asin_eval(const DoubleDouble &u, unsigned &idx,
double y_hi = multiply_add(k, -0x1.0p-5, u.hi); // Exact
DoubleDouble y = fputil::exact_add(y_hi, u.lo);
double y2 = y.hi * y.hi;
- err *= y2 + 0x1.0p-30;
+ // Add double-double errors in addition to the relative errors from y2.
+ err = fputil::multiply_add(err, y2, 0x1.0p-102);
DoubleDouble c0 = fputil::quick_mult(
y, DoubleDouble{ASIN_COEFFS[idx][3], ASIN_COEFFS[idx][2]});
double c1 = multiply_add(y.hi, ASIN_COEFFS[idx][5], ASIN_COEFFS[idx][4]);
@@ -548,6 +551,9 @@ constexpr Float128 ASIN_COEFFS_F128[17][16] = {
constexpr Float128 PI_OVER_TWO_F128 = {
Sign::POS, -127, 0xc90fdaa2'2168c234'c4c6628b'80dc1cd1_u128};
+constexpr Float128 PI_F128 = {Sign::POS, -126,
+ 0xc90fdaa2'2168c234'c4c6628b'80dc1cd1_u128};
+
LIBC_INLINE Float128 asin_eval(const Float128 &u, unsigned idx) {
return fputil::polyeval(u, ASIN_COEFFS_F128[idx][0], ASIN_COEFFS_F128[idx][1],
ASIN_COEFFS_F128[idx][2], ASIN_COEFFS_F128[idx][3],
diff --git a/libc/test/src/math/CMakeLists.txt b/libc/test/src/math/CMakeLists.txt
index fba6b6f133163..9074524403102 100644
--- a/libc/test/src/math/CMakeLists.txt
+++ b/libc/test/src/math/CMakeLists.txt
@@ -2280,6 +2280,17 @@ add_fp_unittest(
libc.src.__support.FPUtil.fp_bits
)
+add_fp_unittest(
+ acos_test
+ NEED_MPFR
+ SUITE
+ libc-math-unittests
+ SRCS
+ acos_test.cpp
+ DEPENDS
+ libc.src.math.acos
+)
+
add_fp_unittest(
acosf16_test
NEED_MPFR
diff --git a/libc/test/src/math/acos_test.cpp b/libc/test/src/math/acos_test.cpp
new file mode 100644
index 0000000000000..140488702f0bc
--- /dev/null
+++ b/libc/test/src/math/acos_test.cpp
@@ -0,0 +1,82 @@
+//===-- Unittests for acos ------------------------------------------------===//
+//
+// 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/acos.h"
+#include "test/UnitTest/FPMatcher.h"
+#include "test/UnitTest/Test.h"
+#include "utils/MPFRWrapper/MPFRUtils.h"
+
+using LlvmLibcAcosTest = LIBC_NAMESPACE::testing::FPTest<double>;
+
+namespace mpfr = LIBC_NAMESPACE::testing::mpfr;
+
+using LIBC_NAMESPACE::testing::tlog;
+
+TEST_F(LlvmLibcAcosTest, InDoubleRange) {
+ constexpr uint64_t COUNT = 123'451;
+ uint64_t START = FPBits(0x1.0p-60).uintval();
+ uint64_t STOP = FPBits(1.0).uintval();
+ uint64_t STEP = (STOP - START) / COUNT;
+
+ auto test = [&](mpfr::RoundingMode rounding_mode) {
+ mpfr::ForceRoundingMode __r(rounding_mode);
+ if (!__r.success)
+ return;
+
+ uint64_t fails = 0;
+ uint64_t count = 0;
+ uint64_t cc = 0;
+ double mx = 0.0, mr = 0.0;
+ double tol = 0.5;
+
+ for (uint64_t i = 0, v = START; i <= COUNT; ++i, v += STEP) {
+ double x = FPBits(v).get_val();
+ if (FPBits(v).is_inf_or_nan())
+ continue;
+ double result = LIBC_NAMESPACE::acos(x);
+ ++cc;
+ if (FPBits(result).is_inf_or_nan())
+ continue;
+
+ ++count;
+
+ if (!TEST_MPFR_MATCH_ROUNDING_SILENTLY(mpfr::Operation::Acos, x, result,
+ 0.5, rounding_mode)) {
+ ++fails;
+ while (!TEST_MPFR_MATCH_ROUNDING_SILENTLY(mpfr::Operation::Acos, x,
+ result, tol, rounding_mode)) {
+ mx = x;
+ mr = result;
+
+ if (tol > 1000.0)
+ break;
+
+ tol *= 2.0;
+ }
+ }
+ }
+ if (fails) {
+ tlog << " Acos failed: " << fails << "/" << count << "/" << cc
+ << " tests.\n";
+ tlog << " Max ULPs is at most: " << static_cast<uint64_t>(tol) << ".\n";
+ EXPECT_MPFR_MATCH(mpfr::Operation::Acos, mx, mr, 0.5, rounding_mode);
+ }
+ };
+
+ tlog << " Test Rounding To Nearest...\n";
+ test(mpfr::RoundingMode::Nearest);
+
+ tlog << " Test Rounding Downward...\n";
+ test(mpfr::RoundingMode::Downward);
+
+ tlog << " Test Rounding Upward...\n";
+ test(mpfr::RoundingMode::Upward);
+
+ tlog << " Test Rounding Toward Zero...\n";
+ test(mpfr::RoundingMode::TowardZero);
+}
diff --git a/libc/test/src/math/smoke/CMakeLists.txt b/libc/test/src/math/smoke/CMakeLists.txt
index ca1ce7fb1f48c..55edffb4f5bdf 100644
--- a/libc/test/src/math/smoke/CMakeLists.txt
+++ b/libc/test/src/math/smoke/CMakeLists.txt
@@ -4051,6 +4051,18 @@ add_fp_unittest(
libc.src.__support.FPUtil.fp_bits
)
+add_fp_unittest(
+ acos_test
+ SUITE
+ libc-math-smoke-tests
+ SRCS
+ acos_test.cpp
+ DEPENDS
+ libc.hdr.fenv_macros
+ libc.src.errno.errno
+ libc.src.math.acos
+)
+
add_fp_unittest(
acosf16_test
SUITE
diff --git a/libc/test/src/math/smoke/acos_test.cpp b/libc/test/src/math/smoke/acos_test.cpp
new file mode 100644
index 0000000000000..3a59bce264077
--- /dev/null
+++ b/libc/test/src/math/smoke/acos_test.cpp
@@ -0,0 +1,64 @@
+//===-- Unittests for acos ------------------------------------------------===//
+//
+// 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 "hdr/fenv_macros.h"
+#include "src/errno/libc_errno.h"
+#include "src/math/acos.h"
+#include "test/UnitTest/FPMatcher.h"
+#include "test/UnitTest/Test.h"
+
+using LlvmLibcAcosTest = LIBC_NAMESPACE::testing::FPTest<double>;
+
+TEST_F(LlvmLibcAcosTest, SpecialNumbers) {
+ EXPECT_FP_EQ_WITH_EXCEPTION_ALL_ROUNDING(aNaN, LIBC_NAMESPACE::acos(sNaN),
+ FE_INVALID);
+ EXPECT_FP_EQ_ALL_ROUNDING(aNaN, LIBC_NAMESPACE::acos(aNaN));
+ EXPECT_FP_EQ(0x1.921fb54442d18p0, LIBC_NAMESPACE::acos(zero));
+ EXPECT_FP_EQ(0x1.921fb54442d18p0, LIBC_NAMESPACE::acos(neg_zero));
+
+ LIBC_NAMESPACE::libc_errno = 0;
+ EXPECT_FP_EQ_WITH_EXCEPTION_ALL_ROUNDING(aNaN, LIBC_NAMESPACE::acos(inf),
+ FE_INVALID);
+ EXPECT_MATH_ERRNO(EDOM);
+ EXPECT_FP_EQ_WITH_EXCEPTION_ALL_ROUNDING(aNaN, LIBC_NAMESPACE::acos(neg_inf),
+ FE_INVALID);
+ EXPECT_MATH_ERRNO(EDOM);
+ EXPECT_FP_EQ_WITH_EXCEPTION_ALL_ROUNDING(aNaN, LIBC_NAMESPACE::acos(2.0),
+ FE_INVALID);
+ EXPECT_MATH_ERRNO(EDOM);
+ EXPECT_FP_EQ_WITH_EXCEPTION_ALL_ROUNDING(aNaN, LIBC_NAMESPACE::acos(-2.0),
+ FE_INVALID);
+ EXPECT_MATH_ERRNO(EDOM);
+ EXPECT_FP_EQ(zero, LIBC_NAMESPACE::acos(1.0));
+ EXPECT_FP_EQ(0x1.921fb54442d18p1, LIBC_NAMESPACE::acos(-1.0));
+ EXPECT_FP_EQ(0x1.921fb54442d18p0, LIBC_NAMESPACE::acos(0x1.0p-54));
+}
+
+#ifdef LIBC_TEST_FTZ_DAZ
+
+using namespace LIBC_NAMESPACE::testing;
+
+TEST_F(LlvmLibcAcosTest, FTZMode) {
+ ModifyMXCSR mxcsr(FTZ);
+
+ EXPECT_FP_EQ(0x1.921fb54442d18p0, LIBC_NAMESPACE::acos(min_denormal));
+}
+
+TEST_F(LlvmLibcAcosTest, DAZMode) {
+ ModifyMXCSR mxcsr(DAZ);
+
+ EXPECT_FP_EQ(0x1.921fb54442d18p0, LIBC_NAMESPACE::acos(min_denormal));
+}
+
+TEST_F(LlvmLibcAcosTest, FTZDAZMode) {
+ ModifyMXCSR mxcsr(FTZ | DAZ);
+
+ EXPECT_FP_EQ(0x1.921fb54442d18p0, LIBC_NAMESPACE::acos(min_denormal));
+}
+
+#endif
More information about the libc-commits
mailing list