1*f3087befSAndrew Turner /*
2*f3087befSAndrew Turner * Double-precision scalar sincospi function.
3*f3087befSAndrew Turner *
4*f3087befSAndrew Turner * Copyright (c) 2024, Arm Limited.
5*f3087befSAndrew Turner * SPDX-License-Identifier: MIT OR Apache-2.0 WITH LLVM-exception
6*f3087befSAndrew Turner */
7*f3087befSAndrew Turner
8*f3087befSAndrew Turner #include "mathlib.h"
9*f3087befSAndrew Turner #include "math_config.h"
10*f3087befSAndrew Turner #include "test_sig.h"
11*f3087befSAndrew Turner #include "test_defs.h"
12*f3087befSAndrew Turner #include "poly_scalar_f64.h"
13*f3087befSAndrew Turner
14*f3087befSAndrew Turner /* Taylor series coefficents for sin(pi * x).
15*f3087befSAndrew Turner C2 coefficient (orginally ~=5.16771278) has been split into two parts:
16*f3087befSAndrew Turner C2_hi = 4, C2_lo = C2 - C2_hi (~=1.16771278)
17*f3087befSAndrew Turner This change in magnitude reduces floating point rounding errors.
18*f3087befSAndrew Turner C2_hi is then reintroduced after the polynomial approxmation. */
19*f3087befSAndrew Turner const static struct sincospi_data
20*f3087befSAndrew Turner {
21*f3087befSAndrew Turner double poly[10];
22*f3087befSAndrew Turner } sincospi_data = {
23*f3087befSAndrew Turner /* Taylor series coefficents for sin(pi * x). */
24*f3087befSAndrew Turner .poly = { 0x1.921fb54442d184p1, -0x1.2aef39896f94bp0, 0x1.466bc6775ab16p1,
25*f3087befSAndrew Turner -0x1.32d2cce62dc33p-1, 0x1.507834891188ep-4, -0x1.e30750a28c88ep-8,
26*f3087befSAndrew Turner 0x1.e8f48308acda4p-12, -0x1.6fc0032b3c29fp-16,
27*f3087befSAndrew Turner 0x1.af86ae521260bp-21, -0x1.012a9870eeb7dp-25 },
28*f3087befSAndrew Turner };
29*f3087befSAndrew Turner
30*f3087befSAndrew Turner /* Top 12 bits of a double (sign and exponent bits). */
31*f3087befSAndrew Turner static inline uint64_t
abstop12(double x)32*f3087befSAndrew Turner abstop12 (double x)
33*f3087befSAndrew Turner {
34*f3087befSAndrew Turner return (asuint64 (x) >> 52) & 0x7ff;
35*f3087befSAndrew Turner }
36*f3087befSAndrew Turner
37*f3087befSAndrew Turner /* Triages special cases into 4 categories:
38*f3087befSAndrew Turner -1 or +1 if iy represents half an integer
39*f3087befSAndrew Turner -1 if round(y) is odd.
40*f3087befSAndrew Turner +1 if round(y) is even.
41*f3087befSAndrew Turner -2 or +2 if iy represents and integer.
42*f3087befSAndrew Turner -2 if iy is odd.
43*f3087befSAndrew Turner +2 if iy is even.
44*f3087befSAndrew Turner The argument is the bit representation of a positive non-zero
45*f3087befSAndrew Turner finite floating-point value which is either a half or an integer. */
46*f3087befSAndrew Turner static inline int
checkint(uint64_t iy)47*f3087befSAndrew Turner checkint (uint64_t iy)
48*f3087befSAndrew Turner {
49*f3087befSAndrew Turner int e = iy >> 52;
50*f3087befSAndrew Turner if (e > 0x3ff + 52)
51*f3087befSAndrew Turner return 2;
52*f3087befSAndrew Turner if (iy & ((1ULL << (0x3ff + 52 - e)) - 1))
53*f3087befSAndrew Turner {
54*f3087befSAndrew Turner if ((iy - 1) & 2)
55*f3087befSAndrew Turner return -1;
56*f3087befSAndrew Turner else
57*f3087befSAndrew Turner return 1;
58*f3087befSAndrew Turner }
59*f3087befSAndrew Turner if (iy & (1 << (0x3ff + 52 - e)))
60*f3087befSAndrew Turner return -2;
61*f3087befSAndrew Turner return 2;
62*f3087befSAndrew Turner }
63*f3087befSAndrew Turner
64*f3087befSAndrew Turner /* Approximation for scalar double-precision sincospi(x).
65*f3087befSAndrew Turner Maximum error for sin: 3.46 ULP:
66*f3087befSAndrew Turner sincospif_sin(0x1.3d8a067cd8961p+14) got 0x1.ffe609a279008p-1 want
67*f3087befSAndrew Turner 0x1.ffe609a27900cp-1.
68*f3087befSAndrew Turner Maximum error for cos: 3.66 ULP:
69*f3087befSAndrew Turner sincospif_cos(0x1.a0ec6997557eep-24) got 0x1.ffffffffffe59p-1 want
70*f3087befSAndrew Turner 0x1.ffffffffffe5dp-1. */
71*f3087befSAndrew Turner void
arm_math_sincospi(double x,double * out_sin,double * out_cos)72*f3087befSAndrew Turner arm_math_sincospi (double x, double *out_sin, double *out_cos)
73*f3087befSAndrew Turner {
74*f3087befSAndrew Turner const struct sincospi_data *d = ptr_barrier (&sincospi_data);
75*f3087befSAndrew Turner uint64_t sign = asuint64 (x) & 0x8000000000000000;
76*f3087befSAndrew Turner
77*f3087befSAndrew Turner if (likely (abstop12 (x) < abstop12 (0x1p51)))
78*f3087befSAndrew Turner {
79*f3087befSAndrew Turner /* ax = |x| - n (range reduction into -1/2 .. 1/2). */
80*f3087befSAndrew Turner double ar_s = x - rint (x);
81*f3087befSAndrew Turner
82*f3087befSAndrew Turner /* We know that cospi(x) = sinpi(0.5 - x)
83*f3087befSAndrew Turner range reduction and offset into sinpi range -1/2 .. 1/2
84*f3087befSAndrew Turner ax = 0.5 - |x - rint(x)|. */
85*f3087befSAndrew Turner double ar_c = 0.5 - fabs (ar_s);
86*f3087befSAndrew Turner
87*f3087befSAndrew Turner /* ss = sin(pi * ax). */
88*f3087befSAndrew Turner double ar2_s = ar_s * ar_s;
89*f3087befSAndrew Turner double ar2_c = ar_c * ar_c;
90*f3087befSAndrew Turner double ar4_s = ar2_s * ar2_s;
91*f3087befSAndrew Turner double ar4_c = ar2_c * ar2_c;
92*f3087befSAndrew Turner
93*f3087befSAndrew Turner uint64_t cc_sign = ((uint64_t) llrint (x)) << 63;
94*f3087befSAndrew Turner uint64_t ss_sign = cc_sign;
95*f3087befSAndrew Turner if (ar_s == 0)
96*f3087befSAndrew Turner ss_sign = sign;
97*f3087befSAndrew Turner
98*f3087befSAndrew Turner double ss = pw_horner_9_f64 (ar2_s, ar4_s, d->poly);
99*f3087befSAndrew Turner double cc = pw_horner_9_f64 (ar2_c, ar4_c, d->poly);
100*f3087befSAndrew Turner
101*f3087befSAndrew Turner /* As all values are reduced to -1/2 .. 1/2, the result of cos(x)
102*f3087befSAndrew Turner always be positive, therefore, the sign must be introduced
103*f3087befSAndrew Turner based upon if x rounds to odd or even. For sin(x) the sign is
104*f3087befSAndrew Turner copied from x. */
105*f3087befSAndrew Turner *out_sin
106*f3087befSAndrew Turner = asdouble (asuint64 (fma (-4 * ar2_s, ar_s, ss * ar_s)) ^ ss_sign);
107*f3087befSAndrew Turner *out_cos
108*f3087befSAndrew Turner = asdouble (asuint64 (fma (-4 * ar2_c, ar_c, cc * ar_c)) ^ cc_sign);
109*f3087befSAndrew Turner }
110*f3087befSAndrew Turner else
111*f3087befSAndrew Turner {
112*f3087befSAndrew Turner /* When abs(x) > 0x1p51, the x will be either
113*f3087befSAndrew Turner - Half integer (relevant if abs(x) in [0x1p51, 0x1p52])
114*f3087befSAndrew Turner - Odd integer (relevant if abs(x) in [0x1p52, 0x1p53])
115*f3087befSAndrew Turner - Even integer (relevant if abs(x) in [0x1p53, inf])
116*f3087befSAndrew Turner - Inf or NaN. */
117*f3087befSAndrew Turner if (abstop12 (x) >= 0x7ff)
118*f3087befSAndrew Turner {
119*f3087befSAndrew Turner double inv_result = __math_invalid (x);
120*f3087befSAndrew Turner *out_sin = inv_result;
121*f3087befSAndrew Turner *out_cos = inv_result;
122*f3087befSAndrew Turner return;
123*f3087befSAndrew Turner }
124*f3087befSAndrew Turner else
125*f3087befSAndrew Turner {
126*f3087befSAndrew Turner uint64_t ax = asuint64 (x) & 0x7fffffffffffffff;
127*f3087befSAndrew Turner int m = checkint (ax);
128*f3087befSAndrew Turner /* The case where ax is half integer. */
129*f3087befSAndrew Turner if (m & 1)
130*f3087befSAndrew Turner {
131*f3087befSAndrew Turner *out_sin = sign ? -m : m;
132*f3087befSAndrew Turner *out_cos = 0;
133*f3087befSAndrew Turner return;
134*f3087befSAndrew Turner }
135*f3087befSAndrew Turner /* The case where ax is integer. */
136*f3087befSAndrew Turner else
137*f3087befSAndrew Turner {
138*f3087befSAndrew Turner *out_sin = asdouble (sign);
139*f3087befSAndrew Turner *out_cos = m >> 1;
140*f3087befSAndrew Turner return;
141*f3087befSAndrew Turner }
142*f3087befSAndrew Turner }
143*f3087befSAndrew Turner }
144*f3087befSAndrew Turner }
145*f3087befSAndrew Turner
146*f3087befSAndrew Turner #if WANT_TRIGPI_TESTS
147*f3087befSAndrew Turner TEST_DISABLE_FENV (arm_math_sincospi_sin)
148*f3087befSAndrew Turner TEST_DISABLE_FENV (arm_math_sincospi_cos)
149*f3087befSAndrew Turner TEST_ULP (arm_math_sincospi_sin, 2.96)
150*f3087befSAndrew Turner TEST_ULP (arm_math_sincospi_cos, 3.16)
151*f3087befSAndrew Turner # define SINCOS_INTERVAL(lo, hi, n) \
152*f3087befSAndrew Turner TEST_SYM_INTERVAL (arm_math_sincospi_sin, lo, hi, n) \
153*f3087befSAndrew Turner TEST_SYM_INTERVAL (arm_math_sincospi_cos, lo, hi, n)
154*f3087befSAndrew Turner SINCOS_INTERVAL (0, 0x1p-63, 10000)
155*f3087befSAndrew Turner SINCOS_INTERVAL (0x1p-63, 0.5, 50000)
156*f3087befSAndrew Turner SINCOS_INTERVAL (0.5, 0x1p51, 50000)
157*f3087befSAndrew Turner SINCOS_INTERVAL (0x1p51, inf, 10000)
158*f3087befSAndrew Turner #endif
159