1*412f47f9SXin Li /*
2*412f47f9SXin Li * Double-precision scalar sinpi function.
3*412f47f9SXin Li *
4*412f47f9SXin Li * Copyright (c) 2023, Arm Limited.
5*412f47f9SXin Li * SPDX-License-Identifier: MIT OR Apache-2.0 WITH LLVM-exception
6*412f47f9SXin Li */
7*412f47f9SXin Li
8*412f47f9SXin Li #define _GNU_SOURCE
9*412f47f9SXin Li #include <math.h>
10*412f47f9SXin Li #include "mathlib.h"
11*412f47f9SXin Li #include "math_config.h"
12*412f47f9SXin Li #include "pl_sig.h"
13*412f47f9SXin Li #include "pl_test.h"
14*412f47f9SXin Li #include "poly_scalar_f64.h"
15*412f47f9SXin Li
16*412f47f9SXin Li /* Taylor series coefficents for sin(pi * x).
17*412f47f9SXin Li C2 coefficient (orginally ~=5.16771278) has been split into two parts:
18*412f47f9SXin Li C2_hi = 4, C2_lo = C2 - C2_hi (~=1.16771278)
19*412f47f9SXin Li This change in magnitude reduces floating point rounding errors.
20*412f47f9SXin Li C2_hi is then reintroduced after the polynomial approxmation. */
21*412f47f9SXin Li static const double poly[]
22*412f47f9SXin Li = { 0x1.921fb54442d184p1, -0x1.2aef39896f94bp0, 0x1.466bc6775ab16p1,
23*412f47f9SXin Li -0x1.32d2cce62dc33p-1, 0x1.507834891188ep-4, -0x1.e30750a28c88ep-8,
24*412f47f9SXin Li 0x1.e8f48308acda4p-12, -0x1.6fc0032b3c29fp-16, 0x1.af86ae521260bp-21,
25*412f47f9SXin Li -0x1.012a9870eeb7dp-25 };
26*412f47f9SXin Li
27*412f47f9SXin Li #define Shift 0x1.8p+52
28*412f47f9SXin Li
29*412f47f9SXin Li /* Approximation for scalar double-precision sinpi(x).
30*412f47f9SXin Li Maximum error: 3.03 ULP:
31*412f47f9SXin Li sinpi(0x1.a90da2818f8b5p+7) got 0x1.fe358f255a4b3p-1
32*412f47f9SXin Li want 0x1.fe358f255a4b6p-1. */
33*412f47f9SXin Li double
sinpi(double x)34*412f47f9SXin Li sinpi (double x)
35*412f47f9SXin Li {
36*412f47f9SXin Li if (isinf (x))
37*412f47f9SXin Li return __math_invalid (x);
38*412f47f9SXin Li
39*412f47f9SXin Li double r = asdouble (asuint64 (x) & ~0x8000000000000000);
40*412f47f9SXin Li uint64_t sign = asuint64 (x) & 0x8000000000000000;
41*412f47f9SXin Li
42*412f47f9SXin Li /* Edge cases for when sinpif should be exactly 0. (Integers)
43*412f47f9SXin Li 0x1p53 is the limit for single precision to store any decimal places. */
44*412f47f9SXin Li if (r >= 0x1p53)
45*412f47f9SXin Li return 0;
46*412f47f9SXin Li
47*412f47f9SXin Li /* If x is an integer, return 0. */
48*412f47f9SXin Li uint64_t m = (uint64_t) r;
49*412f47f9SXin Li if (r == m)
50*412f47f9SXin Li return 0;
51*412f47f9SXin Li
52*412f47f9SXin Li /* For very small inputs, squaring r causes underflow.
53*412f47f9SXin Li Values below this threshold can be approximated via sinpi(x) ≈ pi*x. */
54*412f47f9SXin Li if (r < 0x1p-63)
55*412f47f9SXin Li return M_PI * x;
56*412f47f9SXin Li
57*412f47f9SXin Li /* Any non-integer values >= 0x1x51 will be int + 0.5.
58*412f47f9SXin Li These values should return exactly 1 or -1. */
59*412f47f9SXin Li if (r >= 0x1p51)
60*412f47f9SXin Li {
61*412f47f9SXin Li uint64_t iy = ((m & 1) << 63) ^ asuint64 (1.0);
62*412f47f9SXin Li return asdouble (sign ^ iy);
63*412f47f9SXin Li }
64*412f47f9SXin Li
65*412f47f9SXin Li /* n = rint(|x|). */
66*412f47f9SXin Li double n = r + Shift;
67*412f47f9SXin Li sign ^= (asuint64 (n) << 63);
68*412f47f9SXin Li n = n - Shift;
69*412f47f9SXin Li
70*412f47f9SXin Li /* r = |x| - n (range reduction into -1/2 .. 1/2). */
71*412f47f9SXin Li r = r - n;
72*412f47f9SXin Li
73*412f47f9SXin Li /* y = sin(r). */
74*412f47f9SXin Li double r2 = r * r;
75*412f47f9SXin Li double y = horner_9_f64 (r2, poly);
76*412f47f9SXin Li y = y * r;
77*412f47f9SXin Li
78*412f47f9SXin Li /* Reintroduce C2_hi. */
79*412f47f9SXin Li y = fma (-4 * r2, r, y);
80*412f47f9SXin Li
81*412f47f9SXin Li /* Copy sign of x to sin(|x|). */
82*412f47f9SXin Li return asdouble (asuint64 (y) ^ sign);
83*412f47f9SXin Li }
84*412f47f9SXin Li
85*412f47f9SXin Li PL_SIG (S, D, 1, sinpi, -0.9, 0.9)
86*412f47f9SXin Li PL_TEST_ULP (sinpi, 2.53)
87*412f47f9SXin Li PL_TEST_SYM_INTERVAL (sinpi, 0, 0x1p-63, 5000)
88*412f47f9SXin Li PL_TEST_SYM_INTERVAL (sinpi, 0x1p-63, 0.5, 10000)
89*412f47f9SXin Li PL_TEST_SYM_INTERVAL (sinpi, 0.5, 0x1p51, 10000)
90*412f47f9SXin Li PL_TEST_SYM_INTERVAL (sinpi, 0x1p51, inf, 10000)
91