[libc-commits] [libc] [libc][math][c23] add sinpi function (PR #129379)
via libc-commits
libc-commits at lists.llvm.org
Tue Mar 11 12:18:29 PDT 2025
https://github.com/lara254 updated https://github.com/llvm/llvm-project/pull/129379
>From b4d31d0b5c0303e25054dcb43fcd1db41684f1c7 Mon Sep 17 00:00:00 2001
From: Job Hernandez <lara2993 at proton.me>
Date: Sat, 1 Mar 2025 03:29:20 -0800
Subject: [PATCH 1/3] add draft for sinpi
---
libc/src/math/generic/sinpi.cpp | 63 +++++++++++++++++++++++++++++++++
1 file changed, 63 insertions(+)
create mode 100644 libc/src/math/generic/sinpi.cpp
diff --git a/libc/src/math/generic/sinpi.cpp b/libc/src/math/generic/sinpi.cpp
new file mode 100644
index 0000000000000..69116ced9fbc7
--- /dev/null
+++ b/libc/src/math/generic/sinpi.cpp
@@ -0,0 +1,63 @@
+//===-- double-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/sinpi.h"
+#include "sincos_eval.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/config.h"
+#include "src/__support/macros/optimization.h" // LIBC_UNLIKELY
+
+namespace LIBC_NAMESPACE_DECL {
+
+LLVM_LIBC_FUNCTION(double, sinpi, (double x)) {
+ using FPBits = typename fputil::FPBits<double>;
+ FPBits xbits(x);
+
+ uint64_t x_u = xbits.uintval();
+ uint64_t x_abs = x_u & 0x7fff'ffffU;
+ long double xd = static_cast<long double>(x);
+
+ if (LIBC_UNLIKELY(x_abs >= 0x7c00)) {
+ if (x_abs == 0x7c00) {
+ fputil::set_errno_if_required(EDOM);
+ fputil::raise_except_if_required(FE_INVALID);
+ }
+ return x + FPBits::quiet_nan().get_val();
+ } else {
+ return FPBits::zero(xbits.sign()).get_val();
+ }
+
+ if (LIBC_UNLIKELY(x_abs <= 0x3d80'0000U)) {
+ if (LIBC_UNLIKELY(x_abs < 0x33CD'01D7U)) {
+ if (LIBC_UNLIKELY(x_abs == 0U)) {
+ return x;
+ }
+ long double xdpi = xd * 0x1.921fb5444d18p1;
+ return static_cast<double>(xdpi);
+ }
+ long double xsq = xd * xd;
+
+ long double result = fputil::polyeval(xsq,0x1.921fb54442d183f07b2385653d8p1, -0x1.4abbce625bd95cdc955aeed9abcp2, 0x1.466bc6769ddfdb085486c0ff3ep1, -0x1.32d2c4a48bfd71fa9cdf60a0e4p-1, 0x1.502cbd2c72e3168ff209bc7656cp-4);
+
+ return static_cast<double>(xd * result);
+ }
+
+ long double sin_k, cos_k, sin_y, cosm1_y;
+ sincospi_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<double>(fputil::multiply_add(
+ sin_y, cos_k, fputil::multiply_add(cosm1_y, sin_k, sin_k)));
+}
+}
>From b85e0dcef49a06f3054cea2882f186c3b78a20a7 Mon Sep 17 00:00:00 2001
From: Job Hernandez <lara2993 at proton.me>
Date: Tue, 11 Mar 2025 11:27:01 -0700
Subject: [PATCH 2/3] add outline
---
libc/src/math/generic/sinpi.cpp | 46 ++++++++++++++++++++++++++-------
1 file changed, 36 insertions(+), 10 deletions(-)
diff --git a/libc/src/math/generic/sinpi.cpp b/libc/src/math/generic/sinpi.cpp
index 69116ced9fbc7..f343e0b74aa4d 100644
--- a/libc/src/math/generic/sinpi.cpp
+++ b/libc/src/math/generic/sinpi.cpp
@@ -8,6 +8,7 @@
#include "src/math/sinpi.h"
#include "sincos_eval.h"
+#include "src/__support/FPUtil/double_double.h"
#include "src/__support/FPUtil/FEnvImpl.h"
#include "src/__support/FPUtil/FPBits.h"
#include "src/__support/FPUtil/PolyEval.h"
@@ -15,18 +16,43 @@
#include "src/__support/common.h"
#include "src/__support/macros/config.h"
#include "src/__support/macros/optimization.h" // LIBC_UNLIKELY
+#include "src/math/fmul.h"
namespace LIBC_NAMESPACE_DECL {
+
+static LIB_INLINE int64_t range_reduction_sincospi(double x, double &y) {
+ double b = 32.0;
+
+ //fputil::DoubleDouble prod = fputil::exact_mult(x, b);
+ //using DoubleBits = fputil::FPBits<double>;
+ //using DoubleStorageType = typename DoubleBits::StorageType;
+
+
+ float result = fmul(x, b);
+ double res = static_cast<double>(result);
+ fputil::DoubleDouble sum = fputil::exact_add(res, -res);
+ double sum_result = sum.hi;
+ // do the exceptions here...
+
+ return static_cast<int64_t>(sum_result);
+}
+
+LIBC_INLINE void sincospi_eval(double xd, double &sin_k, double &cos_k, double &sin_k, double &cosm1_y) {
+ double y;
+ int64_t k = range_reduction_sincospi(xd, y);
+ sincospi_eval(xd, sin_k, cos_k, sin_y, cosm1_y);
+ }
+
LLVM_LIBC_FUNCTION(double, sinpi, (double x)) {
using FPBits = typename fputil::FPBits<double>;
FPBits xbits(x);
uint64_t x_u = xbits.uintval();
uint64_t x_abs = x_u & 0x7fff'ffffU;
- long double xd = static_cast<long double>(x);
if (LIBC_UNLIKELY(x_abs >= 0x7c00)) {
+ // If value is equal to infinity
if (x_abs == 0x7c00) {
fputil::set_errno_if_required(EDOM);
fputil::raise_except_if_required(FE_INVALID);
@@ -41,23 +67,23 @@ LLVM_LIBC_FUNCTION(double, sinpi, (double x)) {
if (LIBC_UNLIKELY(x_abs == 0U)) {
return x;
}
- long double xdpi = xd * 0x1.921fb5444d18p1;
- return static_cast<double>(xdpi);
+ double xdpi = xd * 0x1.921fb5444d18p1;
+ return xdpi;
+
}
- long double xsq = xd * xd;
-
- long double result = fputil::polyeval(xsq,0x1.921fb54442d183f07b2385653d8p1, -0x1.4abbce625bd95cdc955aeed9abcp2, 0x1.466bc6769ddfdb085486c0ff3ep1, -0x1.32d2c4a48bfd71fa9cdf60a0e4p-1, 0x1.502cbd2c72e3168ff209bc7656cp-4);
+ double xsq = xd * xd;
+ // todo: generate a new polynomial using double precision
+ double result = fputil::polyeval(xsq,0x1.921fb54442d183f07b2385653d8p1, -0x1.4abbce625bd95cdc955aeed9abcp2, 0x1.466bc6769ddfdb085486c0ff3ep1, -0x1.32d2c4a48bfd71fa9cdf60a0e4p-1, 0x1.502cbd2c72e3168ff209bc7656cp-4);
- return static_cast<double>(xd * result);
+ return (xd * result);
}
- long double sin_k, cos_k, sin_y, cosm1_y;
+ double sin_k, cos_k, sin_y, cosm1_y;
sincospi_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<double>(fputil::multiply_add(
- sin_y, cos_k, fputil::multiply_add(cosm1_y, sin_k, sin_k)));
+ return fputil::multiply_add(sin_y, cos_k, fputil::multiply_add(cosm1_y, sin_k, sin_k)));
}
}
>From a8dae3106a4186095ec24908ac3b0905d3e9f8a3 Mon Sep 17 00:00:00 2001
From: Job Hernandez <lara2993 at proton.me>
Date: Tue, 11 Mar 2025 12:18:14 -0700
Subject: [PATCH 3/3] fix bug in range reduction
---
libc/src/math/generic/sinpi.cpp | 8 ++++----
1 file changed, 4 insertions(+), 4 deletions(-)
diff --git a/libc/src/math/generic/sinpi.cpp b/libc/src/math/generic/sinpi.cpp
index f343e0b74aa4d..55d2a9d350e5d 100644
--- a/libc/src/math/generic/sinpi.cpp
+++ b/libc/src/math/generic/sinpi.cpp
@@ -29,13 +29,13 @@ static LIB_INLINE int64_t range_reduction_sincospi(double x, double &y) {
float result = fmul(x, b);
- double res = static_cast<double>(result);
- fputil::DoubleDouble sum = fputil::exact_add(res, -res);
- double sum_result = sum.hi;
+ double k = static_cast<double>(result);
+ fputil::DoubleDouble y = fputil::exact_add(res, -res);
+ double y = sum.hi;
// do the exceptions here...
- return static_cast<int64_t>(sum_result);
+ return static_cast<int64_t>(k);
}
LIBC_INLINE void sincospi_eval(double xd, double &sin_k, double &cos_k, double &sin_k, double &cosm1_y) {
More information about the libc-commits
mailing list