1*412f47f9SXin Li /*
2*412f47f9SXin Li * Single-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 #include "mathlib.h"
9*412f47f9SXin Li #include "math_config.h"
10*412f47f9SXin Li #include "pl_sig.h"
11*412f47f9SXin Li #include "pl_test.h"
12*412f47f9SXin Li
13*412f47f9SXin Li /* Taylor series coefficents for sin(pi * x). */
14*412f47f9SXin Li #define C0 0x1.921fb6p1f
15*412f47f9SXin Li #define C1 -0x1.4abbcep2f
16*412f47f9SXin Li #define C2 0x1.466bc6p1f
17*412f47f9SXin Li #define C3 -0x1.32d2ccp-1f
18*412f47f9SXin Li #define C4 0x1.50783p-4f
19*412f47f9SXin Li #define C5 -0x1.e30750p-8f
20*412f47f9SXin Li
21*412f47f9SXin Li #define Shift 0x1.0p+23f
22*412f47f9SXin Li
23*412f47f9SXin Li /* Approximation for scalar single-precision sinpi(x) - sinpif.
24*412f47f9SXin Li Maximum error: 2.48 ULP:
25*412f47f9SXin Li sinpif(0x1.d062b6p-2) got 0x1.fa8c06p-1
26*412f47f9SXin Li want 0x1.fa8c02p-1. */
27*412f47f9SXin Li float
sinpif(float x)28*412f47f9SXin Li sinpif (float x)
29*412f47f9SXin Li {
30*412f47f9SXin Li if (isinf (x))
31*412f47f9SXin Li return __math_invalidf (x);
32*412f47f9SXin Li
33*412f47f9SXin Li float r = asfloat (asuint (x) & ~0x80000000);
34*412f47f9SXin Li uint32_t sign = asuint (x) & 0x80000000;
35*412f47f9SXin Li
36*412f47f9SXin Li /* Edge cases for when sinpif should be exactly 0. (Integers)
37*412f47f9SXin Li 0x1p23 is the limit for single precision to store any decimal places. */
38*412f47f9SXin Li if (r >= 0x1p23f)
39*412f47f9SXin Li return 0;
40*412f47f9SXin Li
41*412f47f9SXin Li int32_t m = roundf (r);
42*412f47f9SXin Li if (m == r)
43*412f47f9SXin Li return 0;
44*412f47f9SXin Li
45*412f47f9SXin Li /* For very small inputs, squaring r causes underflow.
46*412f47f9SXin Li Values below this threshold can be approximated via sinpi(x) ~= pi*x. */
47*412f47f9SXin Li if (r < 0x1p-31f)
48*412f47f9SXin Li return C0 * x;
49*412f47f9SXin Li
50*412f47f9SXin Li /* Any non-integer values >= 0x1p22f will be int + 0.5.
51*412f47f9SXin Li These values should return exactly 1 or -1. */
52*412f47f9SXin Li if (r >= 0x1p22f)
53*412f47f9SXin Li {
54*412f47f9SXin Li uint32_t iy = ((m & 1) << 31) ^ asuint (-1.0f);
55*412f47f9SXin Li return asfloat (sign ^ iy);
56*412f47f9SXin Li }
57*412f47f9SXin Li
58*412f47f9SXin Li /* n = rint(|x|). */
59*412f47f9SXin Li float n = r + Shift;
60*412f47f9SXin Li sign ^= (asuint (n) << 31);
61*412f47f9SXin Li n = n - Shift;
62*412f47f9SXin Li
63*412f47f9SXin Li /* r = |x| - n (range reduction into -1/2 .. 1/2). */
64*412f47f9SXin Li r = r - n;
65*412f47f9SXin Li
66*412f47f9SXin Li /* y = sin(pi * r). */
67*412f47f9SXin Li float r2 = r * r;
68*412f47f9SXin Li float y = fmaf (C5, r2, C4);
69*412f47f9SXin Li y = fmaf (y, r2, C3);
70*412f47f9SXin Li y = fmaf (y, r2, C2);
71*412f47f9SXin Li y = fmaf (y, r2, C1);
72*412f47f9SXin Li y = fmaf (y, r2, C0);
73*412f47f9SXin Li
74*412f47f9SXin Li /* Copy sign of x to sin(|x|). */
75*412f47f9SXin Li return asfloat (asuint (y * r) ^ sign);
76*412f47f9SXin Li }
77*412f47f9SXin Li
78*412f47f9SXin Li PL_SIG (S, F, 1, sinpi, -0.9, 0.9)
79*412f47f9SXin Li PL_TEST_ULP (sinpif, 1.99)
80*412f47f9SXin Li PL_TEST_SYM_INTERVAL (sinpif, 0, 0x1p-31, 5000)
81*412f47f9SXin Li PL_TEST_SYM_INTERVAL (sinpif, 0x1p-31, 0.5, 10000)
82*412f47f9SXin Li PL_TEST_SYM_INTERVAL (sinpif, 0.5, 0x1p22f, 10000)
83*412f47f9SXin Li PL_TEST_SYM_INTERVAL (sinpif, 0x1p22f, inf, 10000)
84