diff --git a/libc/config/darwin/arm/entrypoints.txt b/libc/config/darwin/arm/entrypoints.txt index 70c888aec064..308fc49d681d 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 66d7576ffb88..520046f768b5 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 b5e2f59d25a5..7432a7e912e8 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 30d9d00dfefc..c2a31b9f5c96 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 2ac016f065b2..9f447dd0d35d 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 37fa888d6498..09021a08cf73 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 9d4d6698b612..6b0365f481a4 100644 --- a/libc/docs/headers/math/index.rst +++ b/libc/docs/headers/math/index.rst @@ -249,7 +249,7 @@ Higher Math Functions +-----------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+----------------------------+ | | (float) | (double) | (long double) | (float16) | (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 f8c57f7d48ab..0a69653d156a 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 b08486ac6136..e5d1f9075bcb 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 000000000000..c14721faef3c --- /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; + + 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(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(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(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(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(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(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(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(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 ee85cb1988bb..ad77683d1f88 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 0bdf25367f82..44913d573de2 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 fba6b6f13316..907452440310 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 000000000000..140488702f0b --- /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; + +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(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 ca1ce7fb1f48..55edffb4f5bd 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 000000000000..3a59bce26407 --- /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; + +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