From ea93c538c78ed191f83026d372bd71dd4a135cbc Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Hendrik=20H=C3=BCbner?= <117831077+HendrikHuebner@users.noreply.github.com> Date: Mon, 1 Jul 2024 22:38:03 +0200 Subject: [PATCH] [libc][math][c23] Implemented sinpif function correctly rounded for all rounding modes. (#97149) This implements the sinpif function. An exhaustive test shows it's correct for all rounding modes. Issue: #94895 --- libc/config/darwin/arm/entrypoints.txt | 1 + libc/config/linux/aarch64/entrypoints.txt | 1 + libc/config/linux/riscv/entrypoints.txt | 1 + libc/config/linux/x86_64/entrypoints.txt | 1 + libc/docs/math/index.rst | 2 +- libc/src/math/CMakeLists.txt | 1 + libc/src/math/generic/CMakeLists.txt | 19 +++ libc/src/math/generic/range_reduction_fma.h | 1 - libc/src/math/generic/sincosf_utils.h | 48 +++++-- libc/src/math/generic/sinpif.cpp | 111 +++++++++++++++ libc/src/math/sinpif.h | 18 +++ libc/test/src/math/CMakeLists.txt | 16 +++ libc/test/src/math/exhaustive/CMakeLists.txt | 16 +++ libc/test/src/math/exhaustive/sinpif_test.cpp | 35 +++++ libc/test/src/math/sinpif_test.cpp | 126 ++++++++++++++++++ libc/test/src/math/smoke/CMakeLists.txt | 13 ++ libc/test/src/math/smoke/sinpif_test.cpp | 43 ++++++ libc/utils/MPFRWrapper/MPFRUtils.cpp | 22 ++- libc/utils/MPFRWrapper/MPFRUtils.h | 1 + 19 files changed, 459 insertions(+), 17 deletions(-) create mode 100644 libc/src/math/generic/sinpif.cpp create mode 100644 libc/src/math/sinpif.h create mode 100644 libc/test/src/math/exhaustive/sinpif_test.cpp create mode 100644 libc/test/src/math/sinpif_test.cpp create mode 100644 libc/test/src/math/smoke/sinpif_test.cpp diff --git a/libc/config/darwin/arm/entrypoints.txt b/libc/config/darwin/arm/entrypoints.txt index 4d5536318724..cb4603c79c79 100644 --- a/libc/config/darwin/arm/entrypoints.txt +++ b/libc/config/darwin/arm/entrypoints.txt @@ -230,6 +230,7 @@ set(TARGET_LIBM_ENTRYPOINTS libc.src.math.sinhf libc.src.math.sin libc.src.math.sinf + libc.src.math.sinpif libc.src.math.sqrt libc.src.math.sqrtf libc.src.math.sqrtl diff --git a/libc/config/linux/aarch64/entrypoints.txt b/libc/config/linux/aarch64/entrypoints.txt index ea89f8bd138d..07fdbfdbcd7c 100644 --- a/libc/config/linux/aarch64/entrypoints.txt +++ b/libc/config/linux/aarch64/entrypoints.txt @@ -485,6 +485,7 @@ set(TARGET_LIBM_ENTRYPOINTS libc.src.math.sincosf libc.src.math.sinf libc.src.math.sinhf + libc.src.math.sinpif libc.src.math.sqrt libc.src.math.sqrtf libc.src.math.sqrtl diff --git a/libc/config/linux/riscv/entrypoints.txt b/libc/config/linux/riscv/entrypoints.txt index ee8b3d531637..51d85eed9ff1 100644 --- a/libc/config/linux/riscv/entrypoints.txt +++ b/libc/config/linux/riscv/entrypoints.txt @@ -493,6 +493,7 @@ set(TARGET_LIBM_ENTRYPOINTS libc.src.math.sincosf libc.src.math.sinf libc.src.math.sinhf + libc.src.math.sinpif libc.src.math.sqrt libc.src.math.sqrtf libc.src.math.sqrtl diff --git a/libc/config/linux/x86_64/entrypoints.txt b/libc/config/linux/x86_64/entrypoints.txt index 9f67dea00669..89f9011a52b6 100644 --- a/libc/config/linux/x86_64/entrypoints.txt +++ b/libc/config/linux/x86_64/entrypoints.txt @@ -514,6 +514,7 @@ set(TARGET_LIBM_ENTRYPOINTS libc.src.math.sincosf libc.src.math.sinf libc.src.math.sinhf + libc.src.math.sinpif libc.src.math.sqrt libc.src.math.sqrtf libc.src.math.sqrtl diff --git a/libc/docs/math/index.rst b/libc/docs/math/index.rst index 3ca7bf9fdb28..ccafb1f9c73f 100644 --- a/libc/docs/math/index.rst +++ b/libc/docs/math/index.rst @@ -330,7 +330,7 @@ Higher Math Functions +-----------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+----------------------------+ | sinh | |check| | | | | | 7.12.5.5 | F.10.2.5 | +-----------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+----------------------------+ -| sinpi | | | | | | 7.12.4.13 | F.10.1.13 | +| sinpi | |check| | | | | | 7.12.4.13 | F.10.1.13 | +-----------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+----------------------------+ | sqrt | |check| | |check| | |check| | | |check| | 7.12.7.10 | F.10.4.10 | +-----------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+----------------------------+ diff --git a/libc/src/math/CMakeLists.txt b/libc/src/math/CMakeLists.txt index 4777f669f0cb..607051a6ad78 100644 --- a/libc/src/math/CMakeLists.txt +++ b/libc/src/math/CMakeLists.txt @@ -376,6 +376,7 @@ add_math_entrypoint_object(sincosf) add_math_entrypoint_object(sin) add_math_entrypoint_object(sinf) +add_math_entrypoint_object(sinpif) add_math_entrypoint_object(sinh) add_math_entrypoint_object(sinhf) diff --git a/libc/src/math/generic/CMakeLists.txt b/libc/src/math/generic/CMakeLists.txt index 7dd6c4812fa1..395d4f8cc927 100644 --- a/libc/src/math/generic/CMakeLists.txt +++ b/libc/src/math/generic/CMakeLists.txt @@ -283,6 +283,25 @@ add_entrypoint_object( -O3 ) +add_entrypoint_object( + sinpif + SRCS + sinpif.cpp + HDRS + ../sinpif.h + DEPENDS + .sincosf_utils + libc.src.__support.FPUtil.fenv_impl + libc.src.__support.FPUtil.fp_bits + libc.src.__support.FPUtil.fma + libc.src.__support.FPUtil.multiply_add + libc.src.__support.FPUtil.polyeval + libc.src.__support.common + libc.src.__support.macros.optimization + COMPILE_OPTIONS + -O3 +) + add_entrypoint_object( sincosf SRCS diff --git a/libc/src/math/generic/range_reduction_fma.h b/libc/src/math/generic/range_reduction_fma.h index 82b4ae1c705e..97b9e3d0f89d 100644 --- a/libc/src/math/generic/range_reduction_fma.h +++ b/libc/src/math/generic/range_reduction_fma.h @@ -12,7 +12,6 @@ #include "src/__support/FPUtil/FMA.h" #include "src/__support/FPUtil/FPBits.h" #include "src/__support/FPUtil/nearest_integer.h" -#include "src/__support/common.h" namespace LIBC_NAMESPACE { diff --git a/libc/src/math/generic/sincosf_utils.h b/libc/src/math/generic/sincosf_utils.h index 904df4f4ed08..f20fb6a05a32 100644 --- a/libc/src/math/generic/sincosf_utils.h +++ b/libc/src/math/generic/sincosf_utils.h @@ -19,6 +19,7 @@ using LIBC_NAMESPACE::fma::FAST_PASS_BOUND; using LIBC_NAMESPACE::fma::large_range_reduction; using LIBC_NAMESPACE::fma::small_range_reduction; + #else #include "range_reduction.h" // using namespace LIBC_NAMESPACE::generic; @@ -58,18 +59,9 @@ const double SIN_K_PI_OVER_32[64] = { -0x1.917a6bc29b42cp-4, }; -LIBC_INLINE void sincosf_eval(double xd, uint32_t x_abs, double &sin_k, - double &cos_k, double &sin_y, double &cosm1_y) { - int64_t k; - double y; - - if (LIBC_LIKELY(x_abs < FAST_PASS_BOUND)) { - k = small_range_reduction(xd, y); - } else { - fputil::FPBits x_bits(x_abs); - k = large_range_reduction(xd, x_bits.get_exponent(), y); - } - +static LIBC_INLINE void sincosf_poly_eval(int64_t k, double y, double &sin_k, + double &cos_k, double &sin_y, + double &cosm1_y) { // After range reduction, k = round(x * 32 / pi) and y = (x * 32 / pi) - k. // So k is an integer and -0.5 <= y <= 0.5. // Then sin(x) = sin((k + y)*pi/32) @@ -95,6 +87,38 @@ LIBC_INLINE void sincosf_eval(double xd, uint32_t x_abs, double &sin_k, 0x1.03c1f070c2e27p-18, -0x1.55cc84bd942p-30); } +LIBC_INLINE void sincosf_eval(double xd, uint32_t x_abs, double &sin_k, + double &cos_k, double &sin_y, double &cosm1_y) { + int64_t k; + double y; + + if (LIBC_LIKELY(x_abs < FAST_PASS_BOUND)) { + k = small_range_reduction(xd, y); + } else { + fputil::FPBits x_bits(x_abs); + k = large_range_reduction(xd, x_bits.get_exponent(), y); + } + + sincosf_poly_eval(k, y, sin_k, cos_k, sin_y, cosm1_y); +} + +// Return k and y, where +// k = round(x * 32) and y = (x * 32) - k. +// => pi * x = (k + y) * pi / 32 +static LIBC_INLINE int64_t range_reduction_sincospi(double x, double &y) { + double kd = fputil::nearest_integer(x * 32); + y = fputil::multiply_add(x, 32.0, -kd); + + return static_cast(kd); +} + +LIBC_INLINE void sincospif_eval(double xd, double &sin_k, double &cos_k, + double &sin_y, double &cosm1_y) { + double y; + int64_t k = range_reduction_sincospi(xd, y); + sincosf_poly_eval(k, y, sin_k, cos_k, sin_y, cosm1_y); +} + } // namespace LIBC_NAMESPACE #endif // LLVM_LIBC_SRC_MATH_GENERIC_SINCOSF_UTILS_H diff --git a/libc/src/math/generic/sinpif.cpp b/libc/src/math/generic/sinpif.cpp new file mode 100644 index 000000000000..662263c9fc43 --- /dev/null +++ b/libc/src/math/generic/sinpif.cpp @@ -0,0 +1,111 @@ +//===-- Single-precision sinpif 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/sinpif.h" +#include "sincosf_utils.h" +#include "src/__support/FPUtil/FEnvImpl.h" +#include "src/__support/FPUtil/FPBits.h" +#include "src/__support/FPUtil/PolyEval.h" +#include "src/__support/FPUtil/multiply_add.h" +#include "src/__support/common.h" +#include "src/__support/macros/optimization.h" // LIBC_UNLIKELY + +namespace LIBC_NAMESPACE { + +LLVM_LIBC_FUNCTION(float, sinpif, (float x)) { + using FPBits = typename fputil::FPBits; + FPBits xbits(x); + + uint32_t x_u = xbits.uintval(); + uint32_t x_abs = x_u & 0x7fff'ffffU; + double xd = static_cast(x); + + // Range reduction: + // For |x| > pi/32, we perform range reduction as follows: + // Find k and y such that: + // x = (k + y) * 1/32 + // k is an integer + // |y| < 0.5 + // For small range (|x| < 2^45 when FMA instructions are available, 2^22 + // otherwise), this is done by performing: + // k = round(x * 32) + // y = x * 32 - k + // + // Once k and y are computed, we then deduce the answer by the sine of sum + // formula: + // sin(x * pi) = sin((k + y)*pi/32) + // = sin(y*pi/32) * cos(k*pi/32) + cos(y*pi/32) * sin(k*pi/32) + // The values of sin(k*pi/32) and cos(k*pi/32) for k = 0..31 are precomputed + // and stored using a vector of 32 doubles. Sin(y*pi/32) and cos(y*pi/32) are + // computed using degree-7 and degree-6 minimax polynomials generated by + // Sollya respectively. + + // |x| <= 1/16 + if (LIBC_UNLIKELY(x_abs <= 0x3d80'0000U)) { + + if (LIBC_UNLIKELY(x_abs < 0x33CD'01D7U)) { + if (LIBC_UNLIKELY(x_abs == 0U)) { + // For signed zeros. + return x; + } + + // For very small values we can approximate sinpi(x) with x * pi + // An exhaustive test shows that this is accurate for |x| < 9.546391 × + // 10-8 + double xdpi = xd * 0x1.921fb54442d18p1; + return static_cast(xdpi); + } + + // |x| < 1/16. + double xsq = xd * xd; + + // Degree-9 polynomial approximation: + // sinpi(x) ~ x + a_3 x^3 + a_5 x^5 + a_7 x^7 + a_9 x^9 + // = x (1 + a_3 x^2 + ... + a_9 x^8) + // = x * P(x^2) + // generated by Sollya with the following commands: + // > display = hexadecimal; + // > Q = fpminimax(sin(pi * x)/x, [|0, 2, 4, 6, 8|], [|D...|], [0, 1/16]); + double result = fputil::polyeval( + xsq, 0x1.921fb54442d18p1, -0x1.4abbce625bbf2p2, 0x1.466bc675e116ap1, + -0x1.32d2c0b62d41cp-1, 0x1.501ec4497cb7dp-4); + return static_cast(xd * result); + } + + // Numbers greater or equal to 2^23 are always integers or NaN + if (LIBC_UNLIKELY(x_abs >= 0x4B00'0000)) { + + // check for NaN values + if (LIBC_UNLIKELY(x_abs >= 0x7f80'0000U)) { + if (x_abs == 0x7f80'0000U) { + fputil::set_errno_if_required(EDOM); + fputil::raise_except_if_required(FE_INVALID); + } + + return x + FPBits::quiet_nan().get_val(); + } + + return FPBits::zero(xbits.sign()).get_val(); + } + + // Combine the results with the sine of sum formula: + // sin(x * pi) = sin((k + y)*pi/32) + // = sin(y*pi/32) * cos(k*pi/32) + cos(y*pi/32) * sin(k*pi/32) + // = sin_y * cos_k + (1 + cosm1_y) * sin_k + // = sin_y * cos_k + (cosm1_y * sin_k + sin_k) + double sin_k, cos_k, sin_y, cosm1_y; + sincospif_eval(xd, sin_k, cos_k, sin_y, cosm1_y); + + if (LIBC_UNLIKELY(sin_y == 0 && sin_k == 0)) + return FPBits::zero(xbits.sign()).get_val(); + + return static_cast(fputil::multiply_add( + sin_y, cos_k, fputil::multiply_add(cosm1_y, sin_k, sin_k))); +} + +} // namespace LIBC_NAMESPACE diff --git a/libc/src/math/sinpif.h b/libc/src/math/sinpif.h new file mode 100644 index 000000000000..63d481945f47 --- /dev/null +++ b/libc/src/math/sinpif.h @@ -0,0 +1,18 @@ +//===-- Implementation header for sinpif ------------------------*- 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_SINPIF_H +#define LLVM_LIBC_SRC_MATH_SINPIF_H + +namespace LIBC_NAMESPACE { + +float sinpif(float x); + +} // namespace LIBC_NAMESPACE + +#endif // LLVM_LIBC_SRC_MATH_SINPIF_H diff --git a/libc/test/src/math/CMakeLists.txt b/libc/test/src/math/CMakeLists.txt index 263bee65315f..c07c6d77fa23 100644 --- a/libc/test/src/math/CMakeLists.txt +++ b/libc/test/src/math/CMakeLists.txt @@ -44,6 +44,22 @@ add_fp_unittest( libc.src.__support.FPUtil.fp_bits ) +add_fp_unittest( + sinpif_test + NEED_MPFR + SUITE + libc-math-unittests + SRCS + sinpif_test.cpp + HDRS + sdcomp26094.h + DEPENDS + libc.src.errno.errno + libc.src.math.sinpif + libc.src.__support.CPP.array + libc.src.__support.FPUtil.fp_bits +) + add_fp_unittest( sin_test NEED_MPFR diff --git a/libc/test/src/math/exhaustive/CMakeLists.txt b/libc/test/src/math/exhaustive/CMakeLists.txt index fb3596c3378f..412ca031d0e9 100644 --- a/libc/test/src/math/exhaustive/CMakeLists.txt +++ b/libc/test/src/math/exhaustive/CMakeLists.txt @@ -42,6 +42,22 @@ add_fp_unittest( -lpthread ) +add_fp_unittest( + sinpif_test + NO_RUN_POSTBUILD + NEED_MPFR + SUITE + libc_math_exhaustive_tests + SRCS + sinpif_test.cpp + DEPENDS + .exhaustive_test + libc.src.math.sinpif + libc.src.__support.FPUtil.fp_bits + LINK_LIBRARIES + -lpthread +) + add_fp_unittest( cosf_test NO_RUN_POSTBUILD diff --git a/libc/test/src/math/exhaustive/sinpif_test.cpp b/libc/test/src/math/exhaustive/sinpif_test.cpp new file mode 100644 index 000000000000..8bc1d81eb7e3 --- /dev/null +++ b/libc/test/src/math/exhaustive/sinpif_test.cpp @@ -0,0 +1,35 @@ +//===-- Exhaustive test for sinpif ----------------------------------------===// +// +// 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 "mpfr.h" +#include "src/math/sinpif.h" +#include "utils/MPFRWrapper/MPFRUtils.h" +#include + +namespace mpfr = LIBC_NAMESPACE::testing::mpfr; + +using LlvmLibcSinpifExhaustiveTest = + LlvmLibcUnaryOpExhaustiveMathTest; + +static constexpr uint32_t POS_START = 0x0000'0000U; +static constexpr uint32_t POS_STOP = 0x7f80'0000U; + +// Range: [0, Inf] +TEST_F(LlvmLibcSinpifExhaustiveTest, PostiveRange) { + test_full_range_all_roundings(POS_START, POS_STOP); +} + +// Range: [-Inf, 0] +static constexpr uint32_t NEG_START = 0xb000'0000U; +static constexpr uint32_t NEG_STOP = 0xff80'0000U; + +TEST_F(LlvmLibcSinpifExhaustiveTest, NegativeRange) { + test_full_range_all_roundings(NEG_START, NEG_STOP); +} diff --git a/libc/test/src/math/sinpif_test.cpp b/libc/test/src/math/sinpif_test.cpp new file mode 100644 index 000000000000..d00fd77d288c --- /dev/null +++ b/libc/test/src/math/sinpif_test.cpp @@ -0,0 +1,126 @@ +//===-- Unittests for sinpif ----------------------------------------------===// +// +// 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/errno/libc_errno.h" +#include "src/math/sinpif.h" +#include "test/UnitTest/FPMatcher.h" +#include "test/src/math/sdcomp26094.h" +#include "utils/MPFRWrapper/MPFRUtils.h" + +#include + +using LlvmLibcSinpifTest = LIBC_NAMESPACE::testing::FPTest; + +using LIBC_NAMESPACE::testing::SDCOMP26094_VALUES; + +namespace mpfr = LIBC_NAMESPACE::testing::mpfr; + +TEST_F(LlvmLibcSinpifTest, SpecialNumbers) { + LIBC_NAMESPACE::libc_errno = 0; + + EXPECT_FP_EQ(aNaN, LIBC_NAMESPACE::sinpif(aNaN)); + EXPECT_MATH_ERRNO(0); + + EXPECT_FP_EQ(0.0f, LIBC_NAMESPACE::sinpif(0.0f)); + EXPECT_MATH_ERRNO(0); + + EXPECT_FP_EQ(-0.0f, LIBC_NAMESPACE::sinpif(-0.0f)); + EXPECT_MATH_ERRNO(0); + + EXPECT_FP_EQ(aNaN, LIBC_NAMESPACE::sinpif(inf)); + EXPECT_MATH_ERRNO(EDOM); + + EXPECT_FP_EQ(aNaN, LIBC_NAMESPACE::sinpif(neg_inf)); + EXPECT_MATH_ERRNO(EDOM); +} + +TEST_F(LlvmLibcSinpifTest, SpecificBitPatterns) { + constexpr int N = 36; + constexpr uint32_t INPUTS[N] = { + 0x3f06'0a92U, // x = pi/6 + 0x3f3a'dc51U, // x = 0x1.75b8a2p-1f + 0x3f49'0fdbU, // x = pi/4 + 0x3f86'0a92U, // x = pi/3 + 0x3fa7'832aU, // x = 0x1.4f0654p+0f + 0x3fc9'0fdbU, // x = pi/2 + 0x4017'1973U, // x = 0x1.2e32e6p+1f + 0x4049'0fdbU, // x = pi + 0x4096'cbe4U, // x = 0x1.2d97c8p+2f + 0x40c9'0fdbU, // x = 2*pi + 0x433b'7490U, // x = 0x1.76e92p+7f + 0x437c'e5f1U, // x = 0x1.f9cbe2p+7f + 0x4619'9998U, // x = 0x1.33333p+13f + 0x474d'246fU, // x = 0x1.9a48dep+15f + 0x4afd'ece4U, // x = 0x1.fbd9c8p+22f + 0x4c23'32e9U, // x = 0x1.4665d2p+25f + 0x50a3'e87fU, // x = 0x1.47d0fep+34f + 0x5239'47f6U, // x = 0x1.728fecp+37f + 0x53b1'46a6U, // x = 0x1.628d4cp+40f + 0x55ca'fb2aU, // x = 0x1.95f654p+44f + 0x588e'f060U, // x = 0x1.1de0cp+50f + 0x5c07'bcd0U, // x = 0x1.0f79ap+57f + 0x5ebc'fddeU, // x = 0x1.79fbbcp+62f + 0x5fa6'eba7U, // x = 0x1.4dd74ep+64f + 0x61a4'0b40U, // x = 0x1.48168p+68f + 0x6386'134eU, // x = 0x1.0c269cp+72f + 0x6589'8498U, // x = 0x1.13093p+76f + 0x6600'0001U, // x = 0x1.000002p+77f + 0x664e'46e4U, // x = 0x1.9c8dc8p+77f + 0x66b0'14aaU, // x = 0x1.602954p+78f + 0x67a9'242bU, // x = 0x1.524856p+80f + 0x6a19'76f1U, // x = 0x1.32ede2p+85f + 0x6c55'da58U, // x = 0x1.abb4bp+89f + 0x6f79'be45U, // x = 0x1.f37c8ap+95f + 0x7276'69d4U, // x = 0x1.ecd3a8p+101f + 0x7758'4625U, // x = 0x1.b08c4ap+111f + }; + + for (int i = 0; i < N; ++i) { + float x = FPBits(INPUTS[i]).get_val(); + EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Sinpi, x, + LIBC_NAMESPACE::sinpif(x), 0.5); + EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Sinpi, -x, + LIBC_NAMESPACE::sinpif(-x), 0.5); + } +} + +// For small values, sinpi(x) is pi * x. +TEST_F(LlvmLibcSinpifTest, SmallValues) { + float x = FPBits(0x1780'0000U).get_val(); + EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Sinpi, x, + LIBC_NAMESPACE::sinpif(x), 0.5); + + x = FPBits(0x0040'0000U).get_val(); + EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Sinpi, x, + LIBC_NAMESPACE::sinpif(x), 0.5); +} + +// SDCOMP-26094: check sinfpi in the cases for which the range reducer +// returns values furthest beyond its nominal upper bound of pi/4. +TEST_F(LlvmLibcSinpifTest, SDCOMP_26094) { + for (uint32_t v : SDCOMP26094_VALUES) { + float x = FPBits((v)).get_val(); + EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Sinpi, x, + LIBC_NAMESPACE::sinpif(x), 0.5); + } +} + +// sinpi(-n) = -0.0 +// sinpi(+n) = +0.0 +TEST_F(LlvmLibcSinpifTest, SignedZeros) { + EXPECT_FP_EQ(-0.0, LIBC_NAMESPACE::sinpif(-0x99)); + EXPECT_FP_EQ(-0.0, LIBC_NAMESPACE::sinpif(-0x1p+43)); + EXPECT_FP_EQ(-0.0, LIBC_NAMESPACE::sinpif(-0x1.4p+64)); + EXPECT_FP_EQ(-0.0, LIBC_NAMESPACE::sinpif(-0x1.cp+106)); + EXPECT_FP_EQ(-0.0, LIBC_NAMESPACE::sinpif(-0x1.cp+21)); + EXPECT_FP_EQ(-0.0, LIBC_NAMESPACE::sinpif(-0x9999)); + EXPECT_FP_EQ(0.0, LIBC_NAMESPACE::sinpif(0x1p+43)); + EXPECT_FP_EQ(0.0, LIBC_NAMESPACE::sinpif(0x1.4p+64)); + EXPECT_FP_EQ(0.0, LIBC_NAMESPACE::sinpif(0x1.cp+106)); + EXPECT_FP_EQ(0.0, LIBC_NAMESPACE::sinpif(0x1.cp+21)); +} diff --git a/libc/test/src/math/smoke/CMakeLists.txt b/libc/test/src/math/smoke/CMakeLists.txt index a1b4a03eb2e8..a363644e6fa8 100644 --- a/libc/test/src/math/smoke/CMakeLists.txt +++ b/libc/test/src/math/smoke/CMakeLists.txt @@ -27,6 +27,19 @@ add_fp_unittest( libc.src.__support.FPUtil.fp_bits ) +add_fp_unittest( + sinpif_test + SUITE + libc-math-smoke-tests + SRCS + sinpif_test.cpp + DEPENDS + libc.src.errno.errno + libc.src.math.sinpif + libc.src.__support.CPP.array + libc.src.__support.FPUtil.fp_bits +) + add_fp_unittest( sincosf_test SUITE diff --git a/libc/test/src/math/smoke/sinpif_test.cpp b/libc/test/src/math/smoke/sinpif_test.cpp new file mode 100644 index 000000000000..0918294ab361 --- /dev/null +++ b/libc/test/src/math/smoke/sinpif_test.cpp @@ -0,0 +1,43 @@ +//===-- Unittests for sinpif ----------------------------------------------===// +// +// 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/errno/libc_errno.h" +#include "src/math/sinpif.h" +#include "test/UnitTest/FPMatcher.h" + +#include + +using LlvmLibcSinpifTest = LIBC_NAMESPACE::testing::FPTest; + +TEST_F(LlvmLibcSinpifTest, SpecialNumbers) { + LIBC_NAMESPACE::libc_errno = 0; + + EXPECT_FP_EQ(aNaN, LIBC_NAMESPACE::sinpif(aNaN)); + EXPECT_MATH_ERRNO(0); + + EXPECT_FP_EQ(0.0f, LIBC_NAMESPACE::sinpif(0.0f)); + EXPECT_MATH_ERRNO(0); + + EXPECT_FP_EQ(-0.0f, LIBC_NAMESPACE::sinpif(-0.0f)); + EXPECT_MATH_ERRNO(0); + + EXPECT_FP_EQ(aNaN, LIBC_NAMESPACE::sinpif(inf)); + EXPECT_MATH_ERRNO(EDOM); + + EXPECT_FP_EQ(aNaN, LIBC_NAMESPACE::sinpif(neg_inf)); + EXPECT_MATH_ERRNO(EDOM); +} + +TEST_F(LlvmLibcSinpifTest, Integers) { + EXPECT_FP_EQ(-0.0, LIBC_NAMESPACE::sinpif(-0x420)); + EXPECT_FP_EQ(-0.0, LIBC_NAMESPACE::sinpif(-0x1p+43)); + EXPECT_FP_EQ(-0.0, LIBC_NAMESPACE::sinpif(-0x1.4p+64)); + EXPECT_FP_EQ(0.0, LIBC_NAMESPACE::sinpif(0x420)); + EXPECT_FP_EQ(0.0, LIBC_NAMESPACE::sinpif(0x1.cp+106)); + EXPECT_FP_EQ(0.0, LIBC_NAMESPACE::sinpif(0x1.cp+21)); +} diff --git a/libc/utils/MPFRWrapper/MPFRUtils.cpp b/libc/utils/MPFRWrapper/MPFRUtils.cpp index 9a0e498708c0..379a631a356a 100644 --- a/libc/utils/MPFRWrapper/MPFRUtils.cpp +++ b/libc/utils/MPFRWrapper/MPFRUtils.cpp @@ -15,10 +15,7 @@ #include "src/__support/FPUtil/FPBits.h" #include "src/__support/FPUtil/fpbits_str.h" #include "src/__support/macros/properties/types.h" -#include "test/UnitTest/FPMatcher.h" -#include "hdr/math_macros.h" -#include #include #include "mpfr_inc.h" @@ -437,6 +434,23 @@ public: return result; } + MPFRNumber sinpi() const { + MPFRNumber result(*this); + +#if MPFR_VERSION_MAJOR > 4 || \ + (MPFR_VERSION_MAJOR == 4 && MPFR_VERSION_MINOR >= 2) + + mpfr_sinpi(result.value, value, mpfr_rounding); +#else + MPFRNumber value_pi(0.0, 1280); + mpfr_const_pi(value_pi.value, MPFR_RNDN); + mpfr_mul(value_pi.value, value_pi.value, value, MPFR_RNDN); + mpfr_sin(result.value, value_pi.value, mpfr_rounding); +#endif + + return result; + } + MPFRNumber sinh() const { MPFRNumber result(*this); mpfr_sinh(result.value, value, mpfr_rounding); @@ -677,6 +691,8 @@ unary_operation(Operation op, InputType input, unsigned int precision, return mpfrInput.roundeven(); case Operation::Sin: return mpfrInput.sin(); + case Operation::Sinpi: + return mpfrInput.sinpi(); case Operation::Sinh: return mpfrInput.sinh(); case Operation::Sqrt: diff --git a/libc/utils/MPFRWrapper/MPFRUtils.h b/libc/utils/MPFRWrapper/MPFRUtils.h index 46f3375fd4b7..11e323bf6881 100644 --- a/libc/utils/MPFRWrapper/MPFRUtils.h +++ b/libc/utils/MPFRWrapper/MPFRUtils.h @@ -51,6 +51,7 @@ enum class Operation : int { Round, RoundEven, Sin, + Sinpi, Sinh, Sqrt, Tan,