]>
Commit | Line | Data |
---|---|---|
f7eac6eb | 1 | /* |
e4d82761 | 2 | * IBM Accurate Mathematical Library |
aeb25823 | 3 | * written by International Business Machines Corp. |
dff8da6b | 4 | * Copyright (C) 2001-2024 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 */ | |
e4d82761 | 23 | /* */ |
476d692e WD |
24 | /* FILES NEEDED:dla.h endian.h mydefs.h utan.h */ |
25 | /* branred.c */ | |
e4d82761 UD |
26 | /* utan.tbl */ |
27 | /* */ | |
e4d82761 | 28 | /*********************************************************************/ |
337c2708 UD |
29 | |
30 | #include <errno.h> | |
37550cb3 | 31 | #include <float.h> |
e4d82761 | 32 | #include "endian.h" |
c8b3296b | 33 | #include <dla.h> |
476d692e | 34 | #include "mydefs.h" |
1ed0291c RH |
35 | #include <math.h> |
36 | #include <math_private.h> | |
70e2ba33 | 37 | #include <fenv_private.h> |
8f5b00d3 | 38 | #include <math-underflow.h> |
38722448 | 39 | #include <libm-alias-double.h> |
804360ed | 40 | #include <fenv.h> |
15b3c029 | 41 | |
31d3cc00 UD |
42 | #ifndef SECTION |
43 | # define SECTION | |
44 | #endif | |
45 | ||
476d692e | 46 | /* tan with max ULP of ~0.619 based on random sampling. */ |
31d3cc00 UD |
47 | double |
48 | SECTION | |
527cd19c | 49 | __tan (double x) |
27ec37f1 | 50 | { |
e4d82761 UD |
51 | #include "utan.h" |
52 | #include "utan.tbl" | |
f7eac6eb | 53 | |
27ec37f1 | 54 | int ux, i, n; |
476d692e WD |
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; | |
e4d82761 | 59 | |
804360ed JM |
60 | double retval; |
61 | ||
27ec37f1 | 62 | int __branred (double, double *, double *); |
e4d82761 | 63 | |
eb92c487 | 64 | SET_RESTORE_ROUND_53BIT (FE_TONEAREST); |
804360ed | 65 | |
e4d82761 | 66 | /* x=+-INF, x=NaN */ |
27ec37f1 SP |
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 | } | |
e4d82761 | 76 | |
27ec37f1 | 77 | w = (x < 0.0) ? -x : x; |
e4d82761 UD |
78 | |
79 | /* (I) The case abs(x) <= 1.259e-8 */ | |
27ec37f1 SP |
80 | if (w <= g1.d) |
81 | { | |
d96164c3 | 82 | math_check_force_underflow_nonneg (w); |
27ec37f1 SP |
83 | retval = x; |
84 | goto ret; | |
85 | } | |
e4d82761 UD |
86 | |
87 | /* (II) The case 1.259e-8 < abs(x) <= 0.0608 */ | |
27ec37f1 SP |
88 | if (w <= g2.d) |
89 | { | |
27ec37f1 | 90 | x2 = x * x; |
e4d82761 | 91 | |
27ec37f1 SP |
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 | ||
476d692e WD |
98 | y = x + t2; |
99 | retval = y; | |
100 | /* Max ULP is 0.504. */ | |
804360ed | 101 | goto ret; |
e4d82761 UD |
102 | } |
103 | ||
27ec37f1 SP |
104 | /* (III) The case 0.0608 < abs(x) <= 0.787 */ |
105 | if (w <= g3.d) | |
106 | { | |
476d692e | 107 | i = ((int) (mfftnhf.d + 256 * w)); |
27ec37f1 SP |
108 | z = w - xfg[i][0].d; |
109 | z2 = z * z; | |
c2d94018 | 110 | s = (x < 0.0) ? -1 : 1; |
27ec37f1 SP |
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); | |
476d692e WD |
115 | y = fi + t2; |
116 | retval = (s * y); | |
117 | /* Max ULP is 0.60. */ | |
27ec37f1 SP |
118 | goto ret; |
119 | } | |
f7eac6eb | 120 | |
27ec37f1 SP |
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; | |
c2d94018 | 137 | sy = -1; |
27ec37f1 SP |
138 | } |
139 | else | |
140 | { | |
141 | ya = a; | |
142 | yya = da; | |
c2d94018 | 143 | sy = 1; |
27ec37f1 SP |
144 | } |
145 | ||
476d692e | 146 | /* (VI) The case 0.787 < abs(x) <= 25, 0 < abs(y) <= 0.0608 */ |
27ec37f1 SP |
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 | { | |
476d692e | 158 | /* -cot */ |
27ec37f1 | 159 | EADD (a, t2, b, db); |
e93c2643 | 160 | DIV2 (1.0, 0.0, b, db, c, dc, t1, t2, t3, t4); |
476d692e WD |
161 | y = c + dc; |
162 | retval = (-y); | |
163 | /* Max ULP is 0.506. */ | |
164 | goto ret; | |
27ec37f1 SP |
165 | } |
166 | else | |
167 | { | |
476d692e WD |
168 | /* tan */ |
169 | y = a + t2; | |
170 | retval = y; | |
171 | /* Max ULP is 0.506. */ | |
172 | goto ret; | |
27ec37f1 | 173 | } |
27ec37f1 SP |
174 | } |
175 | ||
176 | /* (VII) The case 0.787 < abs(x) <= 25, 0.0608 < abs(y) <= 0.787 */ | |
177 | ||
476d692e | 178 | i = ((int) (mfftnhf.d + 256 * ya)); |
27ec37f1 SP |
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); | |
476d692e WD |
189 | y = gi - t2; |
190 | retval = (-sy * y); | |
191 | /* Max ULP is 0.62. */ | |
192 | goto ret; | |
27ec37f1 SP |
193 | } |
194 | else | |
195 | { | |
196 | /* tan */ | |
197 | t2 = pz * (gi + fi) / (gi - pz); | |
476d692e WD |
198 | y = fi + t2; |
199 | retval = (sy * y); | |
200 | /* Max ULP is 0.62. */ | |
201 | goto ret; | |
27ec37f1 | 202 | } |
27ec37f1 | 203 | } |
e4d82761 UD |
204 | |
205 | /* (---) The case 25 < abs(x) <= 1e8 */ | |
27ec37f1 SP |
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; | |
c2d94018 | 227 | sy = -1; |
27ec37f1 SP |
228 | } |
229 | else | |
230 | { | |
231 | ya = a; | |
232 | yya = da; | |
c2d94018 | 233 | sy = 1; |
27ec37f1 SP |
234 | } |
235 | ||
476d692e | 236 | /* (VIII) The case 25 < abs(x) <= 1e8, 0 < abs(y) <= 0.0608 */ |
27ec37f1 SP |
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 | { | |
476d692e | 248 | /* -cot */ |
27ec37f1 | 249 | EADD (a, t2, b, db); |
e93c2643 | 250 | DIV2 (1.0, 0.0, b, db, c, dc, t1, t2, t3, t4); |
476d692e WD |
251 | y = c + dc; |
252 | retval = (-y); | |
253 | /* Max ULP is 0.506. */ | |
254 | goto ret; | |
27ec37f1 SP |
255 | } |
256 | else | |
257 | { | |
476d692e WD |
258 | /* tan */ |
259 | y = a + t2; | |
260 | retval = y; | |
261 | /* Max ULP is 0.506. */ | |
262 | goto ret; | |
27ec37f1 | 263 | } |
27ec37f1 SP |
264 | } |
265 | ||
266 | /* (IX) The case 25 < abs(x) <= 1e8, 0.0608 < abs(y) <= 0.787 */ | |
476d692e | 267 | i = ((int) (mfftnhf.d + 256 * ya)); |
27ec37f1 SP |
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); | |
476d692e WD |
278 | y = gi - t2; |
279 | retval = (-sy * y); | |
280 | /* Max ULP is 0.62. */ | |
281 | goto ret; | |
27ec37f1 SP |
282 | } |
283 | else | |
284 | { | |
285 | /* tan */ | |
286 | t2 = pz * (gi + fi) / (gi - pz); | |
476d692e WD |
287 | y = fi + t2; |
288 | retval = (sy * y); | |
289 | /* Max ULP is 0.62. */ | |
290 | goto ret; | |
27ec37f1 | 291 | } |
e4d82761 UD |
292 | } |
293 | ||
e4d82761 UD |
294 | /* (---) The case 1e8 < abs(x) < 2**1024 */ |
295 | /* Range reduction by algorithm iii */ | |
27ec37f1 SP |
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; | |
c2d94018 | 304 | sy = -1; |
27ec37f1 SP |
305 | } |
306 | else | |
307 | { | |
308 | ya = a; | |
309 | yya = da; | |
c2d94018 | 310 | sy = 1; |
27ec37f1 | 311 | } |
e4d82761 | 312 | |
476d692e | 313 | /* (X) The case 1e8 < abs(x) < 2**1024, 0 < abs(y) <= 0.0608 */ |
27ec37f1 SP |
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 | { | |
476d692e | 324 | /* -cot */ |
27ec37f1 | 325 | EADD (a, t2, b, db); |
e93c2643 | 326 | DIV2 (1.0, 0.0, b, db, c, dc, t1, t2, t3, t4); |
476d692e WD |
327 | y = c + dc; |
328 | retval = (-y); | |
329 | /* Max ULP is 0.506. */ | |
330 | goto ret; | |
27ec37f1 SP |
331 | } |
332 | else | |
333 | { | |
476d692e WD |
334 | /* tan */ |
335 | y = a + t2; | |
336 | retval = y; | |
337 | /* Max ULP is 0.507. */ | |
338 | goto ret; | |
27ec37f1 | 339 | } |
27ec37f1 | 340 | } |
e4d82761 UD |
341 | |
342 | /* (XI) The case 1e8 < abs(x) < 2**1024, 0.0608 < abs(y) <= 0.787 */ | |
476d692e | 343 | i = ((int) (mfftnhf.d + 256 * ya)); |
27ec37f1 SP |
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); | |
476d692e WD |
354 | y = gi - t2; |
355 | retval = (-sy * y); | |
356 | /* Max ULP is 0.62. */ | |
357 | goto ret; | |
27ec37f1 SP |
358 | } |
359 | else | |
360 | { | |
361 | /* tan */ | |
362 | t2 = pz * (gi + fi) / (gi - pz); | |
476d692e WD |
363 | y = fi + t2; |
364 | retval = (sy * y); | |
365 | /* Max ULP is 0.62. */ | |
366 | goto ret; | |
27ec37f1 | 367 | } |
e4d82761 | 368 | |
27ec37f1 | 369 | ret: |
804360ed JM |
370 | return retval; |
371 | } | |
e4d82761 | 372 | |
527cd19c | 373 | #ifndef __tan |
38722448 | 374 | libm_alias_double (__tan, tan) |
cccda09f | 375 | #endif |