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