diff --git a/libc/shared/math.h b/libc/shared/math.h index 43189e5251d4a..042daf65e45bd 100644 --- a/libc/shared/math.h +++ b/libc/shared/math.h @@ -18,6 +18,7 @@ #include "math/acoshf16.h" #include "math/acospif16.h" #include "math/asin.h" +#include "math/asinf.h" #include "math/erff.h" #include "math/exp.h" #include "math/exp10.h" diff --git a/libc/shared/math/asinf.h b/libc/shared/math/asinf.h new file mode 100644 index 0000000000000..ac051bd2044af --- /dev/null +++ b/libc/shared/math/asinf.h @@ -0,0 +1,23 @@ +//===-- Shared asinf function -----------------------------------*- 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_SHARED_MATH_ASINF_H +#define LLVM_LIBC_SHARED_MATH_ASINF_H + +#include "shared/libc_common.h" +#include "src/__support/math/asinf.h" + +namespace LIBC_NAMESPACE_DECL { +namespace shared { + +using math::asinf; + +} // namespace shared +} // namespace LIBC_NAMESPACE_DECL + +#endif // LLVM_LIBC_SHARED_MATH_ASINF_H diff --git a/libc/src/__support/math/CMakeLists.txt b/libc/src/__support/math/CMakeLists.txt index beeb170f16f70..b096c616b0927 100644 --- a/libc/src/__support/math/CMakeLists.txt +++ b/libc/src/__support/math/CMakeLists.txt @@ -140,6 +140,22 @@ add_header_library( libc.src.__support.macros.properties.cpu_features ) +add_header_library( + asinf + HDRS + asinf.h + DEPENDS + .inv_trigf_utils + libc.src.__support.FPUtil.fenv_impl + libc.src.__support.FPUtil.fp_bits + libc.src.__support.FPUtil.except_value_utils + libc.src.__support.FPUtil.multiply_add + libc.src.__support.FPUtil.sqrt + libc.src.__support.macros.config + libc.src.__support.macros.optimization + libc.src.__support.macros.properties.cpu_features +) + add_header_library( erff HDRS diff --git a/libc/src/__support/math/asinf.h b/libc/src/__support/math/asinf.h new file mode 100644 index 0000000000000..bfa0dc31ecf4c --- /dev/null +++ b/libc/src/__support/math/asinf.h @@ -0,0 +1,175 @@ +//===-- 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___SUPPORT_MATH_ASINF_H +#define LLVM_LIBC_SRC___SUPPORT_MATH_ASINF_H + +#include "inv_trigf_utils.h" +#include "src/__support/FPUtil/FEnvImpl.h" +#include "src/__support/FPUtil/FPBits.h" +#include "src/__support/FPUtil/except_value_utils.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 { + +namespace math { + +LIBC_INLINE static constexpr float asinf(float x) { + using namespace inv_trigf_utils_internal; + +#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS + constexpr size_t N_EXCEPTS = 2; + + // Exceptional values when |x| <= 0.5 + 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 + 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}, + }}; +#endif // !LIBC_MATH_HAS_SKIP_ACCURATE_PASS + + using namespace inv_trigf_utils_internal; + 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 (LIBC_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_CPU_HAS_FMA_FLOAT) + 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_CPU_HAS_FMA_FLOAT + } + +#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS + // Check for exceptional values + if (auto r = ASINF_EXCEPTS_LO.lookup_odd(x_abs, x_sign); + LIBC_UNLIKELY(r.has_value())) + return r.value(); +#endif // !LIBC_MATH_HAS_SKIP_ACCURATE_PASS + + // 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 x3 = xd * xsq; + double r = asin_eval(xsq); + return static_cast(fputil::multiply_add(x3, r, xd)); + } + + // |x| > 1, return NaNs. + if (LIBC_UNLIKELY(x_abs > 0x3f80'0000U)) { + if (xbits.is_signaling_nan()) { + fputil::raise_except_if_required(FE_INVALID); + return FPBits::quiet_nan().get_val(); + } + + if (x_abs <= 0x7f80'0000U) { + fputil::set_errno_if_required(EDOM); + fputil::raise_except_if_required(FE_INVALID); + } + + return FPBits::quiet_nan().get_val(); + } + +#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS + // Check for exceptional values + if (auto r = ASINF_EXCEPTS_HI.lookup_odd(x_abs, x_sign); + LIBC_UNLIKELY(r.has_value())) + return r.value(); +#endif // !LIBC_MATH_HAS_SKIP_ACCURATE_PASS + + // 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( 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 * sqrt(u) * P(u), + + xbits.set_sign(Sign::POS); + double sign = SIGN[x_sign]; + double xd = static_cast(xbits.get_val()); + double u = fputil::multiply_add(-0.5, xd, 0.5); + 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 static_cast(fputil::multiply_add(c3, r, c2)); +} + +} // namespace math + +} // namespace LIBC_NAMESPACE_DECL + +#endif // LLVM_LIBC_SRC___SUPPORT_MATH_ASINF_H diff --git a/libc/src/math/generic/CMakeLists.txt b/libc/src/math/generic/CMakeLists.txt index 17e17a890f70e..ecf09674b4d4c 100644 --- a/libc/src/math/generic/CMakeLists.txt +++ b/libc/src/math/generic/CMakeLists.txt @@ -3958,13 +3958,7 @@ add_entrypoint_object( 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 - libc.src.__support.macros.optimization - libc.src.__support.math.inv_trigf_utils + libc.src.__support.math.asinf ) add_entrypoint_object( diff --git a/libc/src/math/generic/asinf.cpp b/libc/src/math/generic/asinf.cpp index 77d6de910962c..9c6766f29d12b 100644 --- a/libc/src/math/generic/asinf.cpp +++ b/libc/src/math/generic/asinf.cpp @@ -7,161 +7,10 @@ //===----------------------------------------------------------------------===// #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 "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 - -#include "src/__support/math/inv_trigf_utils.h" +#include "src/__support/math/asinf.h" namespace LIBC_NAMESPACE_DECL { -#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS -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}, -}}; -#endif // !LIBC_MATH_HAS_SKIP_ACCURATE_PASS - -LLVM_LIBC_FUNCTION(float, asinf, (float x)) { - using namespace inv_trigf_utils_internal; - 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 (LIBC_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_CPU_HAS_FMA_FLOAT) - 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_CPU_HAS_FMA_FLOAT - } - -#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS - // Check for exceptional values - if (auto r = ASINF_EXCEPTS_LO.lookup_odd(x_abs, x_sign); - LIBC_UNLIKELY(r.has_value())) - return r.value(); -#endif // !LIBC_MATH_HAS_SKIP_ACCURATE_PASS - - // 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 x3 = xd * xsq; - double r = asin_eval(xsq); - return static_cast(fputil::multiply_add(x3, r, xd)); - } - - // |x| > 1, return NaNs. - if (LIBC_UNLIKELY(x_abs > 0x3f80'0000U)) { - if (xbits.is_signaling_nan()) { - fputil::raise_except_if_required(FE_INVALID); - return FPBits::quiet_nan().get_val(); - } - - if (x_abs <= 0x7f80'0000U) { - fputil::set_errno_if_required(EDOM); - fputil::raise_except_if_required(FE_INVALID); - } - - return FPBits::quiet_nan().get_val(); - } - -#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS - // Check for exceptional values - if (auto r = ASINF_EXCEPTS_HI.lookup_odd(x_abs, x_sign); - LIBC_UNLIKELY(r.has_value())) - return r.value(); -#endif // !LIBC_MATH_HAS_SKIP_ACCURATE_PASS - - // 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( 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 * sqrt(u) * P(u), - - xbits.set_sign(Sign::POS); - double sign = SIGN[x_sign]; - double xd = static_cast(xbits.get_val()); - double u = fputil::multiply_add(-0.5, xd, 0.5); - 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 static_cast(fputil::multiply_add(c3, r, c2)); -} +LLVM_LIBC_FUNCTION(float, asinf, (float x)) { return math::asinf(x); } } // namespace LIBC_NAMESPACE_DECL diff --git a/libc/test/shared/CMakeLists.txt b/libc/test/shared/CMakeLists.txt index a38ba76a80325..e5dfac98b5e69 100644 --- a/libc/test/shared/CMakeLists.txt +++ b/libc/test/shared/CMakeLists.txt @@ -14,6 +14,7 @@ add_fp_unittest( libc.src.__support.math.acoshf16 libc.src.__support.math.acospif16 libc.src.__support.math.asin + libc.src.__support.math.asinf libc.src.__support.math.erff libc.src.__support.math.exp libc.src.__support.math.exp10 diff --git a/libc/test/shared/shared_math_test.cpp b/libc/test/shared/shared_math_test.cpp index 1b13fc0e0650f..7881d6841087b 100644 --- a/libc/test/shared/shared_math_test.cpp +++ b/libc/test/shared/shared_math_test.cpp @@ -40,6 +40,7 @@ TEST(LlvmLibcSharedMathTest, AllFloat) { EXPECT_FP_EQ(0x1.921fb6p+0, LIBC_NAMESPACE::shared::acosf(0.0f)); EXPECT_FP_EQ(0x0p+0f, LIBC_NAMESPACE::shared::acoshf(1.0f)); + EXPECT_FP_EQ(0x0p+0f, LIBC_NAMESPACE::shared::asinf(0.0f)); EXPECT_FP_EQ(0x0p+0f, LIBC_NAMESPACE::shared::erff(0.0f)); EXPECT_FP_EQ(0x1p+0f, LIBC_NAMESPACE::shared::exp10f(0.0f)); EXPECT_FP_EQ(0x1p+0f, LIBC_NAMESPACE::shared::expf(0.0f)); diff --git a/utils/bazel/llvm-project-overlay/libc/BUILD.bazel b/utils/bazel/llvm-project-overlay/libc/BUILD.bazel index 5ae999b8a9dcf..a9db35e9dba28 100644 --- a/utils/bazel/llvm-project-overlay/libc/BUILD.bazel +++ b/utils/bazel/llvm-project-overlay/libc/BUILD.bazel @@ -2230,6 +2230,22 @@ libc_support_library( ], ) +libc_support_library( + name = "__support_math_asinf", + hdrs = ["src/__support/math/asinf.h"], + deps = [ + ":__support_math_inv_trigf_utils", + ":__support_fputil_fenv_impl", + ":__support_fputil_fp_bits", + ":__support_fputil_except_value_utils", + ":__support_fputil_multiply_add", + ":__support_fputil_sqrt", + ":__support_macros_config", + ":__support_macros_optimization", + ":__support_macros_properties_cpu_features", + ], +) + libc_support_library( name = "__support_math_erff", hdrs = ["src/__support/math/erff.h"], @@ -2785,14 +2801,7 @@ libc_math_function( libc_math_function( name = "asinf", additional_deps = [ - ":__support_fputil_fma", - ":__support_fputil_multiply_add", - ":__support_fputil_nearest_integer", - ":__support_fputil_polyeval", - ":__support_fputil_sqrt", - ":__support_macros_optimization", - ":__support_macros_properties_cpu_features", - ":__support_math_inv_trigf_utils", + ":__support_math_asinf", ], )