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