From: Tue Ly Date: Fri, 9 Sep 2022 02:30:37 +0000 (-0400) Subject: [libc][math] Implement acosf function correctly rounded for all rounding modes. X-Git-Tag: upstream/17.0.6~33943 X-Git-Url: http://review.tizen.org/git/?a=commitdiff_plain;h=463dcc8749ed713df858dd1738e472725c8b1b35;p=platform%2Fupstream%2Fllvm.git [libc][math] Implement acosf function correctly rounded for all rounding modes. Implement acosf function correctly rounded for all rounding modes. We perform range reduction as follows: - When `|x| < 2^(-10)`, we use cubic Taylor polynomial: ``` acos(x) = pi/2 - asin(x) ~ pi/2 - x - x^3 / 6. ``` - When `2^(-10) <= |x| <= 0.5`, we use the same approximation that is used for `asinf(x)` when `|x| <= 0.5`: ``` acos(x) = pi/2 - asin(x) ~ pi/2 - x - x^3 * P(x^2). ``` - When `0.5 < x <= 1`, we use the double angle formula: `cos(2y) = 1 - 2 * sin^2 (y)` to reduce to: ``` acos(x) = 2 * asin( sqrt( (1 - x)/2 ) ) ``` - When `-1 <= x < -0.5`, we reduce to the positive case above using the formula: ``` acos(x) = pi - acos(-x) ``` Performance benchmark using perf tool from the CORE-MATH project on Ryzen 1700: ``` $ CORE_MATH_PERF_MODE="rdtsc" ./perf.sh acosf GNU libc version: 2.35 GNU libc release: stable CORE-MATH reciprocal throughput : 28.613 System LIBC reciprocal throughput : 29.204 LIBC reciprocal throughput : 24.271 $ CORE_MATH_PERF_MODE="rdtsc" ./perf.sh asinf --latency GNU libc version: 2.35 GNU libc release: stable CORE-MATH latency : 55.554 System LIBC latency : 76.879 LIBC latency : 62.118 ``` Reviewed By: orex, zimmermann6 Differential Revision: https://reviews.llvm.org/D133550 --- diff --git a/libc/config/darwin/arm/entrypoints.txt b/libc/config/darwin/arm/entrypoints.txt index 433c245..02331e1 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.acosf libc.src.math.asinf libc.src.math.atanf libc.src.math.atanhf diff --git a/libc/config/linux/aarch64/entrypoints.txt b/libc/config/linux/aarch64/entrypoints.txt index abc6d77..fc64f1e 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.acosf libc.src.math.asinf libc.src.math.atanf libc.src.math.atanhf diff --git a/libc/config/linux/x86_64/entrypoints.txt b/libc/config/linux/x86_64/entrypoints.txt index 60ec5f07..b1fb34c 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.acosf libc.src.math.asinf libc.src.math.atanf libc.src.math.atanhf diff --git a/libc/config/windows/entrypoints.txt b/libc/config/windows/entrypoints.txt index a679181..4bd46c5 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.acosf libc.src.math.asinf libc.src.math.atanf libc.src.math.atanhf diff --git a/libc/docs/math.rst b/libc/docs/math.rst index 0cbd533..55c8748 100644 --- a/libc/docs/math.rst +++ b/libc/docs/math.rst @@ -116,7 +116,7 @@ Higher Math Functions ============== ================ =============== ====================== (float) (double) (long double) ============== ================ =============== ====================== -acos +acos |check| acosh asin |check| asinh @@ -154,6 +154,7 @@ Accuracy of Higher Math Functions ============== ================ =============== ====================== (float) (double) (long double) ============== ================ =============== ====================== +acos |check| asin |check| atan |check| atanh |check| @@ -204,6 +205,8 @@ Performance | +-----------+-------------------+-----------+-------------------+ +------------+-------------------------+--------------+---------------+ | | LLVM libc | Reference (glibc) | LLVM libc | Reference (glibc) | | CPU | OS | Compiler | Special flags | +==============+===========+===================+===========+===================+=====================================+============+=========================+==============+===============+ +| acosf | 24 | 29 | 62 | 77 | :math:`[-1, 1]` | Ryzen 1700 | Ubuntu 22.04 LTS x86_64 | Clang 14.0.0 | FMA | ++--------------+-----------+-------------------+-----------+-------------------+-------------------------------------+------------+-------------------------+--------------+---------------+ | 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 | diff --git a/libc/spec/stdc.td b/libc/spec/stdc.td index fdb6beb..6ad162f 100644 --- a/libc/spec/stdc.td +++ b/libc/spec/stdc.td @@ -480,6 +480,7 @@ def StdC : StandardSpec<"stdc"> { FunctionSpec<"sinhf", RetValSpec, [ArgSpec]>, FunctionSpec<"tanhf", RetValSpec, [ArgSpec]>, + FunctionSpec<"acosf", RetValSpec, [ArgSpec]>, FunctionSpec<"asinf", RetValSpec, [ArgSpec]>, FunctionSpec<"atanf", RetValSpec, [ArgSpec]>, diff --git a/libc/src/math/CMakeLists.txt b/libc/src/math/CMakeLists.txt index 313ff6a..cb2a507 100644 --- a/libc/src/math/CMakeLists.txt +++ b/libc/src/math/CMakeLists.txt @@ -64,6 +64,7 @@ add_entrypoint_object( -O3 ) +add_math_entrypoint_object(acosf) add_math_entrypoint_object(asinf) add_math_entrypoint_object(atanf) diff --git a/libc/src/math/acosf.h b/libc/src/math/acosf.h new file mode 100644 index 0000000..9604fbf --- /dev/null +++ b/libc/src/math/acosf.h @@ -0,0 +1,18 @@ +//===-- Implementation header for acosf -------------------------*- 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_ACOSF_H +#define LLVM_LIBC_SRC_MATH_ACOSF_H + +namespace __llvm_libc { + +float acosf(float x); + +} // namespace __llvm_libc + +#endif // LLVM_LIBC_SRC_MATH_ACOSF_H diff --git a/libc/src/math/generic/CMakeLists.txt b/libc/src/math/generic/CMakeLists.txt index a946ae4..ffc06c8 100644 --- a/libc/src/math/generic/CMakeLists.txt +++ b/libc/src/math/generic/CMakeLists.txt @@ -1300,6 +1300,24 @@ add_entrypoint_object( libc.src.__support.FPUtil.multiply_add libc.src.__support.FPUtil.polyeval libc.src.__support.FPUtil.sqrt + .inv_trigf_utils + COMPILE_OPTIONS + -O3 +) + +add_entrypoint_object( + acosf + SRCS + acosf.cpp + HDRS + ../acosf.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 + .inv_trigf_utils COMPILE_OPTIONS -O3 ) diff --git a/libc/src/math/generic/acosf.cpp b/libc/src/math/generic/acosf.cpp new file mode 100644 index 0000000..73c93cb --- /dev/null +++ b/libc/src/math/generic/acosf.cpp @@ -0,0 +1,117 @@ +//===-- Single-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/acosf.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 + +#include "inv_trigf_utils.h" + +namespace __llvm_libc { + +static constexpr size_t N_EXCEPTS = 4; + +// Exceptional values when |x| <= 0.5 +static constexpr fputil::ExceptValues ACOSF_EXCEPTS = {{ + // (inputs, RZ output, RU offset, RD offset, RN offset) + // x = 0x1.110b46p-26, acosf(x) = 0x1.921fb4p0 (RZ) + {0x328885a3, 0x3fc90fda, 1, 0, 1}, + // x = -0x1.110b46p-26, acosf(x) = 0x1.921fb4p0 (RZ) + {0xb28885a3, 0x3fc90fda, 1, 0, 1}, + // x = 0x1.04c444p-12, acosf(x) = 0x1.920f68p0 (RZ) + {0x39826222, 0x3fc907b4, 1, 0, 1}, + // x = -0x1.04c444p-12, acosf(x) = 0x1.923p0 (RZ) + {0xb9826222, 0x3fc91800, 1, 0, 1}, +}}; + +LLVM_LIBC_FUNCTION(float, acosf, (float x)) { + using FPBits = typename fputil::FPBits; + FPBits xbits(x); + uint32_t x_uint = xbits.uintval(); + uint32_t x_abs = xbits.uintval() & 0x7fff'ffffU; + uint32_t x_sign = x_uint >> 31; + + // |x| <= 0.5 + if (unlikely(x_abs <= 0x3f00'0000U)) { + // |x| < 0x1p-10 + if (unlikely(x_abs < 0x3a80'0000U)) { + // When |x| < 2^-10, we use the following approximation: + // acos(x) = pi/2 - asin(x) + // ~ pi/2 - x - x^3 / 6 + + // Check for exceptional values + if (auto r = ACOSF_EXCEPTS.lookup(x_uint); unlikely(r.has_value())) + return r.value(); + + double xd = static_cast(x); + return fputil::multiply_add(-0x1.5555555555555p-3 * xd, xd * xd, + M_MATH_PI_2 - xd); + } + + // For |x| <= 0.5, we approximate acosf(x) by: + // acos(x) = pi/2 - asin(x) = pi/2 - 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]); + double xd = static_cast(x); + double xsq = xd * xd; + double x3 = xd * xsq; + double r = asin_eval(xsq); + return fputil::multiply_add(-x3, r, M_MATH_PI_2 - 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::VALUE - 1)); + } + + // When 0.5 < |x| <= 1, we perform range reduction as follow: + // + // Assume further that 0.5 < x <= 1, and let: + // y = acos(x) + // We use the double angle formula: + // 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 ) ) + // 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 use the identity: + // acos(x) = pi - acos(-x) + // which is reduced to the postive case. + + xbits.set_sign(false); + double xd = static_cast(xbits.get_val()); + double u = fputil::multiply_add(-0.5, xd, 0.5); + double cv = 2 * fputil::sqrt(u); + + double r3 = asin_eval(u); + double r = fputil::multiply_add(cv * u, r3, cv); + return x_sign ? M_MATH_PI - r : r; +} + +} // namespace __llvm_libc diff --git a/libc/src/math/generic/asinf.cpp b/libc/src/math/generic/asinf.cpp index 262af0f..7314b38 100644 --- a/libc/src/math/generic/asinf.cpp +++ b/libc/src/math/generic/asinf.cpp @@ -16,10 +16,9 @@ #include -namespace __llvm_libc { +#include "inv_trigf_utils.h" -// PI / 2 -constexpr double M_MATH_PI_2 = 0x1.921fb54442d18p+0; +namespace __llvm_libc { static constexpr size_t N_EXCEPTS = 2; @@ -41,14 +40,6 @@ static constexpr fputil::ExceptValues ASINF_EXCEPTS_HI = {{ {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; FPBits xbits(x); @@ -104,14 +95,9 @@ LLVM_LIBC_FUNCTION(float, asinf, (float x)) { // little more than 0.5. double xd = static_cast(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); + double r = asin_eval(xsq); + return fputil::multiply_add(x3, r, xd); } // |x| > 1, return NaNs. @@ -130,6 +116,7 @@ LLVM_LIBC_FUNCTION(float, asinf, (float x)) { 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: @@ -143,27 +130,24 @@ LLVM_LIBC_FUNCTION(float, asinf, (float x)) { // 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, + // Let u = (1 - x)/2, then: + // asin(x) = pi/2 - 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: - // asin(x) = pi/2 - 2 * u * P(u^2), + // asin(x) ~ pi/2 - 2 * sqrt(u) * P(u), xbits.set_sign(false); + double sign = SIGN[x_sign]; double xd = static_cast(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; + double c1 = sign * (-2 * fputil::sqrt(u)); + double c2 = fputil::multiply_add(sign, M_MATH_PI_2, c1); + double c3 = c1 * u; + + double r = asin_eval(u); + return fputil::multiply_add(c3, r, c2); } } // namespace __llvm_libc diff --git a/libc/src/math/generic/inv_trigf_utils.h b/libc/src/math/generic/inv_trigf_utils.h index 32f0ba7..202a4a6 100644 --- a/libc/src/math/generic/inv_trigf_utils.h +++ b/libc/src/math/generic/inv_trigf_utils.h @@ -21,7 +21,8 @@ namespace __llvm_libc { -// PI / 2 +// PI and PI / 2 +constexpr double M_MATH_PI = 0x1.921fb54442d18p+1; constexpr double M_MATH_PI_2 = 0x1.921fb54442d18p+0; // atan table size @@ -36,7 +37,7 @@ extern const double ATAN_K[5]; // atan(u) + atan(v) = atan((u+v)/(1-uv)) // x should be positive, normal finite value -inline static double atan_eval(double x) { +static inline double atan_eval(double x) { using FPB = fputil::FPBits; // Added some small value to umin and umax mantissa to avoid possible rounding // errors. @@ -89,6 +90,24 @@ inline static double atan_eval(double x) { return sign ? -result : result; } +// > Q = fpminimax(asin(x)/x, [|0, 2, 4, 6, 8, 10, 12, 14, 16, 18, 20|], +// [|1, D...|], [0, 0.5]); +constexpr double ASIN_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}; + +// Evaluate P(x^2) - 1, where P(x^2) ~ asin(x)/x +static inline double asin_eval(double xsq) { + double x4 = xsq * xsq; + double r1 = fputil::polyeval(x4, ASIN_COEFFS[0], ASIN_COEFFS[2], + ASIN_COEFFS[4], ASIN_COEFFS[6], ASIN_COEFFS[8]); + double r2 = fputil::polyeval(x4, ASIN_COEFFS[1], ASIN_COEFFS[3], + ASIN_COEFFS[5], ASIN_COEFFS[7], ASIN_COEFFS[9]); + return fputil::multiply_add(xsq, r2, r1); +} + } // namespace __llvm_libc #endif // LLVM_LIBC_SRC_MATH_GENERIC_INV_TRIGF_UTILS_H diff --git a/libc/test/src/math/CMakeLists.txt b/libc/test/src/math/CMakeLists.txt index b566905..0cce635 100644 --- a/libc/test/src/math/CMakeLists.txt +++ b/libc/test/src/math/CMakeLists.txt @@ -1452,6 +1452,20 @@ add_fp_unittest( ) add_fp_unittest( + acosf_test + NEED_MPFR + SUITE + libc_math_unittests + SRCS + acosf_test.cpp + DEPENDS + libc.include.errno + libc.src.errno.errno + libc.src.math.acosf + libc.src.__support.FPUtil.fp_bits +) + +add_fp_unittest( atanf_test NEED_MPFR SUITE diff --git a/libc/test/src/math/acosf_test.cpp b/libc/test/src/math/acosf_test.cpp new file mode 100644 index 0000000..fa96537 --- /dev/null +++ b/libc/test/src/math/acosf_test.cpp @@ -0,0 +1,75 @@ +//===-- Unittests for acosf -----------------------------------------------===// +// +// 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/acosf.h" +#include "utils/MPFRWrapper/MPFRUtils.h" +#include "utils/UnitTest/FPMatcher.h" +#include "utils/UnitTest/Test.h" +#include + +#include +#include + +using FPBits = __llvm_libc::fputil::FPBits; + +namespace mpfr = __llvm_libc::testing::mpfr; + +DECLARE_SPECIAL_CONSTANTS(float) + +TEST(LlvmLibcAcosfTest, SpecialNumbers) { + errno = 0; + + EXPECT_FP_EQ(aNaN, __llvm_libc::acosf(aNaN)); + EXPECT_MATH_ERRNO(0); + + EXPECT_FP_EQ(aNaN, __llvm_libc::acosf(inf)); + EXPECT_MATH_ERRNO(EDOM); + + EXPECT_FP_EQ(aNaN, __llvm_libc::acosf(neg_inf)); + EXPECT_MATH_ERRNO(EDOM); +} + +TEST(LlvmLibcAcosfTest, 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::Acos, x, + __llvm_libc::acosf(x), 0.5); + } +} + +TEST(LlvmLibcAcosfTest, SpecificBitPatterns) { + constexpr int N = 13; + 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 + 0x328885a3, // x = 0x1.110b46p-26 + 0x39826222, // x = 0x1.04c444p-12 + 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::Acos, x, + __llvm_libc::acosf(x), 0.5); + EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Acos, -x, + __llvm_libc::acosf(-x), 0.5); + } +} diff --git a/libc/test/src/math/exhaustive/CMakeLists.txt b/libc/test/src/math/exhaustive/CMakeLists.txt index 1051e87..819e789 100644 --- a/libc/test/src/math/exhaustive/CMakeLists.txt +++ b/libc/test/src/math/exhaustive/CMakeLists.txt @@ -341,3 +341,20 @@ add_fp_unittest( LINK_LIBRARIES -lpthread ) + +add_fp_unittest( + acosf_test + NO_RUN_POSTBUILD + NEED_MPFR + SUITE + libc_math_exhaustive_tests + SRCS + acosf_test.cpp + DEPENDS + .exhaustive_test + libc.include.math + libc.src.math.acosf + libc.src.__support.FPUtil.fp_bits + LINK_LIBRARIES + -lpthread +) diff --git a/libc/test/src/math/exhaustive/acosf_test.cpp b/libc/test/src/math/exhaustive/acosf_test.cpp new file mode 100644 index 0000000..019232b --- /dev/null +++ b/libc/test/src/math/exhaustive/acosf_test.cpp @@ -0,0 +1,76 @@ +//===-- Exhaustive test for acosf -----------------------------------------===// +// +// 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/acosf.h" +#include "utils/MPFRWrapper/MPFRUtils.h" + +#include + +using FPBits = __llvm_libc::fputil::FPBits; + +namespace mpfr = __llvm_libc::testing::mpfr; + +struct LlvmLibcAcosfExhaustiveTest : public LlvmLibcExhaustiveTest { + 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::Acos, x, + __llvm_libc::acosf(x), 0.5, rounding); + } 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(LlvmLibcAcosfExhaustiveTest, PostiveRangeRoundNearestTieToEven) { + test_full_range(POS_START, POS_STOP, mpfr::RoundingMode::Nearest); +} + +TEST_F(LlvmLibcAcosfExhaustiveTest, PostiveRangeRoundUp) { + test_full_range(POS_START, POS_STOP, mpfr::RoundingMode::Upward); +} + +TEST_F(LlvmLibcAcosfExhaustiveTest, PostiveRangeRoundDown) { + test_full_range(POS_START, POS_STOP, mpfr::RoundingMode::Downward); +} + +TEST_F(LlvmLibcAcosfExhaustiveTest, PostiveRangeRoundTowardZero) { + test_full_range(POS_START, POS_STOP, mpfr::RoundingMode::TowardZero); +} + +// Range: [-Inf, 0]; +static const uint32_t NEG_START = 0xb000'0000U; +static const uint32_t NEG_STOP = 0xff80'0000U; + +TEST_F(LlvmLibcAcosfExhaustiveTest, NegativeRangeRoundNearestTieToEven) { + test_full_range(NEG_START, NEG_STOP, mpfr::RoundingMode::Nearest); +} + +TEST_F(LlvmLibcAcosfExhaustiveTest, NegativeRangeRoundUp) { + test_full_range(NEG_START, NEG_STOP, mpfr::RoundingMode::Upward); +} + +TEST_F(LlvmLibcAcosfExhaustiveTest, NegativeRangeRoundDown) { + test_full_range(NEG_START, NEG_STOP, mpfr::RoundingMode::Downward); +} + +TEST_F(LlvmLibcAcosfExhaustiveTest, NegativeRangeRoundTowardZero) { + test_full_range(NEG_START, NEG_STOP, mpfr::RoundingMode::TowardZero); +} diff --git a/libc/test/src/math/exhaustive/asinf_test.cpp b/libc/test/src/math/exhaustive/asinf_test.cpp index fa16acb..0c6b5c4 100644 --- a/libc/test/src/math/exhaustive/asinf_test.cpp +++ b/libc/test/src/math/exhaustive/asinf_test.cpp @@ -39,9 +39,8 @@ 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_F(LlvmLibcAsinfExhaustiveTest, PostiveRangeRoundNearestTieToEven) { test_full_range(POS_START, POS_STOP, mpfr::RoundingMode::Nearest); } diff --git a/libc/utils/MPFRWrapper/MPFRUtils.cpp b/libc/utils/MPFRWrapper/MPFRUtils.cpp index 910233e..0f71357 100644 --- a/libc/utils/MPFRWrapper/MPFRUtils.cpp +++ b/libc/utils/MPFRWrapper/MPFRUtils.cpp @@ -184,6 +184,12 @@ public: return result; } + MPFRNumber acos() const { + MPFRNumber result(*this); + mpfr_acos(result.value, value, mpfr_rounding); + return result; + } + MPFRNumber asin() const { MPFRNumber result(*this); mpfr_asin(result.value, value, mpfr_rounding); @@ -526,6 +532,8 @@ unary_operation(Operation op, InputType input, unsigned int precision, switch (op) { case Operation::Abs: return mpfrInput.abs(); + case Operation::Acos: + return mpfrInput.acos(); case Operation::Asin: return mpfrInput.asin(); case Operation::Atan: diff --git a/libc/utils/MPFRWrapper/MPFRUtils.h b/libc/utils/MPFRWrapper/MPFRUtils.h index 2f9f6f0..e8d5eca 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, + Acos, Asin, Atan, Atanh,