[libc][math][c23] add sinpi function

libc/config/darwin/arm/entrypoints.txt

@@ -255,6 +255,7 @@ set(TARGET_LIBM_ENTRYPOINTS
     libc.src.math.sin
     libc.src.math.sinf
     libc.src.math.sinpif
+    libc.src.math.sinpi
     libc.src.math.sqrt
     libc.src.math.sqrtf
     libc.src.math.sqrtl

libc/config/linux/aarch64/entrypoints.txt

@@ -610,6 +610,7 @@ set(TARGET_LIBM_ENTRYPOINTS
     libc.src.math.sinf
     libc.src.math.sinhf
     libc.src.math.sinpif
+    libc.src.math.sinpi
     libc.src.math.sqrt
     libc.src.math.sqrtf
     libc.src.math.sqrtl

libc/include/math.yaml

@@ -2541,6 +2541,12 @@ functions:
     arguments:
       - type: const long double *
       - type: const long double *
+  - name: sinpi
+    standards:
+      - stdc
+    return_type: double
+    arguments:
+      - type: double
   - name: totalordermag
     standards:
       - stdc

libc/src/math/CMakeLists.txt

@@ -495,6 +495,8 @@ add_math_entrypoint_object(sincosf)
 add_math_entrypoint_object(sin)
 add_math_entrypoint_object(sinf)
 add_math_entrypoint_object(sinf16)
+
+add_math_entrypoint_object(sinpi)
 add_math_entrypoint_object(sinpif)
 add_math_entrypoint_object(sinpif16)
 

libc/src/math/generic/CMakeLists.txt

@@ -1,5 +1,3 @@
-
-
 add_entrypoint_object(
   canonicalize
   SRCS
@@ -495,6 +493,29 @@ add_entrypoint_object(
     libc.src.__support.macros.optimization
 )
 
+add_entrypoint_object(
+  sinpi
+  SRCS
+    sinpi.cpp
+  HDRS
+    ../sinpi.h
+  DEPENDS
+    .sincosf_utils
+    libc.src.__support.FPUtil.fenv_impl
+    libc.src.__support.FPUtil.fp_bits
+    libc.src.__support.FPUtil.nearest_integer
+    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
+    libc.src.__support.FPUtil.basic_operations
+    libc.src.math.generic.sincos_eval
+    libc.src.__support.FPUtil.double_double
+    libc.src.__support.FPUtil.multiply_add
+    .range_reduction_double
+)
+
 add_entrypoint_object(
   sinpif
   SRCS

libc/src/math/generic/sinpi.cpp

@@ -0,0 +1,95 @@
+#include "src/math/sinpi.h"
+#include "sincos_eval.h"
+#include "src/__support/FPUtil/BasicOperations.h"
+#include "src/__support/FPUtil/FEnvImpl.h"
+#include "src/__support/FPUtil/FPBits.h"
+#include "src/__support/FPUtil/multiply_add.h"
+#include "src/__support/common.h"
+#include "src/__support/macros/config.h"
+#include "src/__support/FPUtil/double_double.h"
+#include "src/__support/FPUtil/generic/mul.h"
+#include "src/__support/FPUtil/nearest_integer.h"
+#include "src/math/pow.h"
+//#include "range_reduction_double_nofma.h"
+//#include "src/__support/FPUtil/multiply_add.h"
+//#include "src/math/generic/range_reduction_double_common.h"
+#include <iostream>
+
+namespace LIBC_NAMESPACE_DECL {
+
+LLVM_LIBC_FUNCTION(double, sinpi, (double x)) {
+  // Given x * pi = y - (k * (pi/128)) 
+  // find y and k such that
+  // y = x * pi - (k * pi/128) = x * pi - kpi/128
+  // k = round(x, 128)
+
+  using FPBits = typename fputil::FPBits<double>;
+  using DoubleDouble = fputil::DoubleDouble;
+
+  FPBits xbits(x);
+
+  double k = fputil::nearest_integer(x * 128);
+  int  k_int = static_cast<int>(k);
+
+  std::cout << "k" << k << std::endl;
+
+  double yk = x - k/128;
+
+  DoubleDouble yy = fputil::exact_mult(yk, 3.141592653589793115997963468544185161590576171875);
+  yy.lo = fputil::multiply_add(yk, 1.2246467991473532071737640294583966046256921246776e-16, yy.lo);
+
+  uint64_t abs_u = xbits.uintval();
+
+  uint64_t x_abs = abs_u & 0xFFFFFFFFFFFFFFFF;
+
+  if (LIBC_UNLIKELY(x_abs == 0U))
+    return x;
+  if (x_abs >= 0x4330000000000000) {
+    if (xbits.is_nan())
+      return x;
+    if (xbits.is_inf()) {
+      fputil::set_errno_if_required(EDOM);
+      fputil::raise_except_if_required(FE_INVALID);
+      return FPBits::quiet_nan().get_val();
+    }
+    return FPBits::zero(xbits.sign()).get_val();
+  }
+
+  DoubleDouble sin_y, cos_y;
+
+  [[maybe_unused]] double err = generic::sincos_eval(yy, sin_y, cos_y);
+  DoubleDouble sin_k = SIN_K_PI_OVER_128[k_int & 255];
+  DoubleDouble cos_k = SIN_K_PI_OVER_128[(k_int + 64) & 255];
+
+  std::cout << "sin_k: " << sin_k.hi << std::endl;
+  std::cout << "sin_klo: " << sin_k.lo << std::endl;
+  std::cout << "sin_y: " << sin_y.hi << std::endl;
+  std::cout << "cos_y: " << cos_y.hi << std::endl;
+  std::cout << "sin_y.lo: " << sin_y.lo << std::endl;
+  std::cout << "cos_y.o: " << cos_y.lo << std::endl;
+
+  double cosm1_y = cos_y.hi - 1.0;
+  DoubleDouble sin_y_cos_k = fputil::quick_mult(sin_y, cos_k);
+
+  std::cout << "cosm1" << cosm1_y << std::endl;
+  DoubleDouble cosm1_yy;
+  cosm1_yy.hi = cosm1_y;
+  cosm1_yy.lo = 0.0;
+
+  DoubleDouble cos_y_sin_k = fputil::quick_mult(cos_y, sin_k);
+  DoubleDouble rr = fputil::exact_add<false>(sin_y_cos_k.hi, cos_y_sin_k.hi);
+
+  std::cout << "r.hi:" << rr.hi << std::endl;
+  std::cout << "r.lo" << rr.lo << std::endl;
+
+  rr.lo += sin_y_cos_k.lo + cos_y_sin_k.lo;
+
+  std::cout << "rrlo2: " << rr.lo << std::endl;  
+  std::cout << "cos_y_sin_k:" << cos_y_sin_k.hi << std::endl;
+  std::cout << "siny*cosk.lo:" << sin_y_cos_k.lo << std::endl;
+  std::cout << "rrhi + rrlo + sink.hi " << rr.hi + rr.lo + sin_k.hi + sin_k.lo << std::endl;
+  std::cout << "rrhi + rrlo " << rr.hi + rr.lo << std::endl;
+
+  return rr.hi + rr.lo;
+ }
+} // namespace LIBC_NAMESPACE_DECL

libc/src/math/sinpi.h

@@ -0,0 +1,20 @@
+//===-- Implementation header for sinpi -------------------------*- C++ -*-===//
+//
+// Part of the LLVM Project, under the Apache Licese 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_SINPI_H
+#define LLVM_LIBC_SRC_MATH_SINPI_H
+
+#include "src/__support/macros/config.h"
+
+namespace LIBC_NAMESPACE_DECL {
+
+double sinpi(double x);
+
+} // namespace LIBC_NAMESPACE_DECL
+
+#endif // LLVM_LIBC_SRC_MATH_SINPI_H

libc/test/src/math/CMakeLists.txt

@@ -106,6 +106,21 @@ add_fp_unittest(
   DEPENDS
     libc.src.math.sinf16
 )
+add_fp_unittest(
+  sinpi_test
+  NEED_MPFR
+  SUITE
+    libc-math-unittests
+  SRCS
+    sinpi_test.cpp
+  HDRS
+    sdcomp26094.h
+  DEPENDS
+    libc.src.errno.errno
+    libc.src.math.sinpi
+    libc.src.__support.CPP.array
+    libc.src.__support.FPUtil.fp_bits
+)
 
 add_fp_unittest(
   sinpif_test

libc/test/src/math/sinpi_test.cpp

@@ -0,0 +1,87 @@
+//===-- Exhaustive test for sinpif16 --------------------------------------===//
+//
+// 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/sinpi.h"
+#include "test/UnitTest/FPMatcher.h"
+#include "test/UnitTest/Test.h"
+#include "utils/MPFRWrapper/MPFRUtils.h"
+
+#include <iostream>
+
+using LlvmLibcSinpiTest = LIBC_NAMESPACE::testing::FPTest<double>;
+using namespace std;
+namespace mpfr = LIBC_NAMESPACE::testing::mpfr;
+
+using LIBC_NAMESPACE::testing::tlog;
+
+TEST_F(LlvmLibcSinpiTest, InDoubleRange) {
+  constexpr uint64_t COUNT = 1'234'51;
+  uint64_t START = LIBC_NAMESPACE::fputil::FPBits<double>(0x1.0p-50).uintval();
+  uint64_t STOP = LIBC_NAMESPACE::fputil::FPBits<double>(0x1.0p200).uintval();
+  uint64_t STEP = (STOP - START) / COUNT;
+
+  auto test = [&](mpfr::RoundingMode rounding_mode) {
+    mpfr::ForceRoundingMode __r(rounding_mode);
+    if (!__r.success)
+      return;
+
+    uint64_t fails = 0;
+    uint64_t count = 0;
+    uint64_t cc = 0;
+    double mx, mr = 0.0;
+    double tol = 2.0;
+
+    for (uint64_t i = 0, v = START; i <= COUNT; ++i, v += STEP) {
+      double x = FPBits(v).get_val();
+      if (FPBits(v).is_nan() || FPBits(v).is_inf())
+        continue;
+      LIBC_NAMESPACE::libc_errno = 0;
+      double result = LIBC_NAMESPACE::sinpi(x);
+      std::cout << "result: " << result << std::endl;
+      ++cc;
+      if (FPBits(result).is_nan() || FPBits(result).is_inf())
+        continue;
+
+      ++count;
+
+
+      if (!TEST_MPFR_MATCH_ROUNDING_SILENTLY(mpfr::Operation::Sinpi, x, result,
+                                             2.0, rounding_mode)) {
+        ++fails;
+        while (!TEST_MPFR_MATCH_ROUNDING_SILENTLY(mpfr::Operation::Sinpi, x,
+                                                  result, tol, rounding_mode)) {
+          mx = x;
+          mr = result;
+
+          if (tol > 1000.0)
+            break;
+
+          tol *= 2.0;
+        }
+      }
+    }
+    if (fails) {
+      tlog << " Sin failed: " << fails << "/" << count << "/" << cc
+           << " tests.\n";
+      tlog << "   Max ULPs is at most: " << static_cast<uint64_t>(tol) << ".\n";
+      EXPECT_MPFR_MATCH(mpfr::Operation::Sinpi, mx, mr, 2.0, rounding_mode);
+    }
+
+  };
+  tlog << " Test Rounding To Nearest...\n";
+  test(mpfr::RoundingMode::Nearest);
+
+  tlog << " Test Rounding Downward...\n";
+  test(mpfr::RoundingMode::Downward);
+
+  tlog << " Test Rounding Upward...\n";
+  test(mpfr::RoundingMode::Upward);
+
+  tlog << " Test Rounding Toward Zero...\n";
+  test(mpfr::RoundingMode::TowardZero);
+}

libc/test/src/math/smoke/CMakeLists.txt

@@ -73,6 +73,18 @@ add_fp_unittest(
     libc.src.__support.FPUtil.cast
 )
 
+add_fp_unittest(
+  sinpi_test
+  SUITE
+    libc-math-smoke-tests
+  SRCS
+    sinpi_test.cpp
+  DEPENDS
+    libc.src.errno.errno
+    libc.src.math.sinpi
+    libc.src.__support.CPP.array
+    libc.src.__support.FPUtil.fp_bits
+)
 add_fp_unittest(
   sinpif_test
   SUITE

libc/test/src/math/smoke/sinpi_test.cpp

@@ -0,0 +1,59 @@
+//===-- Unittests for sinpi -----------------------------------------------===//
+//
+// 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/cast.h"
+#include "src/errno/libc_errno.h"
+#include "src/math/sinpi.h"
+#include "test/UnitTest/FPMatcher.h"
+#include "test/UnitTest/Test.h"
+
+using LlvmLibcSinpiTest = LIBC_NAMESPACE::testing::FPTest<double>;
+/*
+TEST_F(LlvmLibcSinpiTest, SpecialNumbers) {
+  LIBC_NAMESPACE::libc_errno = 0;
+
+  EXPECT_FP_EQ(aNaN, LIBC_NAMESPACE::sinpi(aNaN));
+  EXPECT_MATH_ERRNO(0);
+
+  EXPECT_FP_EQ(zero, LIBC_NAMESPACE::sinpi(zero));
+  EXPECT_MATH_ERRNO(0);
+
+  EXPECT_FP_EQ(neg_zero, LIBC_NAMESPACE::sinpi(neg_zero));
+  EXPECT_MATH_ERRNO(0);
+
+  EXPECT_FP_EQ(aNaN, LIBC_NAMESPACE::sinpi(inf));
+  EXPECT_MATH_ERRNO(EDOM);
+
+  EXPECT_FP_EQ(aNaN, LIBC_NAMESPACE::sinpi(neg_inf));
+  EXPECT_MATH_ERRNO(EDOM);
+}
+
+TEST_F(LlvmLibcSinpiTest, Integers) {
+  EXPECT_FP_EQ(neg_zero, LIBC_NAMESPACE::sinpi(-0x1.0000000000003p52));
+
+  ASSERT_FP_EQ(-1.0, LIBC_NAMESPACE::sinpi(4499003037990983.5));
+
+  EXPECT_FP_EQ(neg_zero, LIBC_NAMESPACE::sinpi(-0x1.0000000000005p52));
+
+  EXPECT_FP_EQ(neg_zero, LIBC_NAMESPACE::sinpi(-0x1.0000000000006p52));
+
+  EXPECT_FP_EQ(zero, LIBC_NAMESPACE::sinpi(0x1.0p1020));
+
+  EXPECT_FP_EQ(zero, LIBC_NAMESPACE::sinpi(0x1.0000000000003p52));
+
+  EXPECT_FP_EQ(zero, LIBC_NAMESPACE::sinpi(0x1.0000000000005p52));
+
+  EXPECT_FP_EQ(zero, LIBC_NAMESPACE::sinpi(0x1.0000000000006p52));
+}
+*/
+
+TEST_F(LlvmLibcSinpiTest, Integers) {
+  //ASSERT_FP_EQ(zero, LIBC_NAMESPACE::sinpi(4503563146482784.00000000000000000000000000000000000000000000000000));
+
+  ASSERT_FP_EQ(-1.0, LIBC_NAMESPACE::sinpi(4499003037990983.50000000000000000000000000000000000000000000000000));
+}