[libc][math][C23] implementation of cospi

libc/config/darwin/arm/entrypoints.txt

@@ -139,6 +139,7 @@ set(TARGET_LIBM_ENTRYPOINTS
     libc.src.math.coshf
     libc.src.math.cos
     libc.src.math.cosf
+    libc.src.math.cospi
     libc.src.math.cospif
     libc.src.math.dfmal
     libc.src.math.dsqrtl

libc/config/linux/aarch64/entrypoints.txt

@@ -424,6 +424,7 @@ set(TARGET_LIBM_ENTRYPOINTS
     libc.src.math.cos
     libc.src.math.cosf
     libc.src.math.coshf
+    libc.src.math.cospi
     libc.src.math.cospif
     libc.src.math.daddl
     libc.src.math.ddivl

libc/include/math.yaml

@@ -249,6 +249,12 @@ functions:
     arguments:
       - type: long double
       - type: long double
+  - name: cospi
+    standards:
+      - stdc
+    return_type: double
+    arguments:
+      - type: double
   - name: dmulf128
     standards:
       - llvm_libc_ext

libc/src/math/CMakeLists.txt

@@ -99,6 +99,7 @@ add_math_entrypoint_object(cosh)
 add_math_entrypoint_object(coshf)
 add_math_entrypoint_object(coshf16)
 
+add_math_entrypoint_object(cospi)
 add_math_entrypoint_object(cospif)
 add_math_entrypoint_object(cospif16)
 

libc/src/math/cospi.h

@@ -0,0 +1,20 @@
+//===-- Implementation header for cospi -------------------------*- 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_COSPI_H
+#define LLVM_LIBC_SRC_MATH_COSPI_H
+
+#include "src/__support/macros/config.h"
+
+namespace LIBC_NAMESPACE_DECL {
+
+double cospi(double x);
+
+} // namespace LIBC_NAMESPACE_DECL
+
+#endif // LLVM_LIBC_SRC_MATH_SINPIF_H

libc/src/math/generic/CMakeLists.txt

@@ -1,5 +1,3 @@
-
-
 add_entrypoint_object(
   canonicalize
   SRCS
@@ -381,6 +379,23 @@ add_entrypoint_object(
     libc.src.__support.macros.optimization
     libc.src.__support.macros.properties.types
 )
+add_entrypoint_object(
+  cospi
+  SRCS
+    cospi.cpp
+  HDRS
+    ../cospi.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.macros.optimization
+    libc.src.__support.FPUtil.polyeval
+    libc.src.math.pow
+    libc.src.__support.FPUtil.basic_operations
+)
 
 add_entrypoint_object(
   cospif

libc/src/math/generic/cospi.cpp

@@ -0,0 +1,97 @@
+//===-- double-precision cospi 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/cospi.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/PolyEval.h"
+#include "src/__support/FPUtil/double_double.h"
+#include "src/__support/FPUtil/multiply_add.h"
+#include "src/__support/common.h"
+#include "src/__support/macros/config.h"
+#include "src/__support/macros/optimization.h" // LIBC_UNLIKELY
+#include "src/math/fmul.h"
+#include "src/math/generic/sincosf_utils.h"
+#include "src/math/pow.h"
+
+namespace LIBC_NAMESPACE_DECL {
+
+static LIBC_INLINE void sincospi_poly_eval(double k, double y, double &sin_k,
+                                           double &cos_k, double &sin_y,
+
+                                           double &cosm1_y) {
+
+  // Q3 = fpminimax(sin(x*pi), 7, [|64...|], [-0.0078125, 0.0078125]);
+  sin_k =
+      k * fputil::polyeval(k, 0x1.59b6a771a45cbab8p-94, 0x1.921fb54442d1846ap1,
+                           -0x1.8633470ba8bd806cp-76, -0x1.4abbce625be56346p2,
+                           0x1.d3e01dfd72e97a92p-61, 0x1.466bc67713dbbfp1,
+                           -0x1.14c2648595e2ad4p-47, -0x1.32d1cc20b89301fcp-1);
+
+  sin_y =
+      y * fputil::polyeval(k, 0x1.59b6a771a45cbab8p-94, 0x1.921fb54442d1846ap1,
+                           -0x1.8633470ba8bd806cp-76, -0x1.4abbce625be56346p2,
+                           0x1.d3e01dfd72e97a92p-61, 0x1.466bc67713dbbfp1,
+                           -0x1.14c2648595e2ad4p-47, -0x1.32d1cc20b89301fcp-1);
+
+  // Q1 = fpminimax(cos(x * pi), 7, [|64...|], [-0.0078125, 0.0078125]);
+  cos_k =
+      k * fputil::polyeval(k, 0x1p0, 0x1.a5b22c564ee1d862p-84,
+                           -0x1.3bd3cc9be45d30e6p2, -0x1.5c2328fefbe60d3ep-66,
+                           0x1.03c1f080a6907a6p2, 0x1.569a4d5c5018eecap-51,
+                           -0x1.55d1f72455a9848ap0, -0x1.6b18e5f7fc6c39a6p-38);
+
+  cosm1_y =
+      y * fputil::polyeval(y, 0x1p0, 0x1.a5b22c564ee1d862p-84,
+                           -0x1.3bd3cc9be45d30e6p2, -0x1.5c2328fefbe60d3ep-66,
+                           0x1.03c1f080a6907a6p2, 0x1.569a4d5c5018eecap-51,
+                           -0x1.55d1f72455a9848ap0, -0x1.6b18e5f7fc6c39a6p-38);
+}
+
+LLVM_LIBC_FUNCTION(double, cospi, (double x)) {
+  using FPBits = typename fputil::FPBits<double>;
+  FPBits xbits(x);
+
+  xbits.set_sign(Sign::POS);
+
+  uint64_t x_abs_ = xbits.uintval();
+  double x_abs = fputil::abs(x);
+  double p = 0x1p52; // 2^p where p is the precision
+  double p2 = 0x1p53;
+  double p3 = 1.0;
+  if (LIBC_UNLIKELY(x_abs == 0U))
+    return p3;
+
+  if (x_abs >= p) {
+    if (x_abs < p2)
+      return ((x_abs_ & 0x1) ? -p3 : p3);
+    if (xbits.is_nan())
+      return x;
+    if (xbits.is_inf()) {
+      fputil::set_errno_if_required(EDOM);
+      fputil::raise_except_if_required(FE_INVALID);
+      return x + FPBits::quiet_nan().get_val();
+    }
+    return p3;
+  }
+  double n = pow(2, -52);
+  double k = fputil::nearest_integer(x * n);
+  double y = x - k;
+  double sin_k, cos_k, sin_y, cosm1_y;
+
+  sincospi_poly_eval(k, y, sin_k, cos_k, sin_y, cosm1_y);
+
+  if (LIBC_UNLIKELY(sin_y == 0 && cos_k == 0))
+    return FPBits::zero(xbits.sign()).get_val();
+
+  return fputil::cast<double>(fputil::multiply_add(
+      cos_k, cosm1_y, fputil::multiply_add(-sin_k, sin_y, cos_k)));
+}
+} // namespace LIBC_NAMESPACE_DECL

libc/test/src/math/CMakeLists.txt

@@ -39,6 +39,22 @@ add_fp_unittest(
     libc.src.math.cosf16
 )
 
+add_fp_unittest(
+  cospi_test
+  NEED_MPFR
+  SUITE
+    libc-math-unittests
+  SRCS
+    cospi_test.cpp
+  HDRS
+    sdcomp26094.h
+  DEPENDS
+    libc.src.errno.errno
+    libc.src.math.cospi
+    libc.src.__support.CPP.array
+    libc.src.__support.FPUtil.fp_bits
+)
+
 add_fp_unittest(
   cospif_test
   NEED_MPFR

libc/test/src/math/cospi_test.cpp

@@ -0,0 +1,30 @@
+//===-- Exhaustive test for cospi -----------------------------------------===//
+//
+// 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/cospi.h"
+#include "test/UnitTest/FPMatcher.h"
+#include "test/UnitTest/Test.h"
+#include "utils/MPFRWrapper/MPFRUtils.h"
+
+#include <iostream>
+
+using LlvmLibcCospiTest = LIBC_NAMESPACE::testing::FPTest<double>;
+using namespace std;
+namespace mpfr = LIBC_NAMESPACE::testing::mpfr;
+
+static constexpr double POS_START = 0;
+static constexpr double POS_STOP = 200;
+
+TEST_F(LlvmLibcCospiTest, PositiveRange) {
+  for (double v = POS_START; v <= POS_STOP; ++v) {
+    EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Cospi, v,
+                                   LIBC_NAMESPACE::cospi(v), 0.5);
+    EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Cospi, -v,
+                                   LIBC_NAMESPACE::cospi(-v), 0.5);
+  }
+}

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

@@ -23,6 +23,18 @@ add_fp_unittest(
     libc.src.math.cosf16
 )
 
+add_fp_unittest(
+  cospi_test
+  SUITE
+    libc-math-smoke-tests
+  SRCS
+    cospi_test.cpp
+  DEPENDS
+    libc.src.errno.errno
+    libc.src.math.cospi
+    libc.src.__support.CPP.array
+    libc.src.__support.FPUtil.fp_bits
+)
 add_fp_unittest(
   cospif_test
   SUITE

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

@@ -0,0 +1,45 @@
+//===-- Unittests for cospi -----------------------------------------------===//
+//
+// 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/cospi.h"
+#include "test/UnitTest/FPMatcher.h"
+#include "test/UnitTest/Test.h"
+
+using LlvmLibcCospiTest = LIBC_NAMESPACE::testing::FPTest<double>;
+
+TEST_F(LlvmLibcCospiTest, SpecialNumbers) {
+  LIBC_NAMESPACE::libc_errno = 0;
+
+  EXPECT_FP_EQ(aNaN, LIBC_NAMESPACE::cospi(aNaN));
+  EXPECT_MATH_ERRNO(0);
+
+  EXPECT_FP_EQ(FPBits::one().get_val(), LIBC_NAMESPACE::cospi(zero));
+  EXPECT_MATH_ERRNO(0);
+
+  EXPECT_FP_EQ(FPBits::one().get_val(), LIBC_NAMESPACE::cospi(neg_zero));
+  EXPECT_MATH_ERRNO(0);
+
+  EXPECT_FP_EQ(aNaN, LIBC_NAMESPACE::cospi(inf));
+  EXPECT_MATH_ERRNO(EDOM);
+
+  EXPECT_FP_EQ(aNaN, LIBC_NAMESPACE::cospi(neg_inf));
+  EXPECT_MATH_ERRNO(EDOM);
+}
+
+TEST_F(LlvmLibcCospiTest, Integers) {
+
+  EXPECT_FP_EQ(FPBits::one(Sign::NEG).get_val(), LIBC_NAMESPACE::cospi(-0x1.8406003b2ae63p52));
+
+  EXPECT_FP_EQ(FPBits::one().get_val(), LIBC_NAMESPACE::cospi(0x1p54));
+
+  EXPECT_FP_EQ(FPBits::one().get_val(), LIBC_NAMESPACE::cospi(0x1p55));
+
+  EXPECT_FP_EQ(FPBits::one().get_val(), LIBC_NAMESPACE::cospi(0x1p56));
+}