]> git.ipfire.org Git - thirdparty/glibc.git/blame - sysdeps/ieee754/dbl-64/s_tan.c
projects / thirdparty / glibc.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Update copyright dates with scripts/update-copyrights
[thirdparty/glibc.git] / sysdeps / ieee754 / dbl-64 / s_tan.c
CommitLineData
f7eac6eb 1/*
e4d82761 2 * IBM Accurate Mathematical Library
aeb25823 3 * written by International Business Machines Corp.
2b778ceb 4 * Copyright (C) 2001-2021 Free Software Foundation, Inc.
f7eac6eb 5 *
e4d82761
UD		 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
cc7375ce 8 * the Free Software Foundation; either version 2.1 of the License, or
e4d82761 9 * (at your option) any later version.
f7eac6eb 10 *
e4d82761
UD		 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
c6c6dd48 14 * GNU Lesser General Public License for more details.
f7eac6eb 15 *
e4d82761 16 * You should have received a copy of the GNU Lesser General Public License
5a82c748 17 * along with this program; if not, see <https://www.gnu.org/licenses/>.
f7eac6eb 18 */
e4d82761
UD		 19/*********************************************************************/
20/* MODULE_NAME: utan.c */
21/* */
22/* FUNCTIONS: utan */
23/* tanMp */
24/* */
25/* FILES NEEDED:dla.h endian.h mpa.h mydefs.h utan.h */
26/* branred.c sincos32.c mptan.c */
27/* utan.tbl */
28/* */
29/* An ultimate tan routine. Given an IEEE double machine number x */
30/* it computes the correctly rounded (to nearest) value of tan(x). */
31/* Assumption: Machine arithmetic operations are performed in */
32/* round to nearest mode of IEEE 754 standard. */
33/* */
34/*********************************************************************/
337c2708
UD		 35
36#include <errno.h>
37550cb3 37#include <float.h>
e4d82761 38#include "endian.h"
c8b3296b 39#include <dla.h>
e4d82761
UD		 40#include "mpa.h"
41#include "MathLib.h"
1ed0291c
RH		 42#include <math.h>
43#include <math_private.h>
70e2ba33 44#include <fenv_private.h>
8f5b00d3 45#include <math-underflow.h>
38722448 46#include <libm-alias-double.h>
804360ed 47#include <fenv.h>
10e1cf6b 48#include <stap-probe.h>
15b3c029 49
31d3cc00
UD		 50#ifndef SECTION
51# define SECTION
52#endif
53
27ec37f1
SP		 54static double tanMp (double);
55void __mptan (double, mp_no *, int);
f7eac6eb 56
31d3cc00
UD		 57double
58SECTION
527cd19c 59__tan (double x)
27ec37f1 60{
e4d82761
UD		 61#include "utan.h"
62#include "utan.tbl"
f7eac6eb 63
27ec37f1
SP		 64 int ux, i, n;
65 double a, da, a2, b, db, c, dc, c1, cc1, c2, cc2, c3, cc3, fi, ffi, gi, pz,
e93c2643 66 s, sy, t, t1, t2, t3, t4, w, x2, xn, xx2, y, ya,
c5d5d574 67 yya, z0, z, zz, z2, zz2;
e4d82761 68 int p;
27ec37f1
SP		 69 number num, v;
70 mp_no mpa, mpt1, mpt2;
e4d82761 71
804360ed
JM		 72 double retval;
73
27ec37f1
SP		 74 int __branred (double, double *, double *);
75 int __mpranred (double, mp_no *, int);
e4d82761 76
eb92c487 77 SET_RESTORE_ROUND_53BIT (FE_TONEAREST);
804360ed 78
e4d82761 79 /* x=+-INF, x=NaN */
27ec37f1
SP		 80 num.d = x;
81 ux = num.i[HIGH_HALF];
82 if ((ux & 0x7ff00000) == 0x7ff00000)
83 {
84 if ((ux & 0x7fffffff) == 0x7ff00000)
85 __set_errno (EDOM);
86 retval = x - x;
87 goto ret;
88 }
e4d82761 89
27ec37f1 90 w = (x < 0.0) ? -x : x;
e4d82761
UD		 91
92 /* (I) The case abs(x) <= 1.259e-8 */
27ec37f1
SP		 93 if (w <= g1.d)
94 {
d96164c3 95 math_check_force_underflow_nonneg (w);
27ec37f1
SP		 96 retval = x;
97 goto ret;
98 }
e4d82761
UD		 99
100 /* (II) The case 1.259e-8 < abs(x) <= 0.0608 */
27ec37f1
SP		 101 if (w <= g2.d)
102 {
27ec37f1
SP		 103 /* First stage */
104 x2 = x * x;
e4d82761 105
27ec37f1
SP		 106 t2 = d9.d + x2 * d11.d;
107 t2 = d7.d + x2 * t2;
108 t2 = d5.d + x2 * t2;
109 t2 = d3.d + x2 * t2;
110 t2 *= x * x2;
111
112 if ((y = x + (t2 - u1.d * t2)) == x + (t2 + u1.d * t2))
113 {
114 retval = y;
115 goto ret;
116 }
e4d82761
UD		 117
118 /* Second stage */
27ec37f1
SP		 119 c1 = a25.d + x2 * a27.d;
120 c1 = a23.d + x2 * c1;
121 c1 = a21.d + x2 * c1;
122 c1 = a19.d + x2 * c1;
123 c1 = a17.d + x2 * c1;
124 c1 = a15.d + x2 * c1;
125 c1 *= x2;
126
e93c2643 127 EMULV (x, x, x2, xx2);
27ec37f1 128 ADD2 (a13.d, aa13.d, c1, 0.0, c2, cc2, t1, t2);
e93c2643 129 MUL2 (x2, xx2, c2, cc2, c1, cc1, t1, t2);
27ec37f1 130 ADD2 (a11.d, aa11.d, c1, cc1, c2, cc2, t1, t2);
e93c2643 131 MUL2 (x2, xx2, c2, cc2, c1, cc1, t1, t2);
27ec37f1 132 ADD2 (a9.d, aa9.d, c1, cc1, c2, cc2, t1, t2);
e93c2643 133 MUL2 (x2, xx2, c2, cc2, c1, cc1, t1, t2);
27ec37f1 134 ADD2 (a7.d, aa7.d, c1, cc1, c2, cc2, t1, t2);
e93c2643 135 MUL2 (x2, xx2, c2, cc2, c1, cc1, t1, t2);
27ec37f1 136 ADD2 (a5.d, aa5.d, c1, cc1, c2, cc2, t1, t2);
e93c2643 137 MUL2 (x2, xx2, c2, cc2, c1, cc1, t1, t2);
27ec37f1 138 ADD2 (a3.d, aa3.d, c1, cc1, c2, cc2, t1, t2);
e93c2643
VG		 139 MUL2 (x2, xx2, c2, cc2, c1, cc1, t1, t2);
140 MUL2 (x, 0.0, c1, cc1, c2, cc2, t1, t2);
27ec37f1
SP		 141 ADD2 (x, 0.0, c2, cc2, c1, cc1, t1, t2);
142 if ((y = c1 + (cc1 - u2.d * c1)) == c1 + (cc1 + u2.d * c1))
143 {
144 retval = y;
145 goto ret;
146 }
147 retval = tanMp (x);
804360ed 148 goto ret;
e4d82761
UD		 149 }
150
27ec37f1
SP		 151 /* (III) The case 0.0608 < abs(x) <= 0.787 */
152 if (w <= g3.d)
153 {
27ec37f1
SP		 154 /* First stage */
155 i = ((int) (mfftnhf.d + TWO8 * w));
156 z = w - xfg[i][0].d;
157 z2 = z * z;
c2d94018 158 s = (x < 0.0) ? -1 : 1;
27ec37f1
SP		 159 pz = z + z * z2 * (e0.d + z2 * e1.d);
160 fi = xfg[i][1].d;
161 gi = xfg[i][2].d;
162 t2 = pz * (gi + fi) / (gi - pz);
163 if ((y = fi + (t2 - fi * u3.d)) == fi + (t2 + fi * u3.d))
164 {
165 retval = (s * y);
166 goto ret;
167 }
168 t3 = (t2 < 0.0) ? -t2 : t2;
169 t4 = fi * ua3.d + t3 * ub3.d;
170 if ((y = fi + (t2 - t4)) == fi + (t2 + t4))
171 {
172 retval = (s * y);
173 goto ret;
174 }
e4d82761 175
27ec37f1
SP		 176 /* Second stage */
177 ffi = xfg[i][3].d;
178 c1 = z2 * (a7.d + z2 * (a9.d + z2 * a11.d));
e93c2643 179 EMULV (z, z, z2, zz2);
27ec37f1 180 ADD2 (a5.d, aa5.d, c1, 0.0, c2, cc2, t1, t2);
e93c2643 181 MUL2 (z2, zz2, c2, cc2, c1, cc1, t1, t2);
27ec37f1 182 ADD2 (a3.d, aa3.d, c1, cc1, c2, cc2, t1, t2);
e93c2643
VG		 183 MUL2 (z2, zz2, c2, cc2, c1, cc1, t1, t2);
184 MUL2 (z, 0.0, c1, cc1, c2, cc2, t1, t2);
27ec37f1
SP		 185 ADD2 (z, 0.0, c2, cc2, c1, cc1, t1, t2);
186
187 ADD2 (fi, ffi, c1, cc1, c2, cc2, t1, t2);
e93c2643 188 MUL2 (fi, ffi, c1, cc1, c3, cc3, t1, t2);
27ec37f1 189 SUB2 (1.0, 0.0, c3, cc3, c1, cc1, t1, t2);
e93c2643 190 DIV2 (c2, cc2, c1, cc1, c3, cc3, t1, t2, t3, t4);
27ec37f1
SP		 191
192 if ((y = c3 + (cc3 - u4.d * c3)) == c3 + (cc3 + u4.d * c3))
193 {
194 retval = (s * y);
195 goto ret;
196 }
197 retval = tanMp (x);
198 goto ret;
199 }
f7eac6eb 200
27ec37f1
SP		 201 /* (---) The case 0.787 < abs(x) <= 25 */
202 if (w <= g4.d)
203 {
204 /* Range reduction by algorithm i */
205 t = (x * hpinv.d + toint.d);
206 xn = t - toint.d;
207 v.d = t;
208 t1 = (x - xn * mp1.d) - xn * mp2.d;
209 n = v.i[LOW_HALF] & 0x00000001;
210 da = xn * mp3.d;
211 a = t1 - da;
212 da = (t1 - a) - da;
213 if (a < 0.0)
214 {
215 ya = -a;
216 yya = -da;
c2d94018 217 sy = -1;
27ec37f1
SP		 218 }
219 else
220 {
221 ya = a;
222 yya = da;
c2d94018 223 sy = 1;
27ec37f1
SP		 224 }
225
226 /* (IV),(V) The case 0.787 < abs(x) <= 25, abs(y) <= 1e-7 */
227 if (ya <= gy1.d)
228 {
229 retval = tanMp (x);
230 goto ret;
231 }
232
233 /* (VI) The case 0.787 < abs(x) <= 25, 1e-7 < abs(y) <= 0.0608 */
234 if (ya <= gy2.d)
235 {
236 a2 = a * a;
237 t2 = d9.d + a2 * d11.d;
238 t2 = d7.d + a2 * t2;
239 t2 = d5.d + a2 * t2;
240 t2 = d3.d + a2 * t2;
241 t2 = da + a * a2 * t2;
242
243 if (n)
244 {
245 /* First stage -cot */
246 EADD (a, t2, b, db);
e93c2643 247 DIV2 (1.0, 0.0, b, db, c, dc, t1, t2, t3, t4);
27ec37f1
SP		 248 if ((y = c + (dc - u6.d * c)) == c + (dc + u6.d * c))
249 {
250 retval = (-y);
251 goto ret;
252 }
253 }
254 else
255 {
256 /* First stage tan */
257 if ((y = a + (t2 - u5.d * a)) == a + (t2 + u5.d * a))
258 {
259 retval = y;
260 goto ret;
261 }
262 }
263 /* Second stage */
264 /* Range reduction by algorithm ii */
265 t = (x * hpinv.d + toint.d);
266 xn = t - toint.d;
267 v.d = t;
268 t1 = (x - xn * mp1.d) - xn * mp2.d;
269 n = v.i[LOW_HALF] & 0x00000001;
270 da = xn * pp3.d;
271 t = t1 - da;
272 da = (t1 - t) - da;
273 t1 = xn * pp4.d;
274 a = t - t1;
275 da = ((t - a) - t1) + da;
276
277 /* Second stage */
278 EADD (a, da, t1, t2);
279 a = t1;
280 da = t2;
e93c2643 281 MUL2 (a, da, a, da, x2, xx2, t1, t2);
27ec37f1
SP		 282
283 c1 = a25.d + x2 * a27.d;
284 c1 = a23.d + x2 * c1;
285 c1 = a21.d + x2 * c1;
286 c1 = a19.d + x2 * c1;
287 c1 = a17.d + x2 * c1;
288 c1 = a15.d + x2 * c1;
289 c1 *= x2;
290
291 ADD2 (a13.d, aa13.d, c1, 0.0, c2, cc2, t1, t2);
e93c2643 292 MUL2 (x2, xx2, c2, cc2, c1, cc1, t1, t2);
27ec37f1 293 ADD2 (a11.d, aa11.d, c1, cc1, c2, cc2, t1, t2);
e93c2643 294 MUL2 (x2, xx2, c2, cc2, c1, cc1, t1, t2);
27ec37f1 295 ADD2 (a9.d, aa9.d, c1, cc1, c2, cc2, t1, t2);
e93c2643 296 MUL2 (x2, xx2, c2, cc2, c1, cc1, t1, t2);
27ec37f1 297 ADD2 (a7.d, aa7.d, c1, cc1, c2, cc2, t1, t2);
e93c2643 298 MUL2 (x2, xx2, c2, cc2, c1, cc1, t1, t2);
27ec37f1 299 ADD2 (a5.d, aa5.d, c1, cc1, c2, cc2, t1, t2);
e93c2643 300 MUL2 (x2, xx2, c2, cc2, c1, cc1, t1, t2);
27ec37f1 301 ADD2 (a3.d, aa3.d, c1, cc1, c2, cc2, t1, t2);
e93c2643
VG		 302 MUL2 (x2, xx2, c2, cc2, c1, cc1, t1, t2);
303 MUL2 (a, da, c1, cc1, c2, cc2, t1, t2);
27ec37f1
SP		 304 ADD2 (a, da, c2, cc2, c1, cc1, t1, t2);
305
306 if (n)
307 {
308 /* Second stage -cot */
e93c2643 309 DIV2 (1.0, 0.0, c1, cc1, c2, cc2, t1, t2, t3, t4);
27ec37f1
SP		 310 if ((y = c2 + (cc2 - u8.d * c2)) == c2 + (cc2 + u8.d * c2))
311 {
312 retval = (-y);
313 goto ret;
314 }
315 }
316 else
317 {
318 /* Second stage tan */
319 if ((y = c1 + (cc1 - u7.d * c1)) == c1 + (cc1 + u7.d * c1))
320 {
321 retval = y;
322 goto ret;
323 }
324 }
325 retval = tanMp (x);
326 goto ret;
327 }
328
329 /* (VII) The case 0.787 < abs(x) <= 25, 0.0608 < abs(y) <= 0.787 */
330
331 /* First stage */
332 i = ((int) (mfftnhf.d + TWO8 * ya));
333 z = (z0 = (ya - xfg[i][0].d)) + yya;
334 z2 = z * z;
335 pz = z + z * z2 * (e0.d + z2 * e1.d);
336 fi = xfg[i][1].d;
337 gi = xfg[i][2].d;
338
339 if (n)
340 {
341 /* -cot */
342 t2 = pz * (fi + gi) / (fi + pz);
343 if ((y = gi - (t2 - gi * u10.d)) == gi - (t2 + gi * u10.d))
344 {
345 retval = (-sy * y);
346 goto ret;
347 }
348 t3 = (t2 < 0.0) ? -t2 : t2;
349 t4 = gi * ua10.d + t3 * ub10.d;
350 if ((y = gi - (t2 - t4)) == gi - (t2 + t4))
351 {
352 retval = (-sy * y);
353 goto ret;
354 }
355 }
356 else
357 {
358 /* tan */
359 t2 = pz * (gi + fi) / (gi - pz);
360 if ((y = fi + (t2 - fi * u9.d)) == fi + (t2 + fi * u9.d))
361 {
362 retval = (sy * y);
363 goto ret;
364 }
365 t3 = (t2 < 0.0) ? -t2 : t2;
366 t4 = fi * ua9.d + t3 * ub9.d;
367 if ((y = fi + (t2 - t4)) == fi + (t2 + t4))
368 {
369 retval = (sy * y);
370 goto ret;
371 }
372 }
e4d82761 373
27ec37f1
SP		 374 /* Second stage */
375 ffi = xfg[i][3].d;
376 EADD (z0, yya, z, zz)
e93c2643 377 MUL2 (z, zz, z, zz, z2, zz2, t1, t2);
27ec37f1
SP		 378 c1 = z2 * (a7.d + z2 * (a9.d + z2 * a11.d));
379 ADD2 (a5.d, aa5.d, c1, 0.0, c2, cc2, t1, t2);
e93c2643 380 MUL2 (z2, zz2, c2, cc2, c1, cc1, t1, t2);
27ec37f1 381 ADD2 (a3.d, aa3.d, c1, cc1, c2, cc2, t1, t2);
e93c2643
VG		 382 MUL2 (z2, zz2, c2, cc2, c1, cc1, t1, t2);
383 MUL2 (z, zz, c1, cc1, c2, cc2, t1, t2);
27ec37f1
SP		 384 ADD2 (z, zz, c2, cc2, c1, cc1, t1, t2);
385
386 ADD2 (fi, ffi, c1, cc1, c2, cc2, t1, t2);
e93c2643 387 MUL2 (fi, ffi, c1, cc1, c3, cc3, t1, t2);
27ec37f1
SP		 388 SUB2 (1.0, 0.0, c3, cc3, c1, cc1, t1, t2);
389
390 if (n)
391 {
392 /* -cot */
e93c2643 393 DIV2 (c1, cc1, c2, cc2, c3, cc3, t1, t2, t3, t4);
27ec37f1
SP		 394 if ((y = c3 + (cc3 - u12.d * c3)) == c3 + (cc3 + u12.d * c3))
395 {
396 retval = (-sy * y);
397 goto ret;
398 }
399 }
400 else
401 {
402 /* tan */
e93c2643 403 DIV2 (c2, cc2, c1, cc1, c3, cc3, t1, t2, t3, t4);
27ec37f1
SP		 404 if ((y = c3 + (cc3 - u11.d * c3)) == c3 + (cc3 + u11.d * c3))
405 {
406 retval = (sy * y);
407 goto ret;
408 }
409 }
410
411 retval = tanMp (x);
412 goto ret;
413 }
e4d82761
UD		 414
415 /* (---) The case 25 < abs(x) <= 1e8 */
27ec37f1
SP		 416 if (w <= g5.d)
417 {
418 /* Range reduction by algorithm ii */
419 t = (x * hpinv.d + toint.d);
420 xn = t - toint.d;
421 v.d = t;
422 t1 = (x - xn * mp1.d) - xn * mp2.d;
423 n = v.i[LOW_HALF] & 0x00000001;
424 da = xn * pp3.d;
425 t = t1 - da;
426 da = (t1 - t) - da;
427 t1 = xn * pp4.d;
428 a = t - t1;
429 da = ((t - a) - t1) + da;
430 EADD (a, da, t1, t2);
431 a = t1;
432 da = t2;
433 if (a < 0.0)
434 {
435 ya = -a;
436 yya = -da;
c2d94018 437 sy = -1;
27ec37f1
SP		 438 }
439 else
440 {
441 ya = a;
442 yya = da;
c2d94018 443 sy = 1;
27ec37f1
SP		 444 }
445
446 /* (+++) The case 25 < abs(x) <= 1e8, abs(y) <= 1e-7 */
447 if (ya <= gy1.d)
448 {
449 retval = tanMp (x);
450 goto ret;
451 }
452
453 /* (VIII) The case 25 < abs(x) <= 1e8, 1e-7 < abs(y) <= 0.0608 */
454 if (ya <= gy2.d)
455 {
456 a2 = a * a;
457 t2 = d9.d + a2 * d11.d;
458 t2 = d7.d + a2 * t2;
459 t2 = d5.d + a2 * t2;
460 t2 = d3.d + a2 * t2;
461 t2 = da + a * a2 * t2;
462
463 if (n)
464 {
465 /* First stage -cot */
466 EADD (a, t2, b, db);
e93c2643 467 DIV2 (1.0, 0.0, b, db, c, dc, t1, t2, t3, t4);
27ec37f1
SP		 468 if ((y = c + (dc - u14.d * c)) == c + (dc + u14.d * c))
469 {
470 retval = (-y);
471 goto ret;
472 }
473 }
474 else
475 {
476 /* First stage tan */
477 if ((y = a + (t2 - u13.d * a)) == a + (t2 + u13.d * a))
478 {
479 retval = y;
480 goto ret;
481 }
482 }
483
484 /* Second stage */
e93c2643 485 MUL2 (a, da, a, da, x2, xx2, t1, t2);
27ec37f1
SP		 486 c1 = a25.d + x2 * a27.d;
487 c1 = a23.d + x2 * c1;
488 c1 = a21.d + x2 * c1;
489 c1 = a19.d + x2 * c1;
490 c1 = a17.d + x2 * c1;
491 c1 = a15.d + x2 * c1;
492 c1 *= x2;
493
494 ADD2 (a13.d, aa13.d, c1, 0.0, c2, cc2, t1, t2);
e93c2643 495 MUL2 (x2, xx2, c2, cc2, c1, cc1, t1, t2);
27ec37f1 496 ADD2 (a11.d, aa11.d, c1, cc1, c2, cc2, t1, t2);
e93c2643 497 MUL2 (x2, xx2, c2, cc2, c1, cc1, t1, t2);
27ec37f1 498 ADD2 (a9.d, aa9.d, c1, cc1, c2, cc2, t1, t2);
e93c2643 499 MUL2 (x2, xx2, c2, cc2, c1, cc1, t1, t2);
27ec37f1 500 ADD2 (a7.d, aa7.d, c1, cc1, c2, cc2, t1, t2);
e93c2643 501 MUL2 (x2, xx2, c2, cc2, c1, cc1, t1, t2);
27ec37f1 502 ADD2 (a5.d, aa5.d, c1, cc1, c2, cc2, t1, t2);
e93c2643 503 MUL2 (x2, xx2, c2, cc2, c1, cc1, t1, t2);
27ec37f1 504 ADD2 (a3.d, aa3.d, c1, cc1, c2, cc2, t1, t2);
e93c2643
VG		 505 MUL2 (x2, xx2, c2, cc2, c1, cc1, t1, t2);
506 MUL2 (a, da, c1, cc1, c2, cc2, t1, t2);
27ec37f1
SP		 507 ADD2 (a, da, c2, cc2, c1, cc1, t1, t2);
508
509 if (n)
510 {
511 /* Second stage -cot */
e93c2643 512 DIV2 (1.0, 0.0, c1, cc1, c2, cc2, t1, t2, t3, t4);
27ec37f1
SP		 513 if ((y = c2 + (cc2 - u16.d * c2)) == c2 + (cc2 + u16.d * c2))
514 {
515 retval = (-y);
516 goto ret;
517 }
518 }
519 else
520 {
521 /* Second stage tan */
522 if ((y = c1 + (cc1 - u15.d * c1)) == c1 + (cc1 + u15.d * c1))
523 {
524 retval = (y);
525 goto ret;
526 }
527 }
528 retval = tanMp (x);
529 goto ret;
530 }
531
532 /* (IX) The case 25 < abs(x) <= 1e8, 0.0608 < abs(y) <= 0.787 */
533 /* First stage */
534 i = ((int) (mfftnhf.d + TWO8 * ya));
535 z = (z0 = (ya - xfg[i][0].d)) + yya;
536 z2 = z * z;
537 pz = z + z * z2 * (e0.d + z2 * e1.d);
538 fi = xfg[i][1].d;
539 gi = xfg[i][2].d;
540
541 if (n)
542 {
543 /* -cot */
544 t2 = pz * (fi + gi) / (fi + pz);
545 if ((y = gi - (t2 - gi * u18.d)) == gi - (t2 + gi * u18.d))
546 {
547 retval = (-sy * y);
548 goto ret;
549 }
550 t3 = (t2 < 0.0) ? -t2 : t2;
551 t4 = gi * ua18.d + t3 * ub18.d;
552 if ((y = gi - (t2 - t4)) == gi - (t2 + t4))
553 {
554 retval = (-sy * y);
555 goto ret;
556 }
557 }
558 else
559 {
560 /* tan */
561 t2 = pz * (gi + fi) / (gi - pz);
562 if ((y = fi + (t2 - fi * u17.d)) == fi + (t2 + fi * u17.d))
563 {
564 retval = (sy * y);
565 goto ret;
566 }
567 t3 = (t2 < 0.0) ? -t2 : t2;
568 t4 = fi * ua17.d + t3 * ub17.d;
569 if ((y = fi + (t2 - t4)) == fi + (t2 + t4))
570 {
571 retval = (sy * y);
572 goto ret;
573 }
574 }
e4d82761
UD		 575
576 /* Second stage */
27ec37f1
SP		 577 ffi = xfg[i][3].d;
578 EADD (z0, yya, z, zz);
e93c2643 579 MUL2 (z, zz, z, zz, z2, zz2, t1, t2);
27ec37f1
SP		 580 c1 = z2 * (a7.d + z2 * (a9.d + z2 * a11.d));
581 ADD2 (a5.d, aa5.d, c1, 0.0, c2, cc2, t1, t2);
e93c2643 582 MUL2 (z2, zz2, c2, cc2, c1, cc1, t1, t2);
27ec37f1 583 ADD2 (a3.d, aa3.d, c1, cc1, c2, cc2, t1, t2);
e93c2643
VG		 584 MUL2 (z2, zz2, c2, cc2, c1, cc1, t1, t2);
585 MUL2 (z, zz, c1, cc1, c2, cc2, t1, t2);
27ec37f1
SP		 586 ADD2 (z, zz, c2, cc2, c1, cc1, t1, t2);
587
588 ADD2 (fi, ffi, c1, cc1, c2, cc2, t1, t2);
e93c2643 589 MUL2 (fi, ffi, c1, cc1, c3, cc3, t1, t2);
27ec37f1
SP		 590 SUB2 (1.0, 0.0, c3, cc3, c1, cc1, t1, t2);
591
592 if (n)
593 {
594 /* -cot */
e93c2643 595 DIV2 (c1, cc1, c2, cc2, c3, cc3, t1, t2, t3, t4);
27ec37f1
SP		 596 if ((y = c3 + (cc3 - u20.d * c3)) == c3 + (cc3 + u20.d * c3))
597 {
598 retval = (-sy * y);
599 goto ret;
600 }
601 }
602 else
603 {
604 /* tan */
e93c2643 605 DIV2 (c2, cc2, c1, cc1, c3, cc3, t1, t2, t3, t4);
27ec37f1
SP		 606 if ((y = c3 + (cc3 - u19.d * c3)) == c3 + (cc3 + u19.d * c3))
607 {
608 retval = (sy * y);
609 goto ret;
610 }
611 }
612 retval = tanMp (x);
804360ed 613 goto ret;
e4d82761
UD		 614 }
615
e4d82761
UD		 616 /* (---) The case 1e8 < abs(x) < 2**1024 */
617 /* Range reduction by algorithm iii */
27ec37f1
SP		 618 n = (__branred (x, &a, &da)) & 0x00000001;
619 EADD (a, da, t1, t2);
620 a = t1;
621 da = t2;
622 if (a < 0.0)
623 {
624 ya = -a;
625 yya = -da;
c2d94018 626 sy = -1;
27ec37f1
SP		 627 }
628 else
629 {
630 ya = a;
631 yya = da;
c2d94018 632 sy = 1;
27ec37f1 633 }
e4d82761
UD		 634
635 /* (+++) The case 1e8 < abs(x) < 2**1024, abs(y) <= 1e-7 */
27ec37f1
SP		 636 if (ya <= gy1.d)
637 {
638 retval = tanMp (x);
639 goto ret;
640 }
e4d82761
UD		 641
642 /* (X) The case 1e8 < abs(x) < 2**1024, 1e-7 < abs(y) <= 0.0608 */
27ec37f1
SP		 643 if (ya <= gy2.d)
644 {
645 a2 = a * a;
646 t2 = d9.d + a2 * d11.d;
647 t2 = d7.d + a2 * t2;
648 t2 = d5.d + a2 * t2;
649 t2 = d3.d + a2 * t2;
650 t2 = da + a * a2 * t2;
651 if (n)
652 {
653 /* First stage -cot */
654 EADD (a, t2, b, db);
e93c2643 655 DIV2 (1.0, 0.0, b, db, c, dc, t1, t2, t3, t4);
27ec37f1
SP		 656 if ((y = c + (dc - u22.d * c)) == c + (dc + u22.d * c))
657 {
658 retval = (-y);
659 goto ret;
660 }
661 }
662 else
663 {
664 /* First stage tan */
665 if ((y = a + (t2 - u21.d * a)) == a + (t2 + u21.d * a))
666 {
667 retval = y;
668 goto ret;
669 }
670 }
671
672 /* Second stage */
673 /* Reduction by algorithm iv */
674 p = 10;
675 n = (__mpranred (x, &mpa, p)) & 0x00000001;
676 __mp_dbl (&mpa, &a, p);
677 __dbl_mp (a, &mpt1, p);
678 __sub (&mpa, &mpt1, &mpt2, p);
679 __mp_dbl (&mpt2, &da, p);
680
e93c2643 681 MUL2 (a, da, a, da, x2, xx2, t1, t2);
27ec37f1
SP		 682
683 c1 = a25.d + x2 * a27.d;
684 c1 = a23.d + x2 * c1;
685 c1 = a21.d + x2 * c1;
686 c1 = a19.d + x2 * c1;
687 c1 = a17.d + x2 * c1;
688 c1 = a15.d + x2 * c1;
689 c1 *= x2;
690
691 ADD2 (a13.d, aa13.d, c1, 0.0, c2, cc2, t1, t2);
e93c2643 692 MUL2 (x2, xx2, c2, cc2, c1, cc1, t1, t2);
27ec37f1 693 ADD2 (a11.d, aa11.d, c1, cc1, c2, cc2, t1, t2);
e93c2643 694 MUL2 (x2, xx2, c2, cc2, c1, cc1, t1, t2);
27ec37f1 695 ADD2 (a9.d, aa9.d, c1, cc1, c2, cc2, t1, t2);
e93c2643 696 MUL2 (x2, xx2, c2, cc2, c1, cc1, t1, t2);
27ec37f1 697 ADD2 (a7.d, aa7.d, c1, cc1, c2, cc2, t1, t2);
e93c2643 698 MUL2 (x2, xx2, c2, cc2, c1, cc1, t1, t2);
27ec37f1 699 ADD2 (a5.d, aa5.d, c1, cc1, c2, cc2, t1, t2);
e93c2643 700 MUL2 (x2, xx2, c2, cc2, c1, cc1, t1, t2);
27ec37f1 701 ADD2 (a3.d, aa3.d, c1, cc1, c2, cc2, t1, t2);
e93c2643
VG		 702 MUL2 (x2, xx2, c2, cc2, c1, cc1, t1, t2);
703 MUL2 (a, da, c1, cc1, c2, cc2, t1, t2);
27ec37f1
SP		 704 ADD2 (a, da, c2, cc2, c1, cc1, t1, t2);
705
706 if (n)
707 {
708 /* Second stage -cot */
e93c2643 709 DIV2 (1.0, 0.0, c1, cc1, c2, cc2, t1, t2, t3, t4);
27ec37f1
SP		 710 if ((y = c2 + (cc2 - u24.d * c2)) == c2 + (cc2 + u24.d * c2))
711 {
712 retval = (-y);
713 goto ret;
714 }
715 }
716 else
717 {
718 /* Second stage tan */
719 if ((y = c1 + (cc1 - u23.d * c1)) == c1 + (cc1 + u23.d * c1))
720 {
721 retval = y;
722 goto ret;
723 }
724 }
725 retval = tanMp (x);
726 goto ret;
727 }
e4d82761
UD		 728
729 /* (XI) The case 1e8 < abs(x) < 2**1024, 0.0608 < abs(y) <= 0.787 */
730 /* First stage */
27ec37f1
SP		 731 i = ((int) (mfftnhf.d + TWO8 * ya));
732 z = (z0 = (ya - xfg[i][0].d)) + yya;
733 z2 = z * z;
734 pz = z + z * z2 * (e0.d + z2 * e1.d);
735 fi = xfg[i][1].d;
736 gi = xfg[i][2].d;
737
738 if (n)
739 {
740 /* -cot */
741 t2 = pz * (fi + gi) / (fi + pz);
742 if ((y = gi - (t2 - gi * u26.d)) == gi - (t2 + gi * u26.d))
743 {
744 retval = (-sy * y);
745 goto ret;
746 }
747 t3 = (t2 < 0.0) ? -t2 : t2;
748 t4 = gi * ua26.d + t3 * ub26.d;
749 if ((y = gi - (t2 - t4)) == gi - (t2 + t4))
750 {
751 retval = (-sy * y);
752 goto ret;
753 }
754 }
755 else
756 {
757 /* tan */
758 t2 = pz * (gi + fi) / (gi - pz);
759 if ((y = fi + (t2 - fi * u25.d)) == fi + (t2 + fi * u25.d))
760 {
761 retval = (sy * y);
762 goto ret;
763 }
764 t3 = (t2 < 0.0) ? -t2 : t2;
765 t4 = fi * ua25.d + t3 * ub25.d;
766 if ((y = fi + (t2 - t4)) == fi + (t2 + t4))
767 {
768 retval = (sy * y);
769 goto ret;
770 }
771 }
e4d82761
UD		 772
773 /* Second stage */
774 ffi = xfg[i][3].d;
27ec37f1 775 EADD (z0, yya, z, zz);
e93c2643 776 MUL2 (z, zz, z, zz, z2, zz2, t1, t2);
27ec37f1
SP		 777 c1 = z2 * (a7.d + z2 * (a9.d + z2 * a11.d));
778 ADD2 (a5.d, aa5.d, c1, 0.0, c2, cc2, t1, t2);
e93c2643 779 MUL2 (z2, zz2, c2, cc2, c1, cc1, t1, t2);
27ec37f1 780 ADD2 (a3.d, aa3.d, c1, cc1, c2, cc2, t1, t2);
e93c2643
VG		 781 MUL2 (z2, zz2, c2, cc2, c1, cc1, t1, t2);
782 MUL2 (z, zz, c1, cc1, c2, cc2, t1, t2);
27ec37f1
SP		 783 ADD2 (z, zz, c2, cc2, c1, cc1, t1, t2);
784
785 ADD2 (fi, ffi, c1, cc1, c2, cc2, t1, t2);
e93c2643 786 MUL2 (fi, ffi, c1, cc1, c3, cc3, t1, t2);
27ec37f1
SP		 787 SUB2 (1.0, 0.0, c3, cc3, c1, cc1, t1, t2);
788
789 if (n)
790 {
791 /* -cot */
e93c2643 792 DIV2 (c1, cc1, c2, cc2, c3, cc3, t1, t2, t3, t4);
27ec37f1
SP		 793 if ((y = c3 + (cc3 - u28.d * c3)) == c3 + (cc3 + u28.d * c3))
794 {
795 retval = (-sy * y);
796 goto ret;
797 }
798 }
799 else
800 {
801 /* tan */
e93c2643 802 DIV2 (c2, cc2, c1, cc1, c3, cc3, t1, t2, t3, t4);
27ec37f1
SP		 803 if ((y = c3 + (cc3 - u27.d * c3)) == c3 + (cc3 + u27.d * c3))
804 {
805 retval = (sy * y);
806 goto ret;
807 }
808 }
809 retval = tanMp (x);
804360ed 810 goto ret;
e4d82761 811
27ec37f1 812ret:
804360ed
JM		 813 return retval;
814}
e4d82761
UD		 815
816/* multiple precision stage */
817/* Convert x to multi precision number,compute tan(x) by mptan() routine */
818/* and converts result back to double */
31d3cc00
UD		 819static double
820SECTION
27ec37f1 821tanMp (double x)
e4d82761
UD		 822{
823 int p;
824 double y;
825 mp_no mpy;
27ec37f1
SP		 826 p = 32;
827 __mptan (x, &mpy, p);
828 __mp_dbl (&mpy, &y, p);
10e1cf6b 829 LIBC_PROBE (slowtan, 2, &x, &y);
e4d82761 830 return y;
f7eac6eb 831}
e4d82761 832
527cd19c 833#ifndef __tan
38722448 834libm_alias_double (__tan, tan)
cccda09f 835#endif
A mirror of the official glibc repository