Skip to content

[libc][math][c23] Implemented sinpif function correctly rounded for all rounding modes. - #97149

Merged
lntue merged 19 commits into
llvm:mainfrom
HendrikHuebner:dev/libc-trigpi
Jul 1, 2024
Merged

[libc][math][c23] Implemented sinpif function correctly rounded for all rounding modes.#97149
lntue merged 19 commits into
llvm:mainfrom
HendrikHuebner:dev/libc-trigpi

Conversation

@HendrikHuebner

@HendrikHuebner HendrikHuebner commented Jun 29, 2024

Copy link
Copy Markdown
Contributor

This implements the sinpif function. An exhaustive test shows it's correct for all rounding modes.

Issue: #94895

CC: @lntue

@llvmbot llvmbot added the libc label Jun 29, 2024
@llvmbot

llvmbot commented Jun 29, 2024

Copy link
Copy Markdown
Member

@llvm/pr-subscribers-libc

Author: Hendrik Hübner (HendrikHuebner)

Changes

This implements the sinpif function. An exhaustive test shows it's correct for all rounding modes.

Issue: #94895


Patch is 22.44 KiB, truncated to 20.00 KiB below, full version: https://github.com/llvm/llvm-project/pull/97149.diff

16 Files Affected:

  • (modified) libc/config/linux/x86_64/entrypoints.txt (+1)
  • (modified) libc/docs/math/index.rst (+1-1)
  • (modified) libc/src/math/CMakeLists.txt (+1)
  • (modified) libc/src/math/generic/CMakeLists.txt (+22)
  • (modified) libc/src/math/generic/range_reduction_fma.h (+10)
  • (modified) libc/src/math/generic/sincosf_utils.h (+26-11)
  • (added) libc/src/math/generic/sinpif.cpp (+132)
  • (added) libc/src/math/sinpif.h (+18)
  • (modified) libc/test/src/math/CMakeLists.txt (+16)
  • (modified) libc/test/src/math/exhaustive/CMakeLists.txt (+16)
  • (added) libc/test/src/math/exhaustive/sinpif_test.cpp (+35)
  • (added) libc/test/src/math/sinpif_test.cpp (+115)
  • (modified) libc/test/src/math/smoke/CMakeLists.txt (+13)
  • (added) libc/test/src/math/smoke/sinpif_test.cpp (+38)
  • (modified) libc/utils/MPFRWrapper/MPFRUtils.cpp (+8)
  • (modified) libc/utils/MPFRWrapper/MPFRUtils.h (+1)
diff --git a/libc/config/linux/x86_64/entrypoints.txt b/libc/config/linux/x86_64/entrypoints.txt
index 9581f7f2604c4..9a7ad32b66470 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 56cc8d658257d..9adc93b5e104d 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 3dfc4ac94827d..3037c39c85506 100644
--- a/libc/src/math/CMakeLists.txt
+++ b/libc/src/math/CMakeLists.txt
@@ -370,6 +370,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 3773a2b49c416..de9701dec1594 100644
--- a/libc/src/math/generic/CMakeLists.txt
+++ b/libc/src/math/generic/CMakeLists.txt
@@ -283,6 +283,28 @@ add_entrypoint_object(
     -O3
 )
 
