[libc-commits] [libc] [libc][math][C23] Implement double precision sinpi correctly rounded for all rounding modes (PR #134921)
via libc-commits
libc-commits at lists.llvm.org
Tue Apr 8 12:48:17 PDT 2025
https://github.com/laraabcd created https://github.com/llvm/llvm-project/pull/134921
This is the implementation of `sinpi` using double-double arithmetic.
>From fdfdd7ed5936de6567d7dd565e831bed3e9bad7c Mon Sep 17 00:00:00 2001
From: Job Hernandez Lara <lararicciab at proton.me>
Date: Tue, 8 Apr 2025 12:45:44 -0700
Subject: [PATCH] add sinpi and tests
---
libc/src/math/generic/sinpi.cpp | 95 +++++++++++++++++++++++++
libc/test/src/math/sinpi_tests.cpp | 88 +++++++++++++++++++++++
libc/test/src/math/smoke/sinpi_test.cpp | 47 ++++++++++++
3 files changed, 230 insertions(+)
create mode 100644 libc/src/math/generic/sinpi.cpp
create mode 100644 libc/test/src/math/sinpi_tests.cpp
create mode 100644 libc/test/src/math/smoke/sinpi_test.cpp
diff --git a/libc/src/math/generic/sinpi.cpp b/libc/src/math/generic/sinpi.cpp
new file mode 100644
index 0000000000000..31aa5f105617c
--- /dev/null
+++ b/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_DE
diff --git a/libc/test/src/math/sinpi_tests.cpp b/libc/test/src/math/sinpi_tests.cpp
new file mode 100644
index 0000000000000..514c437743f48
--- /dev/null
+++ b/libc/test/src/math/sinpi_tests.cpp
@@ -0,0 +1,88 @@
+@@ -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);
+}
diff --git a/libc/test/src/math/smoke/sinpi_test.cpp b/libc/test/src/math/smoke/sinpi_test.cpp
new file mode 100644
index 0000000000000..2b39189a365d4
--- /dev/null
+++ b/libc/test/src/math/smoke/sinpi_test.cpp
@@ -0,0 +1,47 @@
+//===-- 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));
+}
More information about the libc-commits
mailing list