1 /*
2 * IBM Accurate Mathematical Library
3 * written by International Business Machines Corp.
4 * Copyright (C) 2001-2023 Free Software Foundation, Inc.
5 *
6 * This program is free software; you can redistribute it and/or modify
7 * it under the terms of the GNU Lesser General Public License as published by
8 * the Free Software Foundation; either version 2.1 of the License, or
9 * (at your option) any later version.
10 *
11 * This program is distributed in the hope that it will be useful,
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 * GNU Lesser General Public License for more details.
15 *
16 * You should have received a copy of the GNU Lesser General Public License
17 * along with this program; if not, see <https://www.gnu.org/licenses/>.
18 */
19 /*********************************************************************/
20 /* MODULE_NAME: utan.c */
21 /* */
22 /* FUNCTIONS: utan */
23 /* */
24 /* FILES NEEDED:dla.h endian.h mydefs.h utan.h */
25 /* branred.c */
26 /* utan.tbl */
27 /* */
28 /*********************************************************************/
29
30 #include <errno.h>
31 #include <float.h>
32 #include "endian.h"
33 #include <dla.h>
34 #include "mydefs.h"
35 #include <math.h>
36 #include <math_private.h>
37 #include <fenv_private.h>
38 #include <math-underflow.h>
39 #include <libm-alias-double.h>
40 #include <fenv.h>
41
42 #ifndef SECTION
43 # define SECTION
44 #endif
45
46 /* tan with max ULP of ~0.619 based on random sampling. */
47 double
48 SECTION
49 __tan (double x)
50 {
51 #include "utan.h"
52 #include "utan.tbl"
53
54 int ux, i, n;
55 double a, da, a2, b, db, c, dc, fi, gi, pz,
56 s, sy, t, t1, t2, t3, t4, w, x2, xn, y, ya,
57 yya, z0, z, z2;
58 mynumber num, v;
59
60 double retval;
61
62 int __branred (double, double *, double *);
63
64 SET_RESTORE_ROUND_53BIT (FE_TONEAREST);
65
66 /* x=+-INF, x=NaN */
67 num.d = x;
68 ux = num.i[HIGH_HALF];
69 if ((ux & 0x7ff00000) == 0x7ff00000)
70 {
71 if ((ux & 0x7fffffff) == 0x7ff00000)
72 __set_errno (EDOM);
73 retval = x - x;
74 goto ret;
75 }
76
77 w = (x < 0.0) ? -x : x;
78
79 /* (I) The case abs(x) <= 1.259e-8 */
80 if (w <= g1.d)
81 {
82 math_check_force_underflow_nonneg (w);
83 retval = x;
84 goto ret;
85 }
86
87 /* (II) The case 1.259e-8 < abs(x) <= 0.0608 */
88 if (w <= g2.d)
89 {
90 x2 = x * x;
91
92 t2 = d9.d + x2 * d11.d;
93 t2 = d7.d + x2 * t2;
94 t2 = d5.d + x2 * t2;
95 t2 = d3.d + x2 * t2;
96 t2 *= x * x2;
97
98 y = x + t2;
99 retval = y;
100 /* Max ULP is 0.504. */
101 goto ret;
102 }
103
104 /* (III) The case 0.0608 < abs(x) <= 0.787 */
105 if (w <= g3.d)
106 {
107 i = ((int) (mfftnhf.d + 256 * w));
108 z = w - xfg[i][0].d;
109 z2 = z * z;
110 s = (x < 0.0) ? -1 : 1;
111 pz = z + z * z2 * (e0.d + z2 * e1.d);
112 fi = xfg[i][1].d;
113 gi = xfg[i][2].d;
114 t2 = pz * (gi + fi) / (gi - pz);
115 y = fi + t2;
116 retval = (s * y);
117 /* Max ULP is 0.60. */
118 goto ret;
119 }
120
121 /* (---) The case 0.787 < abs(x) <= 25 */
122 if (w <= g4.d)
123 {
124 /* Range reduction by algorithm i */
125 t = (x * hpinv.d + toint.d);
126 xn = t - toint.d;
127 v.d = t;
128 t1 = (x - xn * mp1.d) - xn * mp2.d;
129 n = v.i[LOW_HALF] & 0x00000001;
130 da = xn * mp3.d;
131 a = t1 - da;
132 da = (t1 - a) - da;
133 if (a < 0.0)
134 {
135 ya = -a;
136 yya = -da;
137 sy = -1;
138 }
139 else
140 {
141 ya = a;
142 yya = da;
143 sy = 1;
144 }
145
146 /* (VI) The case 0.787 < abs(x) <= 25, 0 < abs(y) <= 0.0608 */
147 if (ya <= gy2.d)
148 {
149 a2 = a * a;
150 t2 = d9.d + a2 * d11.d;
151 t2 = d7.d + a2 * t2;
152 t2 = d5.d + a2 * t2;
153 t2 = d3.d + a2 * t2;
154 t2 = da + a * a2 * t2;
155
156 if (n)
157 {
158 /* -cot */
159 EADD (a, t2, b, db);
160 DIV2 (1.0, 0.0, b, db, c, dc, t1, t2, t3, t4);
161 y = c + dc;
162 retval = (-y);
163 /* Max ULP is 0.506. */
164 goto ret;
165 }
166 else
167 {
168 /* tan */
169 y = a + t2;
170 retval = y;
171 /* Max ULP is 0.506. */
172 goto ret;
173 }
174 }
175
176 /* (VII) The case 0.787 < abs(x) <= 25, 0.0608 < abs(y) <= 0.787 */
177
178 i = ((int) (mfftnhf.d + 256 * ya));
179 z = (z0 = (ya - xfg[i][0].d)) + yya;
180 z2 = z * z;
181 pz = z + z * z2 * (e0.d + z2 * e1.d);
182 fi = xfg[i][1].d;
183 gi = xfg[i][2].d;
184
185 if (n)
186 {
187 /* -cot */
188 t2 = pz * (fi + gi) / (fi + pz);
189 y = gi - t2;
190 retval = (-sy * y);
191 /* Max ULP is 0.62. */
192 goto ret;
193 }
194 else
195 {
196 /* tan */
197 t2 = pz * (gi + fi) / (gi - pz);
198 y = fi + t2;
199 retval = (sy * y);
200 /* Max ULP is 0.62. */
201 goto ret;
202 }
203 }
204
205 /* (---) The case 25 < abs(x) <= 1e8 */
206 if (w <= g5.d)
207 {
208 /* Range reduction by algorithm ii */
209 t = (x * hpinv.d + toint.d);
210 xn = t - toint.d;
211 v.d = t;
212 t1 = (x - xn * mp1.d) - xn * mp2.d;
213 n = v.i[LOW_HALF] & 0x00000001;
214 da = xn * pp3.d;
215 t = t1 - da;
216 da = (t1 - t) - da;
217 t1 = xn * pp4.d;
218 a = t - t1;
219 da = ((t - a) - t1) + da;
220 EADD (a, da, t1, t2);
221 a = t1;
222 da = t2;
223 if (a < 0.0)
224 {
225 ya = -a;
226 yya = -da;
227 sy = -1;
228 }
229 else
230 {
231 ya = a;
232 yya = da;
233 sy = 1;
234 }
235
236 /* (VIII) The case 25 < abs(x) <= 1e8, 0 < abs(y) <= 0.0608 */
237 if (ya <= gy2.d)
238 {
239 a2 = a * a;
240 t2 = d9.d + a2 * d11.d;
241 t2 = d7.d + a2 * t2;
242 t2 = d5.d + a2 * t2;
243 t2 = d3.d + a2 * t2;
244 t2 = da + a * a2 * t2;
245
246 if (n)
247 {
248 /* -cot */
249 EADD (a, t2, b, db);
250 DIV2 (1.0, 0.0, b, db, c, dc, t1, t2, t3, t4);
251 y = c + dc;
252 retval = (-y);
253 /* Max ULP is 0.506. */
254 goto ret;
255 }
256 else
257 {
258 /* tan */
259 y = a + t2;
260 retval = y;
261 /* Max ULP is 0.506. */
262 goto ret;
263 }
264 }
265
266 /* (IX) The case 25 < abs(x) <= 1e8, 0.0608 < abs(y) <= 0.787 */
267 i = ((int) (mfftnhf.d + 256 * ya));
268 z = (z0 = (ya - xfg[i][0].d)) + yya;
269 z2 = z * z;
270 pz = z + z * z2 * (e0.d + z2 * e1.d);
271 fi = xfg[i][1].d;
272 gi = xfg[i][2].d;
273
274 if (n)
275 {
276 /* -cot */
277 t2 = pz * (fi + gi) / (fi + pz);
278 y = gi - t2;
279 retval = (-sy * y);
280 /* Max ULP is 0.62. */
281 goto ret;
282 }
283 else
284 {
285 /* tan */
286 t2 = pz * (gi + fi) / (gi - pz);
287 y = fi + t2;
288 retval = (sy * y);
289 /* Max ULP is 0.62. */
290 goto ret;
291 }
292 }
293
294 /* (---) The case 1e8 < abs(x) < 2**1024 */
295 /* Range reduction by algorithm iii */
296 n = (__branred (x, &a, &da)) & 0x00000001;
297 EADD (a, da, t1, t2);
298 a = t1;
299 da = t2;
300 if (a < 0.0)
301 {
302 ya = -a;
303 yya = -da;
304 sy = -1;
305 }
306 else
307 {
308 ya = a;
309 yya = da;
310 sy = 1;
311 }
312
313 /* (X) The case 1e8 < abs(x) < 2**1024, 0 < abs(y) <= 0.0608 */
314 if (ya <= gy2.d)
315 {
316 a2 = a * a;
317 t2 = d9.d + a2 * d11.d;
318 t2 = d7.d + a2 * t2;
319 t2 = d5.d + a2 * t2;
320 t2 = d3.d + a2 * t2;
321 t2 = da + a * a2 * t2;
322 if (n)
323 {
324 /* -cot */
325 EADD (a, t2, b, db);
326 DIV2 (1.0, 0.0, b, db, c, dc, t1, t2, t3, t4);
327 y = c + dc;
328 retval = (-y);
329 /* Max ULP is 0.506. */
330 goto ret;
331 }
332 else
333 {
334 /* tan */
335 y = a + t2;
336 retval = y;
337 /* Max ULP is 0.507. */
338 goto ret;
339 }
340 }
341
342 /* (XI) The case 1e8 < abs(x) < 2**1024, 0.0608 < abs(y) <= 0.787 */
343 i = ((int) (mfftnhf.d + 256 * ya));
344 z = (z0 = (ya - xfg[i][0].d)) + yya;
345 z2 = z * z;
346 pz = z + z * z2 * (e0.d + z2 * e1.d);
347 fi = xfg[i][1].d;
348 gi = xfg[i][2].d;
349
350 if (n)
351 {
352 /* -cot */
353 t2 = pz * (fi + gi) / (fi + pz);
354 y = gi - t2;
355 retval = (-sy * y);
356 /* Max ULP is 0.62. */
357 goto ret;
358 }
359 else
360 {
361 /* tan */
362 t2 = pz * (gi + fi) / (gi - pz);
363 y = fi + t2;
364 retval = (sy * y);
365 /* Max ULP is 0.62. */
366 goto ret;
367 }
368
369 ret:
370 return retval;
371 }
372
373 #ifndef __tan
374 libm_alias_double (__tan, tan)
375 #endif