Skip to content

Commit ea93c53

Browse files
HendrikHuebnerHendrikHuebner
authored
[libc][math][c23] Implemented sinpif function correctly rounded for all rounding modes. (#97149)
This implements the sinpif function. An exhaustive test shows it's correct for all rounding modes. Issue: #94895
1 parent 7840c00 commit ea93c53

File tree

19 files changed

+459
-17
lines changed

19 files changed

+459
-17
lines changed

‎libc/config/darwin/arm/entrypoints.txt

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -230,6 +230,7 @@ set(TARGET_LIBM_ENTRYPOINTS
230230
libc.src.math.sinhf
231231
libc.src.math.sin
232232
libc.src.math.sinf
233+
libc.src.math.sinpif
233234
libc.src.math.sqrt
234235
libc.src.math.sqrtf
235236
libc.src.math.sqrtl

‎libc/config/linux/aarch64/entrypoints.txt

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -485,6 +485,7 @@ set(TARGET_LIBM_ENTRYPOINTS
485485
libc.src.math.sincosf
486486
libc.src.math.sinf
487487
libc.src.math.sinhf
488+
libc.src.math.sinpif
488489
libc.src.math.sqrt
489490
libc.src.math.sqrtf
490491
libc.src.math.sqrtl

‎libc/config/linux/riscv/entrypoints.txt

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -493,6 +493,7 @@ set(TARGET_LIBM_ENTRYPOINTS
493493
libc.src.math.sincosf
494494
libc.src.math.sinf
495495
libc.src.math.sinhf
496+
libc.src.math.sinpif
496497
libc.src.math.sqrt
497498
libc.src.math.sqrtf
498499
libc.src.math.sqrtl

‎libc/config/linux/x86_64/entrypoints.txt

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -514,6 +514,7 @@ set(TARGET_LIBM_ENTRYPOINTS
514514
libc.src.math.sincosf
515515
libc.src.math.sinf
516516
libc.src.math.sinhf
517+
libc.src.math.sinpif
517518
libc.src.math.sqrt
518519
libc.src.math.sqrtf
519520
libc.src.math.sqrtl

‎libc/docs/math/index.rst

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -330,7 +330,7 @@ Higher Math Functions
330330
+-----------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+----------------------------+
331331
| sinh | |check| | | | | | 7.12.5.5 | F.10.2.5 |
332332
+-----------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+----------------------------+
333-
| sinpi | | | | | | 7.12.4.13 | F.10.1.13 |
333+
| sinpi | |check| | | | | | 7.12.4.13 | F.10.1.13 |
334334
+-----------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+----------------------------+
335335
| sqrt | |check| | |check| | |check| | | |check| | 7.12.7.10 | F.10.4.10 |
336336
+-----------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+----------------------------+

‎libc/src/math/CMakeLists.txt

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -376,6 +376,7 @@ add_math_entrypoint_object(sincosf)
376376

377377
add_math_entrypoint_object(sin)
378378
add_math_entrypoint_object(sinf)
379+
add_math_entrypoint_object(sinpif)
379380

380381
add_math_entrypoint_object(sinh)
381382
add_math_entrypoint_object(sinhf)

‎libc/src/math/generic/CMakeLists.txt

Lines changed: 19 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -283,6 +283,25 @@ add_entrypoint_object(
283283
-O3
284284
)
285285

286+
add_entrypoint_object(
287+
sinpif
288+
SRCS
289+
sinpif.cpp
290+
HDRS
291+
../sinpif.h
292+
DEPENDS
293+
.sincosf_utils
294+
libc.src.__support.FPUtil.fenv_impl
295+
libc.src.__support.FPUtil.fp_bits
296+
libc.src.__support.FPUtil.fma
297+
libc.src.__support.FPUtil.multiply_add
298+
libc.src.__support.FPUtil.polyeval
299+
libc.src.__support.common
300+
libc.src.__support.macros.optimization
301+
COMPILE_OPTIONS
302+
-O3
303+
)
304+
286305
add_entrypoint_object(
287306
sincosf
288307
SRCS

‎libc/src/math/generic/range_reduction_fma.h

Lines changed: 0 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -12,7 +12,6 @@
1212
#include "src/__support/FPUtil/FMA.h"
1313
#include "src/__support/FPUtil/FPBits.h"
1414
#include "src/__support/FPUtil/nearest_integer.h"
15-
#include "src/__support/common.h"
1615

1716
namespace LIBC_NAMESPACE {
1817

‎libc/src/math/generic/sincosf_utils.h

Lines changed: 36 additions & 12 deletions
Original file line numberDiff line numberDiff line change
@@ -19,6 +19,7 @@
1919
using LIBC_NAMESPACE::fma::FAST_PASS_BOUND;
2020
using LIBC_NAMESPACE::fma::large_range_reduction;
2121
using LIBC_NAMESPACE::fma::small_range_reduction;
22+
2223
#else
2324
#include "range_reduction.h"
2425
// using namespace LIBC_NAMESPACE::generic;
@@ -58,18 +59,9 @@ const double SIN_K_PI_OVER_32[64] = {
5859
-0x1.917a6bc29b42cp-4,
5960
};
6061

61-
LIBC_INLINE void sincosf_eval(double xd, uint32_t x_abs, double &sin_k,
62-
double &cos_k, double &sin_y, double &cosm1_y) {
63-
int64_t k;
64-
double y;
65-
66-
if (LIBC_LIKELY(x_abs < FAST_PASS_BOUND)) {
67-
k = small_range_reduction(xd, y);
68-
} else {
69-
fputil::FPBits<float> x_bits(x_abs);
70-
k = large_range_reduction(xd, x_bits.get_exponent(), y);
71-
}
72-
62+
static LIBC_INLINE void sincosf_poly_eval(int64_t k, double y, double &sin_k,
63+
double &cos_k, double &sin_y,
64+
double &cosm1_y) {
7365
// After range reduction, k = round(x * 32 / pi) and y = (x * 32 / pi) - k.
7466
// So k is an integer and -0.5 <= y <= 0.5.
7567
// Then sin(x) = sin((k + y)*pi/32)
@@ -95,6 +87,38 @@ LIBC_INLINE void sincosf_eval(double xd, uint32_t x_abs, double &sin_k,
9587
0x1.03c1f070c2e27p-18, -0x1.55cc84bd942p-30);
9688
}
9789

90+
LIBC_INLINE void sincosf_eval(double xd, uint32_t x_abs, double &sin_k,
91+
double &cos_k, double &sin_y, double &cosm1_y) {
92+
int64_t k;
93+
double y;
94+
95+
if (LIBC_LIKELY(x_abs < FAST_PASS_BOUND)) {
96+
k = small_range_reduction(xd, y);
97+
} else {
98+
fputil::FPBits<float> x_bits(x_abs);
99+
k = large_range_reduction(xd, x_bits.get_exponent(), y);
100+
}
101+
102+
sincosf_poly_eval(k, y, sin_k, cos_k, sin_y, cosm1_y);
103+
}
104+
105+
// Return k and y, where
106+
// k = round(x * 32) and y = (x * 32) - k.
107+
// => pi * x = (k + y) * pi / 32
108+
static LIBC_INLINE int64_t range_reduction_sincospi(double x, double &y) {
109+
double kd = fputil::nearest_integer(x * 32);
110+
y = fputil::multiply_add<double>(x, 32.0, -kd);
111+
112+
return static_cast<int64_t>(kd);
113+
}
114+
115+
LIBC_INLINE void sincospif_eval(double xd, double &sin_k, double &cos_k,
116+
double &sin_y, double &cosm1_y) {
117+
double y;
118+
int64_t k = range_reduction_sincospi(xd, y);
119+
sincosf_poly_eval(k, y, sin_k, cos_k, sin_y, cosm1_y);
120+
}
121+
98122
} // namespace LIBC_NAMESPACE
99123

100124
#endif // LLVM_LIBC_SRC_MATH_GENERIC_SINCOSF_UTILS_H

‎libc/src/math/generic/sinpif.cpp

Lines changed: 111 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,111 @@
1+
//===-- Single-precision sinpif function ----------------------------------===//
2+
//
3+
// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
4+
// See https://llvm.org/LICENSE.txt for license information.
5+
// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
6+
//
7+
//===----------------------------------------------------------------------===//
8+
9+
#include "src/math/sinpif.h"
10+
#include "sincosf_utils.h"
11+
#include "src/__support/FPUtil/FEnvImpl.h"
12+
#include "src/__support/FPUtil/FPBits.h"
13+
#include "src/__support/FPUtil/PolyEval.h"
14+
#include "src/__support/FPUtil/multiply_add.h"
15+
#include "src/__support/common.h"
16+
#include "src/__support/macros/optimization.h" // LIBC_UNLIKELY
17+
18+
namespace LIBC_NAMESPACE {
19+
20+
LLVM_LIBC_FUNCTION(float, sinpif, (float x)) {
21+
using FPBits = typename fputil::FPBits<float>;
22+
FPBits xbits(x);
23+
24+
uint32_t x_u = xbits.uintval();
25+
uint32_t x_abs = x_u & 0x7fff'ffffU;
26+
double xd = static_cast<double>(x);
27+
28+
// Range reduction:
29+
// For |x| > pi/32, we perform range reduction as follows:
30+
// Find k and y such that:
31+
// x = (k + y) * 1/32
32+
// k is an integer
33+
// |y| < 0.5
34+
// For small range (|x| < 2^45 when FMA instructions are available, 2^22
35+
// otherwise), this is done by performing:
36+
// k = round(x * 32)
37+
// y = x * 32 - k
38+
//
39+
// Once k and y are computed, we then deduce the answer by the sine of sum
40+
// formula:
41+
// sin(x * pi) = sin((k + y)*pi/32)
42+
// = sin(y*pi/32) * cos(k*pi/32) + cos(y*pi/32) * sin(k*pi/32)
43+
// The values of sin(k*pi/32) and cos(k*pi/32) for k = 0..31 are precomputed
44+
// and stored using a vector of 32 doubles. Sin(y*pi/32) and cos(y*pi/32) are
45+
// computed using degree-7 and degree-6 minimax polynomials generated by
46+
// Sollya respectively.
47+
48+
// |x| <= 1/16
49+
if (LIBC_UNLIKELY(x_abs <= 0x3d80'0000U)) {
50+
51+
if (LIBC_UNLIKELY(x_abs < 0x33CD'01D7U)) {
52+
if (LIBC_UNLIKELY(x_abs == 0U)) {
53+
// For signed zeros.
54+
return x;
55+
}
56+
57+
// For very small values we can approximate sinpi(x) with x * pi
58+
// An exhaustive test shows that this is accurate for |x| < 9.546391 ×
59+
// 10-8
60+
double xdpi = xd * 0x1.921fb54442d18p1;
61+
return static_cast<float>(xdpi);
62+
}
63+
64+
// |x| < 1/16.
65+
double xsq = xd * xd;
66+
67+
// Degree-9 polynomial approximation:
68+
// sinpi(x) ~ x + a_3 x^3 + a_5 x^5 + a_7 x^7 + a_9 x^9
69+
// = x (1 + a_3 x^2 + ... + a_9 x^8)
70+
// = x * P(x^2)
71+
// generated by Sollya with the following commands:
72+
// > display = hexadecimal;
73+
// > Q = fpminimax(sin(pi * x)/x, [|0, 2, 4, 6, 8|], [|D...|], [0, 1/16]);
74+
double result = fputil::polyeval(
75+
xsq, 0x1.921fb54442d18p1, -0x1.4abbce625bbf2p2, 0x1.466bc675e116ap1,
76+
-0x1.32d2c0b62d41cp-1, 0x1.501ec4497cb7dp-4);
77+
return static_cast<float>(xd * result);
78+
}
79+
80+
// Numbers greater or equal to 2^23 are always integers or NaN
81+
if (LIBC_UNLIKELY(x_abs >= 0x4B00'0000)) {
82+
83+
// check for NaN values
84+
if (LIBC_UNLIKELY(x_abs >= 0x7f80'0000U)) {
85+
if (x_abs == 0x7f80'0000U) {
86+
fputil::set_errno_if_required(EDOM);
87+
fputil::raise_except_if_required(FE_INVALID);
88+
}
89+
90+
return x + FPBits::quiet_nan().get_val();
91+
}
92+
93+
return FPBits::zero(xbits.sign()).get_val();
94+
}
95+
96+
// Combine the results with the sine of sum formula:
97+
// sin(x * pi) = sin((k + y)*pi/32)
98+
// = sin(y*pi/32) * cos(k*pi/32) + cos(y*pi/32) * sin(k*pi/32)
99+
// = sin_y * cos_k + (1 + cosm1_y) * sin_k
100+
// = sin_y * cos_k + (cosm1_y * sin_k + sin_k)
101+
double sin_k, cos_k, sin_y, cosm1_y;
102+
sincospif_eval(xd, sin_k, cos_k, sin_y, cosm1_y);
103+
104+
if (LIBC_UNLIKELY(sin_y == 0 && sin_k == 0))
105+
return FPBits::zero(xbits.sign()).get_val();
106+
107+
return static_cast<float>(fputil::multiply_add(
108+
sin_y, cos_k, fputil::multiply_add(cosm1_y, sin_k, sin_k)));
109+
}
110+
111+
} // namespace LIBC_NAMESPACE

‎libc/src/math/sinpif.h

Lines changed: 18 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,18 @@
1+
//===-- Implementation header for sinpif ------------------------*- C++ -*-===//
2+
//
3+
// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
4+
// See https://llvm.org/LICENSE.txt for license information.
5+
// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
6+
//
7+
//===----------------------------------------------------------------------===//
8+
9+
#ifndef LLVM_LIBC_SRC_MATH_SINPIF_H
10+
#define LLVM_LIBC_SRC_MATH_SINPIF_H
11+
12+
namespace LIBC_NAMESPACE {
13+
14+
float sinpif(float x);
15+
16+
} // namespace LIBC_NAMESPACE
17+
18+
#endif // LLVM_LIBC_SRC_MATH_SINPIF_H

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

Lines changed: 16 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -44,6 +44,22 @@ add_fp_unittest(
4444
libc.src.__support.FPUtil.fp_bits
4545
)
4646

47+
add_fp_unittest(
48+
sinpif_test
49+
NEED_MPFR
50+
SUITE
51+
libc-math-unittests
52+
SRCS
53+
sinpif_test.cpp
54+
HDRS
55+
sdcomp26094.h
56+
DEPENDS
57+
libc.src.errno.errno
58+
libc.src.math.sinpif
59+
libc.src.__support.CPP.array
60+
libc.src.__support.FPUtil.fp_bits
61+
)
62+
4763
add_fp_unittest(
4864
sin_test
4965
NEED_MPFR

‎libc/test/src/math/exhaustive/CMakeLists.txt

Lines changed: 16 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -42,6 +42,22 @@ add_fp_unittest(
4242
-lpthread
4343
)
4444

45+
add_fp_unittest(
46+
sinpif_test
47+
NO_RUN_POSTBUILD
48+
NEED_MPFR
49+
SUITE
50+
libc_math_exhaustive_tests
51+
SRCS
52+
sinpif_test.cpp
53+
DEPENDS
54+
.exhaustive_test
55+
libc.src.math.sinpif
56+
libc.src.__support.FPUtil.fp_bits
57+
LINK_LIBRARIES
58+
-lpthread
59+
)
60+
4561
add_fp_unittest(
4662
cosf_test
4763
NO_RUN_POSTBUILD

‎libc/test/src/math/exhaustive/sinpif_test.cpp

Lines changed: 35 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,35 @@
1+
//===-- Exhaustive test for sinpif ----------------------------------------===//
2+
//
3+
// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
4+
// See https://llvm.org/LICENSE.txt for license information.
5+
// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
6+
//
7+
//===----------------------------------------------------------------------===//
8+
9+
#include "exhaustive_test.h"
10+
#include "mpfr.h"
11+
#include "src/math/sinpif.h"
12+
#include "utils/MPFRWrapper/MPFRUtils.h"
13+
#include <sys/types.h>
14+
15+
namespace mpfr = LIBC_NAMESPACE::testing::mpfr;
16+
17+
using LlvmLibcSinpifExhaustiveTest =
18+
LlvmLibcUnaryOpExhaustiveMathTest<float, mpfr::Operation::Sinpi,
19+
LIBC_NAMESPACE::sinpif>;
20+
21+
static constexpr uint32_t POS_START = 0x0000'0000U;
22+
static constexpr uint32_t POS_STOP = 0x7f80'0000U;
23+
24+
// Range: [0, Inf]
25+
TEST_F(LlvmLibcSinpifExhaustiveTest, PostiveRange) {
26+
test_full_range_all_roundings(POS_START, POS_STOP);
27+
}
28+
29+
// Range: [-Inf, 0]
30+
static constexpr uint32_t NEG_START = 0xb000'0000U;
31+
static constexpr uint32_t NEG_STOP = 0xff80'0000U;
32+
33+
TEST_F(LlvmLibcSinpifExhaustiveTest, NegativeRange) {
34+
test_full_range_all_roundings(NEG_START, NEG_STOP);
35+
}

0 commit comments

Comments
 (0)