diff --git a/libc/shared/math.h b/libc/shared/math.h index 70c6d375c22de..6f7805bc5e7ac 100644 --- a/libc/shared/math.h +++ b/libc/shared/math.h @@ -64,5 +64,6 @@ #include "math/ldexpf16.h" #include "math/rsqrtf.h" #include "math/rsqrtf16.h" +#include "math/sin.h" #endif // LLVM_LIBC_SHARED_MATH_H diff --git a/libc/shared/math/sin.h b/libc/shared/math/sin.h new file mode 100644 index 0000000000000..011cec6d8e2ed --- /dev/null +++ b/libc/shared/math/sin.h @@ -0,0 +1,23 @@ +//===-- Shared sin 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_SIN_H +#define LLVM_LIBC_SHARED_MATH_SIN_H + +#include "shared/libc_common.h" +#include "src/__support/math/sin.h" + +namespace LIBC_NAMESPACE_DECL { +namespace shared { + +using math::sin; + +} // namespace shared +} // namespace LIBC_NAMESPACE_DECL + +#endif // LLVM_LIBC_SHARED_MATH_SIN_H diff --git a/libc/src/__support/math/CMakeLists.txt b/libc/src/__support/math/CMakeLists.txt index f8622da52d983..e41cfec980ff9 100644 --- a/libc/src/__support/math/CMakeLists.txt +++ b/libc/src/__support/math/CMakeLists.txt @@ -1008,3 +1008,19 @@ add_header_library( .exp10f_utils libc.src.__support.FPUtil.multiply_add ) + +add_header_library( + sin + HDRS + sin.h + DEPENDS + libc.hdr.errno_macros + libc.src.errno.errno + 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.math.range_reduction_double + libc.src.__support.math.sincos_eval + libc.src.__support.macros.optimization +) diff --git a/libc/src/__support/math/sin.h b/libc/src/__support/math/sin.h new file mode 100644 index 0000000000000..af1eb59118c2d --- /dev/null +++ b/libc/src/__support/math/sin.h @@ -0,0 +1,181 @@ +//===-- Implementation header for sin ---------------------------*- 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 LIBC_SRC___SUPPORT_MATH_SIN_H +#define LIBC_SRC___SUPPORT_MATH_SIN_H + +#include "range_reduction_double_common.h" +#include "sincos_eval.h" +#include "src/__support/FPUtil/FEnvImpl.h" +#include "src/__support/FPUtil/FPBits.h" +#include "src/__support/FPUtil/double_double.h" +#include "src/__support/FPUtil/dyadic_float.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 + +#ifdef LIBC_TARGET_CPU_HAS_FMA_DOUBLE +#include "range_reduction_double_fma.h" +#else +#include "range_reduction_double_nofma.h" +#endif // LIBC_TARGET_CPU_HAS_FMA_DOUBLE + +namespace LIBC_NAMESPACE_DECL { + +namespace math { + +LIBC_INLINE static constexpr double sin(double x) { + using namespace math::range_reduction_double_internal; + using FPBits = typename fputil::FPBits; + FPBits xbits(x); + + uint16_t x_e = xbits.get_biased_exponent(); + + DoubleDouble y; + unsigned k = 0; + LargeRangeReduction range_reduction_large{}; + + // |x| < 2^16 + if (LIBC_LIKELY(x_e < FPBits::EXP_BIAS + FAST_PASS_EXPONENT)) { + // |x| < 2^-7 + if (LIBC_UNLIKELY(x_e < FPBits::EXP_BIAS - 7)) { + // |x| < 2^-26, |sin(x) - x| < ulp(x)/2. + if (LIBC_UNLIKELY(x_e < FPBits::EXP_BIAS - 26)) { + // Signed zeros. + if (LIBC_UNLIKELY(x == 0.0)) + return x + x; // Make sure it works with FTZ/DAZ. + +#ifdef LIBC_TARGET_CPU_HAS_FMA_DOUBLE + return fputil::multiply_add(x, -0x1.0p-54, x); +#else + if (LIBC_UNLIKELY(x_e < 4)) { + int rounding_mode = fputil::quick_get_round(); + if (rounding_mode == FE_TOWARDZERO || + (xbits.sign() == Sign::POS && rounding_mode == FE_DOWNWARD) || + (xbits.sign() == Sign::NEG && rounding_mode == FE_UPWARD)) + return FPBits(xbits.uintval() - 1).get_val(); + } + return fputil::multiply_add(x, -0x1.0p-54, x); +#endif // LIBC_TARGET_CPU_HAS_FMA_DOUBLE + } + // No range reduction needed. + k = 0; + y.lo = 0.0; + y.hi = x; + } else { + // Small range reduction. + k = range_reduction_small(x, y); + } + } else { + // Inf or NaN + if (LIBC_UNLIKELY(x_e > 2 * FPBits::EXP_BIAS)) { + // sin(+-Inf) = NaN + if (xbits.is_signaling_nan()) { + fputil::raise_except_if_required(FE_INVALID); + return FPBits::quiet_nan().get_val(); + } + + if (xbits.get_mantissa() == 0) { + fputil::set_errno_if_required(EDOM); + fputil::raise_except_if_required(FE_INVALID); + } + return x + FPBits::quiet_nan().get_val(); + } + + // Large range reduction. + k = range_reduction_large.fast(x, y); + } + + DoubleDouble sin_y, cos_y; + + [[maybe_unused]] double err = + math::sincos_eval_internal::sincos_eval(y, sin_y, cos_y); + + // Look up sin(k * pi/128) and cos(k * pi/128) +#ifdef LIBC_MATH_HAS_SMALL_TABLES + // Memory saving versions. Use 65-entry table. + auto get_idx_dd = [](unsigned kk) -> DoubleDouble { + unsigned idx = (kk & 64) ? 64 - (kk & 63) : (kk & 63); + DoubleDouble ans = SIN_K_PI_OVER_128[idx]; + if (kk & 128) { + ans.hi = -ans.hi; + ans.lo = -ans.lo; + } + return ans; + }; + DoubleDouble sin_k = get_idx_dd(k); + DoubleDouble cos_k = get_idx_dd(k + 64); +#else + // Fast look up version, but needs 256-entry table. + // cos(k * pi/128) = sin(k * pi/128 + pi/2) = sin((k + 64) * pi/128). + DoubleDouble sin_k = SIN_K_PI_OVER_128[k & 255]; + DoubleDouble cos_k = SIN_K_PI_OVER_128[(k + 64) & 255]; +#endif + + // After range reduction, k = round(x * 128 / pi) and y = x - k * (pi / 128). + // So k is an integer and -pi / 256 <= y <= pi / 256. + // Then sin(x) = sin((k * pi/128 + y) + // = sin(y) * cos(k*pi/128) + cos(y) * sin(k*pi/128) + DoubleDouble sin_k_cos_y = fputil::quick_mult(cos_y, sin_k); + DoubleDouble cos_k_sin_y = fputil::quick_mult(sin_y, cos_k); + + DoubleDouble rr = fputil::exact_add(sin_k_cos_y.hi, cos_k_sin_y.hi); + rr.lo += sin_k_cos_y.lo + cos_k_sin_y.lo; + +#ifdef LIBC_MATH_HAS_SKIP_ACCURATE_PASS + return rr.hi + rr.lo; +#else + // Accurate test and pass for correctly rounded implementation. + + double rlp = rr.lo + err; + double rlm = rr.lo - err; + + double r_upper = rr.hi + rlp; // (rr.lo + ERR); + double r_lower = rr.hi + rlm; // (rr.lo - ERR); + + // Ziv's rounding test. + if (LIBC_LIKELY(r_upper == r_lower)) + return r_upper; + + Float128 u_f128, sin_u, cos_u; + if (LIBC_LIKELY(x_e < FPBits::EXP_BIAS + FAST_PASS_EXPONENT)) + u_f128 = range_reduction_small_f128(x); + else + u_f128 = range_reduction_large.accurate(); + + math::sincos_eval_internal::sincos_eval(u_f128, sin_u, cos_u); + + auto get_sin_k = [](unsigned kk) -> Float128 { + unsigned idx = (kk & 64) ? 64 - (kk & 63) : (kk & 63); + Float128 ans = SIN_K_PI_OVER_128_F128[idx]; + if (kk & 128) + ans.sign = Sign::NEG; + return ans; + }; + + // cos(k * pi/128) = sin(k * pi/128 + pi/2) = sin((k + 64) * pi/128). + Float128 sin_k_f128 = get_sin_k(k); + Float128 cos_k_f128 = get_sin_k(k + 64); + + // sin(x) = sin(k * pi/128 + u) + // = sin(u) * cos(k*pi/128) + cos(u) * sin(k*pi/128) + Float128 r = fputil::quick_add(fputil::quick_mul(sin_k_f128, cos_u), + fputil::quick_mul(cos_k_f128, sin_u)); + + // TODO: Add assertion if Ziv's accuracy tests fail in debug mode. + // https://github.com/llvm/llvm-project/issues/96452. + + return static_cast(r); +#endif // !LIBC_MATH_HAS_SKIP_ACCURATE_PASS +} + +} // namespace math + +} // namespace LIBC_NAMESPACE_DECL + +#endif // LIBC_SRC___SUPPORT_MATH_SIN_H diff --git a/libc/src/math/generic/CMakeLists.txt b/libc/src/math/generic/CMakeLists.txt index e343c1f15c3f1..97c41c47f70f1 100644 --- a/libc/src/math/generic/CMakeLists.txt +++ b/libc/src/math/generic/CMakeLists.txt @@ -363,17 +363,7 @@ add_entrypoint_object( HDRS ../sin.h DEPENDS - libc.src.__support.math.range_reduction_double - libc.src.__support.math.sincos_eval - libc.hdr.errno_macros - libc.src.errno.errno - 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.rounding_mode - libc.src.__support.macros.optimization + libc.src.__support.math.sin ) add_entrypoint_object( diff --git a/libc/src/math/generic/sin.cpp b/libc/src/math/generic/sin.cpp index 1b6310f99366c..6d758d8450bc2 100644 --- a/libc/src/math/generic/sin.cpp +++ b/libc/src/math/generic/sin.cpp @@ -7,174 +7,10 @@ //===----------------------------------------------------------------------===// #include "src/math/sin.h" -#include "hdr/errno_macros.h" -#include "src/__support/FPUtil/FEnvImpl.h" -#include "src/__support/FPUtil/FPBits.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/rounding_mode.h" -#include "src/__support/common.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/range_reduction_double_common.h" -#include "src/__support/math/sincos_eval.h" - -#ifdef LIBC_TARGET_CPU_HAS_FMA_DOUBLE -#include "src/__support/math/range_reduction_double_fma.h" -#else -#include "src/__support/math/range_reduction_double_nofma.h" -#endif // LIBC_TARGET_CPU_HAS_FMA_DOUBLE +#include "src/__support/math/sin.h" namespace LIBC_NAMESPACE_DECL { -using DoubleDouble = fputil::DoubleDouble; -using Float128 = typename fputil::DyadicFloat<128>; - -LLVM_LIBC_FUNCTION(double, sin, (double x)) { - using namespace math::range_reduction_double_internal; - using FPBits = typename fputil::FPBits; - FPBits xbits(x); - - uint16_t x_e = xbits.get_biased_exponent(); - - DoubleDouble y; - unsigned k; - LargeRangeReduction range_reduction_large{}; - - // |x| < 2^16 - if (LIBC_LIKELY(x_e < FPBits::EXP_BIAS + FAST_PASS_EXPONENT)) { - // |x| < 2^-7 - if (LIBC_UNLIKELY(x_e < FPBits::EXP_BIAS - 7)) { - // |x| < 2^-26, |sin(x) - x| < ulp(x)/2. - if (LIBC_UNLIKELY(x_e < FPBits::EXP_BIAS - 26)) { - // Signed zeros. - if (LIBC_UNLIKELY(x == 0.0)) - return x + x; // Make sure it works with FTZ/DAZ. - -#ifdef LIBC_TARGET_CPU_HAS_FMA_DOUBLE - return fputil::multiply_add(x, -0x1.0p-54, x); -#else - if (LIBC_UNLIKELY(x_e < 4)) { - int rounding_mode = fputil::quick_get_round(); - if (rounding_mode == FE_TOWARDZERO || - (xbits.sign() == Sign::POS && rounding_mode == FE_DOWNWARD) || - (xbits.sign() == Sign::NEG && rounding_mode == FE_UPWARD)) - return FPBits(xbits.uintval() - 1).get_val(); - } - return fputil::multiply_add(x, -0x1.0p-54, x); -#endif // LIBC_TARGET_CPU_HAS_FMA_DOUBLE - } - // No range reduction needed. - k = 0; - y.lo = 0.0; - y.hi = x; - } else { - // Small range reduction. - k = range_reduction_small(x, y); - } - } else { - // Inf or NaN - if (LIBC_UNLIKELY(x_e > 2 * FPBits::EXP_BIAS)) { - // sin(+-Inf) = NaN - if (xbits.is_signaling_nan()) { - fputil::raise_except_if_required(FE_INVALID); - return FPBits::quiet_nan().get_val(); - } - - if (xbits.get_mantissa() == 0) { - fputil::set_errno_if_required(EDOM); - fputil::raise_except_if_required(FE_INVALID); - } - return x + FPBits::quiet_nan().get_val(); - } - - // Large range reduction. - k = range_reduction_large.fast(x, y); - } - - DoubleDouble sin_y, cos_y; - - [[maybe_unused]] double err = - math::sincos_eval_internal::sincos_eval(y, sin_y, cos_y); - - // Look up sin(k * pi/128) and cos(k * pi/128) -#ifdef LIBC_MATH_HAS_SMALL_TABLES - // Memory saving versions. Use 65-entry table. - auto get_idx_dd = [](unsigned kk) -> DoubleDouble { - unsigned idx = (kk & 64) ? 64 - (kk & 63) : (kk & 63); - DoubleDouble ans = SIN_K_PI_OVER_128[idx]; - if (kk & 128) { - ans.hi = -ans.hi; - ans.lo = -ans.lo; - } - return ans; - }; - DoubleDouble sin_k = get_idx_dd(k); - DoubleDouble cos_k = get_idx_dd(k + 64); -#else - // Fast look up version, but needs 256-entry table. - // cos(k * pi/128) = sin(k * pi/128 + pi/2) = sin((k + 64) * pi/128). - DoubleDouble sin_k = SIN_K_PI_OVER_128[k & 255]; - DoubleDouble cos_k = SIN_K_PI_OVER_128[(k + 64) & 255]; -#endif - - // After range reduction, k = round(x * 128 / pi) and y = x - k * (pi / 128). - // So k is an integer and -pi / 256 <= y <= pi / 256. - // Then sin(x) = sin((k * pi/128 + y) - // = sin(y) * cos(k*pi/128) + cos(y) * sin(k*pi/128) - DoubleDouble sin_k_cos_y = fputil::quick_mult(cos_y, sin_k); - DoubleDouble cos_k_sin_y = fputil::quick_mult(sin_y, cos_k); - - DoubleDouble rr = fputil::exact_add(sin_k_cos_y.hi, cos_k_sin_y.hi); - rr.lo += sin_k_cos_y.lo + cos_k_sin_y.lo; - -#ifdef LIBC_MATH_HAS_SKIP_ACCURATE_PASS - return rr.hi + rr.lo; -#else - // Accurate test and pass for correctly rounded implementation. - - double rlp = rr.lo + err; - double rlm = rr.lo - err; - - double r_upper = rr.hi + rlp; // (rr.lo + ERR); - double r_lower = rr.hi + rlm; // (rr.lo - ERR); - - // Ziv's rounding test. - if (LIBC_LIKELY(r_upper == r_lower)) - return r_upper; - - Float128 u_f128, sin_u, cos_u; - if (LIBC_LIKELY(x_e < FPBits::EXP_BIAS + FAST_PASS_EXPONENT)) - u_f128 = range_reduction_small_f128(x); - else - u_f128 = range_reduction_large.accurate(); - - math::sincos_eval_internal::sincos_eval(u_f128, sin_u, cos_u); - - auto get_sin_k = [](unsigned kk) -> Float128 { - unsigned idx = (kk & 64) ? 64 - (kk & 63) : (kk & 63); - Float128 ans = SIN_K_PI_OVER_128_F128[idx]; - if (kk & 128) - ans.sign = Sign::NEG; - return ans; - }; - - // cos(k * pi/128) = sin(k * pi/128 + pi/2) = sin((k + 64) * pi/128). - Float128 sin_k_f128 = get_sin_k(k); - Float128 cos_k_f128 = get_sin_k(k + 64); - - // sin(x) = sin(k * pi/128 + u) - // = sin(u) * cos(k*pi/128) + cos(u) * sin(k*pi/128) - Float128 r = fputil::quick_add(fputil::quick_mul(sin_k_f128, cos_u), - fputil::quick_mul(cos_k_f128, sin_u)); - - // TODO: Add assertion if Ziv's accuracy tests fail in debug mode. - // https://github.com/llvm/llvm-project/issues/96452. - - return static_cast(r); -#endif // !LIBC_MATH_HAS_SKIP_ACCURATE_PASS -} +LLVM_LIBC_FUNCTION(double, sin, (double x)) { return math::sin(x); } } // namespace LIBC_NAMESPACE_DECL diff --git a/libc/test/shared/CMakeLists.txt b/libc/test/shared/CMakeLists.txt index 4482d6344ae03..70849de9fb59c 100644 --- a/libc/test/shared/CMakeLists.txt +++ b/libc/test/shared/CMakeLists.txt @@ -60,4 +60,5 @@ add_fp_unittest( libc.src.__support.math.ldexpf16 libc.src.__support.math.rsqrtf libc.src.__support.math.rsqrtf16 + libc.src.__support.math.sin ) diff --git a/libc/test/shared/shared_math_test.cpp b/libc/test/shared/shared_math_test.cpp index f09b0dcd4ef9f..fc05e535f9266 100644 --- a/libc/test/shared/shared_math_test.cpp +++ b/libc/test/shared/shared_math_test.cpp @@ -89,6 +89,7 @@ TEST(LlvmLibcSharedMathTest, AllDouble) { EXPECT_FP_EQ(0x1p+0, LIBC_NAMESPACE::shared::exp2(0.0)); EXPECT_FP_EQ(0x1p+0, LIBC_NAMESPACE::shared::exp10(0.0)); EXPECT_FP_EQ(0x0p+0, LIBC_NAMESPACE::shared::expm1(0.0)); + EXPECT_FP_EQ(0.0, LIBC_NAMESPACE::shared::sin(0.0)); } #ifdef LIBC_TYPES_HAS_FLOAT128 diff --git a/utils/bazel/llvm-project-overlay/libc/BUILD.bazel b/utils/bazel/llvm-project-overlay/libc/BUILD.bazel index 913e913572b14..ac41bd359f0d4 100644 --- a/utils/bazel/llvm-project-overlay/libc/BUILD.bazel +++ b/utils/bazel/llvm-project-overlay/libc/BUILD.bazel @@ -3176,6 +3176,18 @@ libc_support_library( ], ) +libc_support_library( + name = "__support_math_sin", + hdrs = ["src/__support/math/sin.h"], + deps = [ + ":__support_fputil_multiply_add", + ":__support_macros_optimization", + ":__support_macros_properties_cpu_features", + ":__support_range_reduction_double", + ":__support_sincos_eval", + ], +) + libc_support_library( name = "__support_math_sinhfcoshf_utils", hdrs = ["src/__support/math/sinhfcoshf_utils.h"], @@ -4732,11 +4744,7 @@ libc_math_function(name = "setpayloadsigf16") libc_math_function( name = "sin", additional_deps = [ - ":__support_fputil_multiply_add", - ":__support_macros_optimization", - ":__support_macros_properties_cpu_features", - ":__support_range_reduction_double", - ":__support_sincos_eval", + ":__support_math_sin", ], )