]>
sourceware.org Git - glibc.git/blob - sysdeps/ieee754/dbl-64/s_tan.c
2 * IBM Accurate Mathematical Library
3 * written by International Business Machines Corp.
4 * Copyright (C) 2001-2024 Free Software Foundation, Inc.
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.
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.
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/>.
19 /*********************************************************************/
20 /* MODULE_NAME: utan.c */
24 /* FILES NEEDED:dla.h endian.h mydefs.h utan.h */
28 /*********************************************************************/
36 #include <math_private.h>
37 #include <fenv_private.h>
38 #include <math-underflow.h>
39 #include <libm-alias-double.h>
46 /* tan with max ULP of ~0.619 based on random sampling. */
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
,
62 int __branred (double, double *, double *);
64 SET_RESTORE_ROUND_53BIT (FE_TONEAREST
);
68 ux
= num
.i
[HIGH_HALF
];
69 if ((ux
& 0x7ff00000) == 0x7ff00000)
71 if ((ux
& 0x7fffffff) == 0x7ff00000)
77 w
= (x
< 0.0) ? -x
: x
;
79 /* (I) The case abs(x) <= 1.259e-8 */
82 math_check_force_underflow_nonneg (w
);
87 /* (II) The case 1.259e-8 < abs(x) <= 0.0608 */
92 t2
= d9
.d
+ x2
* d11
.d
;
100 /* Max ULP is 0.504. */
104 /* (III) The case 0.0608 < abs(x) <= 0.787 */
107 i
= ((int) (mfftnhf
.d
+ 256 * w
));
110 s
= (x
< 0.0) ? -1 : 1;
111 pz
= z
+ z
* z2
* (e0
.d
+ z2
* e1
.d
);
114 t2
= pz
* (gi
+ fi
) / (gi
- pz
);
117 /* Max ULP is 0.60. */
121 /* (---) The case 0.787 < abs(x) <= 25 */
124 /* Range reduction by algorithm i */
125 t
= (x
* hpinv
.d
+ toint
.d
);
128 t1
= (x
- xn
* mp1
.d
) - xn
* mp2
.d
;
129 n
= v
.i
[LOW_HALF
] & 0x00000001;
146 /* (VI) The case 0.787 < abs(x) <= 25, 0 < abs(y) <= 0.0608 */
150 t2
= d9
.d
+ a2
* d11
.d
;
154 t2
= da
+ a
* a2
* t2
;
160 DIV2 (1.0, 0.0, b
, db
, c
, dc
, t1
, t2
, t3
, t4
);
163 /* Max ULP is 0.506. */
171 /* Max ULP is 0.506. */
176 /* (VII) The case 0.787 < abs(x) <= 25, 0.0608 < abs(y) <= 0.787 */
178 i
= ((int) (mfftnhf
.d
+ 256 * ya
));
179 z
= (z0
= (ya
- xfg
[i
][0].d
)) + yya
;
181 pz
= z
+ z
* z2
* (e0
.d
+ z2
* e1
.d
);
188 t2
= pz
* (fi
+ gi
) / (fi
+ pz
);
191 /* Max ULP is 0.62. */
197 t2
= pz
* (gi
+ fi
) / (gi
- pz
);
200 /* Max ULP is 0.62. */
205 /* (---) The case 25 < abs(x) <= 1e8 */
208 /* Range reduction by algorithm ii */
209 t
= (x
* hpinv
.d
+ toint
.d
);
212 t1
= (x
- xn
* mp1
.d
) - xn
* mp2
.d
;
213 n
= v
.i
[LOW_HALF
] & 0x00000001;
219 da
= ((t
- a
) - t1
) + da
;
220 EADD (a
, da
, t1
, t2
);
236 /* (VIII) The case 25 < abs(x) <= 1e8, 0 < abs(y) <= 0.0608 */
240 t2
= d9
.d
+ a2
* d11
.d
;
244 t2
= da
+ a
* a2
* t2
;
250 DIV2 (1.0, 0.0, b
, db
, c
, dc
, t1
, t2
, t3
, t4
);
253 /* Max ULP is 0.506. */
261 /* Max ULP is 0.506. */
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
;
270 pz
= z
+ z
* z2
* (e0
.d
+ z2
* e1
.d
);
277 t2
= pz
* (fi
+ gi
) / (fi
+ pz
);
280 /* Max ULP is 0.62. */
286 t2
= pz
* (gi
+ fi
) / (gi
- pz
);
289 /* Max ULP is 0.62. */
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
);
313 /* (X) The case 1e8 < abs(x) < 2**1024, 0 < abs(y) <= 0.0608 */
317 t2
= d9
.d
+ a2
* d11
.d
;
321 t2
= da
+ a
* a2
* t2
;
326 DIV2 (1.0, 0.0, b
, db
, c
, dc
, t1
, t2
, t3
, t4
);
329 /* Max ULP is 0.506. */
337 /* Max ULP is 0.507. */
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
;
346 pz
= z
+ z
* z2
* (e0
.d
+ z2
* e1
.d
);
353 t2
= pz
* (fi
+ gi
) / (fi
+ pz
);
356 /* Max ULP is 0.62. */
362 t2
= pz
* (gi
+ fi
) / (gi
- pz
);
365 /* Max ULP is 0.62. */
374 libm_alias_double (__tan
, tan
)
This page took 0.061589 seconds and 6 git commands to generate.