+add_entrypoint_object(
+  sinpif
+  SRCS
+    sinpif.cpp
+  HDRS
+    ../sinpif.h
+  DEPENDS
+    .range_reduction
+    .sincosf_utils
+    libc.src.errno.errno
+    libc.src.__support.FPUtil.basic_operations
+    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.FPUtil.rounding_mode
+    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 82b4ae1c705e1..a270fdc0033f2 100644
--- a/libc/src/math/generic/range_reduction_fma.h
+++ b/libc/src/math/generic/range_reduction_fma.h
@@ -38,6 +38,16 @@ LIBC_INLINE int64_t small_range_reduction(double x, double &y) {
   return static_cast<int64_t>(kd);
 }
 
+// Return k and y, where
+//   k = round(x * 32) and y = (x * 32) - k.
+//   => pi * x = (k + y) * pi / 32
+LIBC_INLINE int64_t small_range_reduction_mul_pi(double x, double &y) {
+  double kd = fputil::nearest_integer(x * 32);
+  y = fputil::fma<double>(x, 32.0, -kd);
+  
+  return static_cast<int64_t>(kd);
+}
+
 // Return k and y, where
 //   k = round(x * 32 / pi) and y = (x * 32 / pi) - k.
 // This is used for sinf, cosf, sincosf.
diff --git a/libc/src/math/generic/sincosf_utils.h b/libc/src/math/generic/sincosf_utils.h
index 904df4f4ed08e..5d93b287057a1 100644
--- a/libc/src/math/generic/sincosf_utils.h
+++ b/libc/src/math/generic/sincosf_utils.h
@@ -18,7 +18,9 @@
 // using namespace LIBC_NAMESPACE::fma;
 using LIBC_NAMESPACE::fma::FAST_PASS_BOUND;
 using LIBC_NAMESPACE::fma::large_range_reduction;
+using LIBC_NAMESPACE::fma::small_range_reduction_mul_pi;
 using LIBC_NAMESPACE::fma::small_range_reduction;
+
 #else
 #include "range_reduction.h"
 // using namespace LIBC_NAMESPACE::generic;
@@ -58,18 +60,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<float> 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 +88,28 @@ 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<float> 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);
+}
+
+LIBC_INLINE void sincospif_eval(double xd, double &sin_k, double &cos_k,
+                                double &sin_y, double &cosm1_y) {
+  double y;
+  int64_t k = small_range_reduction_mul_pi(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 0000000000000..2cedb6f29021b
--- /dev/null
+++ b/libc/src/math/generic/sinpif.cpp
@@ -0,0 +1,132 @@
+//===-- Single-precision sinpi 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/BasicOperations.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/FPUtil/rounding_mode.h"
+#include "src/__support/common.h"
+#include "src/__support/macros/optimization.h"            // LIBC_UNLIKELY
+#include "src/__support/macros/properties/cpu_features.h" // LIBC_TARGET_CPU_HAS_FMA
+
+#include <cstdio>
+#include <errno.h>
+#include <iostream>
+
+#if defined(LIBC_TARGET_CPU_HAS_FMA)
+#include "range_reduction_fma.h"
+#else
+#include "range_reduction.h"
+#endif
+
+namespace LIBC_NAMESPACE {
+
+LLVM_LIBC_FUNCTION(float, sinpif, (float x)) {
+  using FPBits = typename fputil::FPBits<float>;
+  FPBits xbits(x);
+
+  uint32_t x_u = xbits.uintval();
+  uint32_t x_abs = x_u & 0x7fff'ffffU;
+  double xd = static_cast<double>(x);
+
+  // Range reduction:
+  // For |x| > pi/32, we perform range reduction as follows:
+  // Find k and y such that:
+  //   x = (k + y) * pi/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/pi)
+  //   y = x * 32/pi - k
+  // For large range, we will omit all the higher parts of 32/pi such that the
+  // least significant bits of their full products with x are larger than 63,
+  // since sin((k + y + 64*i) * pi/32) = sin(x + i * 2pi) = sin(x).
+  //
+  // When FMA instructions are not available, we store the digits of 32/pi in
+  // chunks of 28-bit precision.  This will make sure that the products:
+  //   x * THIRTYTWO_OVER_PI_28[i] are all exact.
+  // When FMA instructions are available, we simply store the digits of 32/pi in
+  // chunks of doubles (53-bit of precision).
+  // So when multiplying by the largest values of single precision, the
+  // resulting output should be correct up to 2^(-208 + 128) ~ 2^-80.  By the
+  // worst-case analysis of range reduction, |y| >= 2^-38, so this should give
+  // us more than 40 bits of accuracy. For the worst-case estimation of range
+  // reduction, see for instances:
+  //   Elementary Functions by J-M. Muller, Chapter 11,
+  //   Handbook of Floating-Point Arithmetic by J-M. Muller et. al.,
+  //   Chapter 10.2.
+  //
+  // Once k and y are computed, we then deduce the answer by the sine of sum
+  // formula:
+  //   sin(x) = 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<float>(xdpi);
+    }
+
+    // |x| < 1/16.
+    double xsq = xd * xd;
+
+    // Degree-9 polynomial approximation:
+    //   sin(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<float>(xd * result);
+  }
+
+
+  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();
+  }
+
+  // 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);
+
+  return static_cast<float>(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 0000000000000..a5bef435efae6
--- /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 36d2a2fbfd302..6e9ea9bd0466b 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 fb3596c3378ff..412ca031d0e99 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 0000000000000..3c74c2872bef0
--- /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 "src/math/sinpif.h"
+#include "utils/MPFRWrapper/MPFRUtils.h"
+#include "mpfr.h"
+#include <sys/types.h>
+
+namespace mpfr = LIBC_NAMESPACE::testing::mpfr;
+
+using LlvmLibcSinpifExhaustiveTest =
+    LlvmLibcUnaryOpExhaustiveMathTest<float, mpfr::Operation::Sinpi,
+                                      LIBC_NAMESPACE::sinpif>;
+
+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 0000000000000..15f25b86574ac
--- /dev/null
+++ b/libc/test/src/math/sinpif_test.cpp
@@ -0,0 +1,115 @@
+//===-- 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 "hdr/math_macros.h"
+#include "src/__support/FPUtil/FPBits.h"
+#include "src/errno/libc_errno.h"
+#include "src/math/sinpif.h"
+#include "test/UnitTest/FPMatcher.h"
+#include "test/UnitTest/Test.h"
+#include "test/src/math/sdcomp26094.h"
+#include "utils/MPFRWrapper/MPFRUtils.h"
+
+#include <errno.h>
+#include <stdint.h>
+
+using LlvmLibcSinpifTest = LIBC_NAMESPACE::testing::FPTest<float>;
+
+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);
+  }
+}
diff --git a/libc/test/src/math/smoke/CMakeLists.txt b/libc/test/src/math/smoke/CMakeLists.txt
index ee6f1595fe0fd..c58167a150adc 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 m...
[truncated]

