diff --git a/libc/config/darwin/arm/entrypoints.txt b/libc/config/darwin/arm/entrypoints.txt --- a/libc/config/darwin/arm/entrypoints.txt +++ b/libc/config/darwin/arm/entrypoints.txt @@ -105,6 +105,7 @@ 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 --- a/libc/config/linux/aarch64/entrypoints.txt +++ b/libc/config/linux/aarch64/entrypoints.txt @@ -149,6 +149,7 @@ 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 --- a/libc/config/linux/x86_64/entrypoints.txt +++ b/libc/config/linux/x86_64/entrypoints.txt @@ -149,6 +149,7 @@ 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 --- a/libc/config/windows/entrypoints.txt +++ b/libc/config/windows/entrypoints.txt @@ -106,6 +106,7 @@ 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 --- a/libc/docs/math.rst +++ b/libc/docs/math.rst @@ -118,7 +118,7 @@ ============== ================ =============== ====================== acos acosh -asin +asin |check| asinh atan |check| atan2 @@ -154,6 +154,7 @@ ============== ================ =============== ====================== (float) (double) (long double) ============== ================ =============== ====================== +asin |check| atan |check| atanh |check| cos |check| large @@ -203,6 +204,8 @@ | +-----------+-------------------+-----------+-------------------+ +------------+-------------------------+--------------+---------------+ | | 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 --- a/libc/spec/stdc.td +++ b/libc/spec/stdc.td @@ -480,6 +480,7 @@ FunctionSpec<"sinhf", RetValSpec, [ArgSpec]>, FunctionSpec<"tanhf", RetValSpec, [ArgSpec]>, + FunctionSpec<"asinf", RetValSpec, [ArgSpec]>, FunctionSpec<"atanf", RetValSpec, [ArgSpec]>, FunctionSpec<"atanhf", RetValSpec, [ArgSpec]>, diff --git a/libc/src/math/CMakeLists.txt b/libc/src/math/CMakeLists.txt --- a/libc/src/math/CMakeLists.txt +++ b/libc/src/math/CMakeLists.txt @@ -64,7 +64,9 @@ -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 --- /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 --- a/libc/src/math/generic/CMakeLists.txt +++ b/libc/src/math/generic/CMakeLists.txt @@ -1288,6 +1288,22 @@ 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 --- /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 + +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 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 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; + 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(x); + return static_cast(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(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::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(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 --- a/libc/test/src/math/CMakeLists.txt +++ b/libc/test/src/math/CMakeLists.txt @@ -1437,6 +1437,20 @@ 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 --- /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 + +#include +#include + +using FPBits = __llvm_libc::fputil::FPBits; + +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 --- a/libc/test/src/math/exhaustive/CMakeLists.txt +++ b/libc/test/src/math/exhaustive/CMakeLists.txt @@ -325,3 +325,19 @@ -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 --- /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 + +using FPBits = __llvm_libc::fputil::FPBits; + +namespace mpfr = __llvm_libc::testing::mpfr; + +struct LlvmLibcAsinfExhaustiveTest : 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::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.h b/libc/utils/MPFRWrapper/MPFRUtils.h --- a/libc/utils/MPFRWrapper/MPFRUtils.h +++ b/libc/utils/MPFRWrapper/MPFRUtils.h @@ -25,6 +25,7 @@ // and output floating point numbers are of the same kind. BeginUnaryOperationsSingleOutput, Abs, + Asin, Atan, Atanh, Ceil, diff --git a/libc/utils/MPFRWrapper/MPFRUtils.cpp b/libc/utils/MPFRWrapper/MPFRUtils.cpp --- a/libc/utils/MPFRWrapper/MPFRUtils.cpp +++ b/libc/utils/MPFRWrapper/MPFRUtils.cpp @@ -176,12 +176,20 @@ 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 @@ if (thisAsT == input) return T(0.0); + if (is_nan()) { + if (fputil::FPBits(input).is_nan()) + return T(0.0); + return T(fputil::FPBits::inf()); + } + int thisExponent = fputil::FPBits(thisAsT).get_exponent(); int inputExponent = fputil::FPBits(input).get_exponent(); // Adjust the exponents for denormal numbers. @@ -512,6 +526,8 @@ switch (op) { case Operation::Abs: return mpfrInput.abs(); + case Operation::Asin: + return mpfrInput.asin(); case Operation::Atan: return mpfrInput.atan(); case Operation::Atanh: