1*f3087befSAndrew Turner /*
2*f3087befSAndrew Turner * Double-precision scalar tanpi(x) 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 #include "mathlib.h"
8*f3087befSAndrew Turner #include "math_config.h"
9*f3087befSAndrew Turner #include "test_sig.h"
10*f3087befSAndrew Turner #include "test_defs.h"
11*f3087befSAndrew Turner #include "poly_scalar_f64.h"
12*f3087befSAndrew Turner
13*f3087befSAndrew Turner #define SIGN_MASK 0x8000000000000000
14*f3087befSAndrew Turner
15*f3087befSAndrew Turner const static struct tanpi_data
16*f3087befSAndrew Turner {
17*f3087befSAndrew Turner double tan_poly[14], cot_poly[9], pi, invpi;
18*f3087befSAndrew Turner } tanpi_data = {
19*f3087befSAndrew Turner /* Coefficents for tan(pi * x). */
20*f3087befSAndrew Turner .tan_poly = {
21*f3087befSAndrew Turner 0x1.4abbce625be52p3,
22*f3087befSAndrew Turner 0x1.466bc6775b0f9p5,
23*f3087befSAndrew Turner 0x1.45fff9b426f5ep7,
24*f3087befSAndrew Turner 0x1.45f4730dbca5cp9,
25*f3087befSAndrew Turner 0x1.45f3265994f85p11,
26*f3087befSAndrew Turner 0x1.45f4234b330cap13,
27*f3087befSAndrew Turner 0x1.45dca11be79ebp15,
28*f3087befSAndrew Turner 0x1.47283fc5eea69p17,
29*f3087befSAndrew Turner 0x1.3a6d958cdefaep19,
30*f3087befSAndrew Turner 0x1.927896baee627p21,
31*f3087befSAndrew Turner -0x1.89333f6acd922p19,
32*f3087befSAndrew Turner 0x1.5d4e912bb8456p27,
33*f3087befSAndrew Turner -0x1.a854d53ab6874p29,
34*f3087befSAndrew Turner 0x1.1b76de7681424p32,
35*f3087befSAndrew Turner },
36*f3087befSAndrew Turner /* Coefficents for cot(pi * x). */
37*f3087befSAndrew Turner .cot_poly = {
38*f3087befSAndrew Turner -0x1.0c152382d7366p0,
39*f3087befSAndrew Turner -0x1.60c8539c1d316p-1,
40*f3087befSAndrew Turner -0x1.4b9a2f3516354p-1,
41*f3087befSAndrew Turner -0x1.47474060b6ba8p-1,
42*f3087befSAndrew Turner -0x1.464633ad9dcb1p-1,
43*f3087befSAndrew Turner -0x1.45ff229d7edd6p-1,
44*f3087befSAndrew Turner -0x1.46d8dbf492923p-1,
45*f3087befSAndrew Turner -0x1.3873892311c6bp-1,
46*f3087befSAndrew Turner -0x1.b2f3d0ff96d73p-1,
47*f3087befSAndrew Turner },
48*f3087befSAndrew Turner .pi = 0x1.921fb54442d18p1,
49*f3087befSAndrew Turner .invpi = 0x1.45f306dc9c883p-2,
50*f3087befSAndrew Turner };
51*f3087befSAndrew Turner
52*f3087befSAndrew Turner /* Double-precision scalar tanpi(x) implementation.
53*f3087befSAndrew Turner Maximum error 2.19 ULP:
54*f3087befSAndrew Turner tanpi(0x1.68847e177a855p-2) got 0x1.fe9a0ff9bb9d7p+0
55*f3087befSAndrew Turner want 0x1.fe9a0ff9bb9d5p+0. */
56*f3087befSAndrew Turner double
arm_math_tanpi(double x)57*f3087befSAndrew Turner arm_math_tanpi (double x)
58*f3087befSAndrew Turner {
59*f3087befSAndrew Turner uint64_t xabs_12 = asuint64 (x) >> 52 & 0x7ff;
60*f3087befSAndrew Turner
61*f3087befSAndrew Turner /* x >= 0x1p54. */
62*f3087befSAndrew Turner if (unlikely (xabs_12 >= 0x434))
63*f3087befSAndrew Turner {
64*f3087befSAndrew Turner /* tanpi(+/-inf) and tanpi(+/-nan) = nan. */
65*f3087befSAndrew Turner if (unlikely (xabs_12 == 0x7ff))
66*f3087befSAndrew Turner {
67*f3087befSAndrew Turner return __math_invalid (x);
68*f3087befSAndrew Turner }
69*f3087befSAndrew Turner
70*f3087befSAndrew Turner uint64_t x_sign = asuint64 (x) & SIGN_MASK;
71*f3087befSAndrew Turner return asdouble (x_sign);
72*f3087befSAndrew Turner }
73*f3087befSAndrew Turner
74*f3087befSAndrew Turner const struct tanpi_data *d = ptr_barrier (&tanpi_data);
75*f3087befSAndrew Turner
76*f3087befSAndrew Turner double rounded = round (x);
77*f3087befSAndrew Turner if (unlikely (rounded == x))
78*f3087befSAndrew Turner {
79*f3087befSAndrew Turner /* If x == 0, return with sign. */
80*f3087befSAndrew Turner if (x == 0)
81*f3087befSAndrew Turner {
82*f3087befSAndrew Turner return x;
83*f3087befSAndrew Turner }
84*f3087befSAndrew Turner /* Otherwise, return zero with alternating sign. */
85*f3087befSAndrew Turner int64_t m = (int64_t) rounded;
86*f3087befSAndrew Turner if (x < 0)
87*f3087befSAndrew Turner {
88*f3087befSAndrew Turner return m & 1 ? 0.0 : -0.0;
89*f3087befSAndrew Turner }
90*f3087befSAndrew Turner else
91*f3087befSAndrew Turner {
92*f3087befSAndrew Turner return m & 1 ? -0.0 : 0.0;
93*f3087befSAndrew Turner }
94*f3087befSAndrew Turner }
95*f3087befSAndrew Turner
96*f3087befSAndrew Turner double x_reduced = x - rounded;
97*f3087befSAndrew Turner double abs_x_reduced = 0.5 - fabs (x_reduced);
98*f3087befSAndrew Turner
99*f3087befSAndrew Turner /* Prevent underflow exceptions. x <= 0x1p-63. */
100*f3087befSAndrew Turner if (unlikely (xabs_12 < 0x3c0))
101*f3087befSAndrew Turner {
102*f3087befSAndrew Turner return d->pi * x;
103*f3087befSAndrew Turner }
104*f3087befSAndrew Turner
105*f3087befSAndrew Turner double result, offset, scale;
106*f3087befSAndrew Turner
107*f3087befSAndrew Turner /* Test 0.25 < abs_x < 0.5 independent from abs_x_reduced. */
108*f3087befSAndrew Turner double x2 = x + x;
109*f3087befSAndrew Turner int64_t rounded_x2 = (int64_t) round (x2);
110*f3087befSAndrew Turner if (rounded_x2 & 1)
111*f3087befSAndrew Turner {
112*f3087befSAndrew Turner double r_x = abs_x_reduced;
113*f3087befSAndrew Turner
114*f3087befSAndrew Turner double r_x2 = r_x * r_x;
115*f3087befSAndrew Turner double r_x4 = r_x2 * r_x2;
116*f3087befSAndrew Turner
117*f3087befSAndrew Turner uint64_t sign = asuint64 (x_reduced) & SIGN_MASK;
118*f3087befSAndrew Turner r_x = asdouble (asuint64 (r_x) ^ sign);
119*f3087befSAndrew Turner
120*f3087befSAndrew Turner // calculate sign for half-fractional inf values
121*f3087befSAndrew Turner uint64_t is_finite = asuint64 (abs_x_reduced);
122*f3087befSAndrew Turner uint64_t is_odd = (rounded_x2 & 2) << 62;
123*f3087befSAndrew Turner uint64_t is_neg = rounded_x2 & SIGN_MASK;
124*f3087befSAndrew Turner uint64_t keep_sign = is_finite | (is_odd ^ is_neg);
125*f3087befSAndrew Turner offset = d->invpi / (keep_sign ? r_x : -r_x);
126*f3087befSAndrew Turner scale = r_x;
127*f3087befSAndrew Turner
128*f3087befSAndrew Turner result = pw_horner_8_f64 (r_x2, r_x4, d->cot_poly);
129*f3087befSAndrew Turner }
130*f3087befSAndrew Turner else
131*f3087befSAndrew Turner {
132*f3087befSAndrew Turner double r_x2 = x_reduced * x_reduced;
133*f3087befSAndrew Turner double r_x4 = r_x2 * r_x2;
134*f3087befSAndrew Turner
135*f3087befSAndrew Turner offset = d->pi * x_reduced;
136*f3087befSAndrew Turner scale = x_reduced * r_x2;
137*f3087befSAndrew Turner
138*f3087befSAndrew Turner result = pw_horner_13_f64 (r_x2, r_x4, d->tan_poly);
139*f3087befSAndrew Turner }
140*f3087befSAndrew Turner
141*f3087befSAndrew Turner return fma (scale, result, offset);
142*f3087befSAndrew Turner }
143*f3087befSAndrew Turner
144*f3087befSAndrew Turner #if WANT_EXPERIMENTAL_MATH
145*f3087befSAndrew Turner double
tanpi(double x)146*f3087befSAndrew Turner tanpi (double x)
147*f3087befSAndrew Turner {
148*f3087befSAndrew Turner return arm_math_tanpi (x);
149*f3087befSAndrew Turner }
150*f3087befSAndrew Turner #endif
151*f3087befSAndrew Turner
152*f3087befSAndrew Turner #if WANT_TRIGPI_TESTS
153*f3087befSAndrew Turner TEST_ULP (arm_math_tanpi, 1.69)
154*f3087befSAndrew Turner TEST_SYM_INTERVAL (arm_math_tanpi, 0, 0x1p-63, 50000)
155*f3087befSAndrew Turner TEST_SYM_INTERVAL (arm_math_tanpi, 0x1p-63, 0.5, 100000)
156*f3087befSAndrew Turner TEST_SYM_INTERVAL (arm_math_tanpi, 0.5, 0x1p53, 100000)
157*f3087befSAndrew Turner TEST_SYM_INTERVAL (arm_math_tanpi, 0x1p53, inf, 100000)
158*f3087befSAndrew Turner #endif
159