Comment thread libc/src/math/generic/range_reduction_fma.h Outdated
@github-actions

github-actions Bot commented Jun 29, 2024

Copy link
Copy Markdown

✅ With the latest revision this PR passed the C/C++ code formatter.

@lntue
lntue self-requested a review June 29, 2024 12:05
@lntue lntue changed the title [libc][math] Implemented sinpif function [libc][math][c23] Implemented sinpif function Jun 29, 2024
@lntue lntue changed the title [libc][math][c23] Implemented sinpif function [libc][math][c23] Implemented sinpif function correctly rounded for all rounding modes. Jun 29, 2024
Comment thread libc/config/linux/x86_64/entrypoints.txt
Comment thread libc/src/math/generic/range_reduction_fma.h Outdated
Comment thread libc/src/math/generic/sinpif.cpp Outdated
Comment thread libc/src/math/generic/sinpif.cpp Outdated
Comment thread libc/test/src/math/sinpif_test.cpp Outdated
Comment thread libc/test/src/math/smoke/sinpif_test.cpp Outdated
@lntue

lntue commented Jun 29, 2024

Copy link
Copy Markdown
Contributor

⚠️ C/C++ code formatter, clang-format found issues in your code. ⚠️

You can test this locally with the following command:
View the diff from clang-format here.

Please fix formatting.

@HendrikHuebner

Copy link
Copy Markdown
Contributor Author

Alright I made the changes you suggested and removed all unnecessary includes. I noticed a bunch of unused includes in some other tests but that should probably be a separate PR

Comment thread libc/test/src/math/smoke/sinpif_test.cpp Outdated
Comment thread libc/test/src/math/sinpif_test.cpp Outdated
Comment thread libc/test/src/math/exhaustive/sinpif_test.cpp Outdated
Comment thread libc/src/math/generic/sinpif.cpp Outdated
Comment thread libc/utils/MPFRWrapper/MPFRUtils.cpp Outdated
@lntue

lntue commented Jun 29, 2024

Copy link
Copy Markdown
Contributor

@zimmermann6

@zimmermann6

Copy link
Copy Markdown

I get some failures with the CORE-MATH test suite for integer n. IEEE 754-2019 says that sin(+n) = +0 and sin(-n) = -0.

Running exhaustive check in --rndn mode...
FAIL x=-0x1.cp+21 ref=-0x0p+0 y=0x0p+0
FAIL x=-0x1.2p+96 ref=-0x0p+0 y=0x0p+0
FAIL x=-0x1.dp+5 ref=-0x0p+0 y=0x0p+0
Running exhaustive check in --rndz mode...
FAIL x=-0x1.1p+112 ref=-0x0p+0 y=0x0p+0
FAIL x=-0x1.ap+53 ref=-0x0p+0 y=0x0p+0
Running exhaustive check in --rndu mode...
FAIL x=-0x1.fp+58 ref=-0x0p+0 y=0x0p+0
FRunning exhaustive check in --rndd mode...
FAIL x=0x1.7p+101 ref=0x0p+0 y=-0x0p+0

@HendrikHuebner

Copy link
Copy Markdown
Contributor Author

@zimmermann6 How do I reproduce this?

Comment thread libc/utils/MPFRWrapper/MPFRUtils.cpp Outdated
Comment thread libc/utils/MPFRWrapper/MPFRUtils.cpp Outdated
@zimmermann6

Copy link
Copy Markdown

Hendrik you can reproduce with the following command within CORE-MATH:

CORE_MATH_CHECK_STD=true LIBM=$L ./check.sh sinpif

where L points to libllvmlibc.a (you need MPFR 4.2.0 or later)

Comment thread libc/src/math/generic/sinpif.cpp Outdated
sincospif_eval(xd, sin_k, cos_k, sin_y, cosm1_y);

if (LIBC_UNLIKELY(sin_y == 0 && sin_k == 0))
return fputil::copysign(0.0f, x);

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This handles signed/unsigned zeros

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Another option if you don't want to bring in src/__support/FPUtils/ManipulationFunctions.h dependency is to use:

  return FPBits::zero(xbits.sign());

Comment thread libc/src/math/generic/sinpif.cpp Outdated
return static_cast<float>(xd * result);
}

if (LIBC_UNLIKELY(x_abs >= 0x7f80'0000U)) {

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think we can speed thing up a bit by doing quick return for |x| >= 2^23, since it will then always be an integer.

@zimmermann6

Copy link
Copy Markdown

all tests do pass now, thanks!

@lntue lntue left a comment

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Everything looks good to me. Thanks!

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

Projects

None yet

Development

Successfully merging this pull request may close these issues.

4 participants