]>
Commit | Line | Data |
---|---|---|
f7eac6eb | 1 | /* |
e4d82761 | 2 | * IBM Accurate Mathematical Library |
aeb25823 | 3 | * written by International Business Machines Corp. |
bfff8b1b | 4 | * Copyright (C) 2001-2017 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 |
59ba27a6 | 17 | * along with this program; if not, see <http://www.gnu.org/licenses/>. |
f7eac6eb | 18 | */ |
e4d82761 UD |
19 | /****************************************************************************/ |
20 | /* */ | |
21 | /* MODULE_NAME:usncs.c */ | |
22 | /* */ | |
23 | /* FUNCTIONS: usin */ | |
24 | /* ucos */ | |
25 | /* slow */ | |
26 | /* slow1 */ | |
27 | /* slow2 */ | |
28 | /* sloww */ | |
29 | /* sloww1 */ | |
30 | /* sloww2 */ | |
31 | /* bsloww */ | |
32 | /* bsloww1 */ | |
33 | /* bsloww2 */ | |
34 | /* cslow2 */ | |
e4d82761 UD |
35 | /* FILES NEEDED: dla.h endian.h mpa.h mydefs.h usncs.h */ |
36 | /* branred.c sincos32.c dosincos.c mpa.c */ | |
37 | /* sincos.tbl */ | |
38 | /* */ | |
39 | /* An ultimate sin and routine. Given an IEEE double machine number x */ | |
40 | /* it computes the correctly rounded (to nearest) value of sin(x) or cos(x) */ | |
41 | /* Assumption: Machine arithmetic operations are performed in */ | |
42 | /* round to nearest mode of IEEE 754 standard. */ | |
43 | /* */ | |
44 | /****************************************************************************/ | |
f7eac6eb | 45 | |
f7eac6eb | 46 | |
0c59a196 | 47 | #include <errno.h> |
ad39cce0 | 48 | #include <float.h> |
e4d82761 UD |
49 | #include "endian.h" |
50 | #include "mydefs.h" | |
51 | #include "usncs.h" | |
52 | #include "MathLib.h" | |
0e9be4db | 53 | #include <math.h> |
1ed0291c | 54 | #include <math_private.h> |
804360ed | 55 | #include <fenv.h> |
e4d82761 | 56 | |
4aafb73c | 57 | /* Helper macros to compute sin of the input values. */ |
196f7f5d | 58 | #define POLYNOMIAL2(xx) ((((s5 * (xx) + s4) * (xx) + s3) * (xx) + s2) * (xx)) |
4aafb73c | 59 | |
196f7f5d | 60 | #define POLYNOMIAL(xx) (POLYNOMIAL2 (xx) + s1) |
4aafb73c SP |
61 | |
62 | /* The computed polynomial is a variation of the Taylor series expansion for | |
63 | sin(a): | |
64 | ||
65 | a - a^3/3! + a^5/5! - a^7/7! + a^9/9! + (1 - a^2) * da / 2 | |
66 | ||
67 | The constants s1, s2, s3, etc. are pre-computed values of 1/3!, 1/5! and so | |
68 | on. The result is returned to LHS and correction in COR. */ | |
8d561986 | 69 | #define TAYLOR_SIN(xx, a, da, cor) \ |
4aafb73c SP |
70 | ({ \ |
71 | double t = ((POLYNOMIAL (xx) * (a) - 0.5 * (da)) * (xx) + (da)); \ | |
72 | double res = (a) + t; \ | |
73 | (cor) = ((a) - res) + t; \ | |
74 | res; \ | |
75 | }) | |
76 | ||
77 | /* This is again a variation of the Taylor series expansion with the term | |
78 | x^3/3! expanded into the following for better accuracy: | |
79 | ||
80 | bb * x ^ 3 + 3 * aa * x * x1 * x2 + aa * x1 ^ 3 + aa * x2 ^ 3 | |
81 | ||
82 | The correction term is dx and bb + aa = -1/3! | |
83 | */ | |
84 | #define TAYLOR_SLOW(x0, dx, cor) \ | |
85 | ({ \ | |
86 | static const double th2_36 = 206158430208.0; /* 1.5*2**37 */ \ | |
87 | double xx = (x0) * (x0); \ | |
88 | double x1 = ((x0) + th2_36) - th2_36; \ | |
196f7f5d | 89 | double y = aa * x1 * x1 * x1; \ |
4aafb73c SP |
90 | double r = (x0) + y; \ |
91 | double x2 = ((x0) - x1) + (dx); \ | |
196f7f5d SP |
92 | double t = (((POLYNOMIAL2 (xx) + bb) * xx + 3.0 * aa * x1 * x2) \ |
93 | * (x0) + aa * x2 * x2 * x2 + (dx)); \ | |
4aafb73c SP |
94 | t = (((x0) - r) + y) + t; \ |
95 | double res = r + t; \ | |
96 | (cor) = (r - res) + t; \ | |
97 | res; \ | |
98 | }) | |
99 | ||
b348e1e3 SP |
100 | #define SINCOS_TABLE_LOOKUP(u, sn, ssn, cs, ccs) \ |
101 | ({ \ | |
102 | int4 k = u.i[LOW_HALF] << 2; \ | |
103 | sn = __sincostab.x[k]; \ | |
104 | ssn = __sincostab.x[k + 1]; \ | |
105 | cs = __sincostab.x[k + 2]; \ | |
106 | ccs = __sincostab.x[k + 3]; \ | |
107 | }) | |
108 | ||
31d3cc00 UD |
109 | #ifndef SECTION |
110 | # define SECTION | |
111 | #endif | |
112 | ||
af968f62 UD |
113 | extern const union |
114 | { | |
115 | int4 i[880]; | |
116 | double x[440]; | |
117 | } __sincostab attribute_hidden; | |
118 | ||
e4d82761 | 119 | static const double |
6dbe713d SP |
120 | sn3 = -1.66666666666664880952546298448555E-01, |
121 | sn5 = 8.33333214285722277379541354343671E-03, | |
122 | cs2 = 4.99999999999999999999950396842453E-01, | |
123 | cs4 = -4.16666666666664434524222570944589E-02, | |
124 | cs6 = 1.38888874007937613028114285595617E-03; | |
125 | ||
7a74607f SP |
126 | static const double t22 = 0x1.8p22; |
127 | ||
6dbe713d SP |
128 | void __dubsin (double x, double dx, double w[]); |
129 | void __docos (double x, double dx, double w[]); | |
09544cbc SP |
130 | double __mpsin (double x, double dx, bool reduce_range); |
131 | double __mpcos (double x, double dx, bool reduce_range); | |
6dbe713d SP |
132 | static double slow (double x); |
133 | static double slow1 (double x); | |
134 | static double slow2 (double x); | |
b8b7e5e6 SP |
135 | static double sloww (double x, double dx, double orig, bool shift_quadrant); |
136 | static double sloww1 (double x, double dx, double orig, bool shift_quadrant); | |
6dbe713d SP |
137 | static double sloww2 (double x, double dx, double orig, int n); |
138 | static double bsloww (double x, double dx, double orig, int n); | |
139 | static double bsloww1 (double x, double dx, double orig, int n); | |
140 | static double bsloww2 (double x, double dx, double orig, int n); | |
141 | int __branred (double x, double *a, double *aa); | |
142 | static double cslow2 (double x); | |
6dbe713d | 143 | |
758e79ec SP |
144 | /* Given a number partitioned into X and DX, this function computes the cosine |
145 | of the number by combining the sin and cos of X (as computed by a variation | |
146 | of the Taylor series) with the values looked up from the sin/cos table to | |
147 | get the result in RES and a correction value in COR. */ | |
54c86cca SP |
148 | static inline double |
149 | __always_inline | |
758e79ec | 150 | do_cos (double x, double dx, double *corp) |
392dd2de | 151 | { |
758e79ec SP |
152 | mynumber u; |
153 | ||
154 | if (x < 0) | |
155 | dx = -dx; | |
156 | ||
157 | u.x = big + fabs (x); | |
158 | x = fabs (x) - (u.x - big) + dx; | |
159 | ||
392dd2de SP |
160 | double xx, s, sn, ssn, c, cs, ccs, res, cor; |
161 | xx = x * x; | |
162 | s = x + x * xx * (sn3 + xx * sn5); | |
163 | c = xx * (cs2 + xx * (cs4 + xx * cs6)); | |
164 | SINCOS_TABLE_LOOKUP (u, sn, ssn, cs, ccs); | |
165 | cor = (ccs - s * ssn - cs * c) - sn * s; | |
166 | res = cs + cor; | |
167 | cor = (cs - res) + cor; | |
168 | *corp = cor; | |
169 | return res; | |
170 | } | |
171 | ||
758e79ec SP |
172 | /* A more precise variant of DO_COS. EPS is the adjustment to the correction |
173 | COR. */ | |
54c86cca SP |
174 | static inline double |
175 | __always_inline | |
758e79ec | 176 | do_cos_slow (double x, double dx, double eps, double *corp) |
392dd2de | 177 | { |
758e79ec SP |
178 | mynumber u; |
179 | ||
180 | if (x <= 0) | |
181 | dx = -dx; | |
182 | ||
183 | u.x = big + fabs (x); | |
184 | x = fabs (x) - (u.x - big); | |
185 | ||
392dd2de SP |
186 | double xx, y, x1, x2, e1, e2, res, cor; |
187 | double s, sn, ssn, c, cs, ccs; | |
188 | xx = x * x; | |
189 | s = x * xx * (sn3 + xx * sn5); | |
190 | c = x * dx + xx * (cs2 + xx * (cs4 + xx * cs6)); | |
191 | SINCOS_TABLE_LOOKUP (u, sn, ssn, cs, ccs); | |
192 | x1 = (x + t22) - t22; | |
193 | x2 = (x - x1) + dx; | |
194 | e1 = (sn + t22) - t22; | |
195 | e2 = (sn - e1) + ssn; | |
196 | cor = (ccs - cs * c - e1 * x2 - e2 * x) - sn * s; | |
197 | y = cs - e1 * x1; | |
198 | cor = cor + ((cs - y) - e1 * x1); | |
199 | res = y + cor; | |
200 | cor = (y - res) + cor; | |
3459931a | 201 | cor = 1.0005 * cor + __copysign (eps, cor); |
392dd2de SP |
202 | *corp = cor; |
203 | return res; | |
204 | } | |
205 | ||
758e79ec SP |
206 | /* Given a number partitioned into X and DX, this function computes the sine of |
207 | the number by combining the sin and cos of X (as computed by a variation of | |
208 | the Taylor series) with the values looked up from the sin/cos table to get | |
209 | the result in RES and a correction value in COR. */ | |
54c86cca SP |
210 | static inline double |
211 | __always_inline | |
758e79ec | 212 | do_sin (double x, double dx, double *corp) |
392dd2de | 213 | { |
758e79ec SP |
214 | mynumber u; |
215 | ||
216 | if (x <= 0) | |
217 | dx = -dx; | |
218 | u.x = big + fabs (x); | |
219 | x = fabs (x) - (u.x - big); | |
220 | ||
392dd2de SP |
221 | double xx, s, sn, ssn, c, cs, ccs, cor, res; |
222 | xx = x * x; | |
223 | s = x + (dx + x * xx * (sn3 + xx * sn5)); | |
224 | c = x * dx + xx * (cs2 + xx * (cs4 + xx * cs6)); | |
225 | SINCOS_TABLE_LOOKUP (u, sn, ssn, cs, ccs); | |
226 | cor = (ssn + s * ccs - sn * c) + cs * s; | |
227 | res = sn + cor; | |
228 | cor = (sn - res) + cor; | |
229 | *corp = cor; | |
230 | return res; | |
231 | } | |
232 | ||
758e79ec SP |
233 | /* A more precise variant of DO_SIN. EPS is the adjustment to the correction |
234 | COR. */ | |
54c86cca SP |
235 | static inline double |
236 | __always_inline | |
758e79ec | 237 | do_sin_slow (double x, double dx, double eps, double *corp) |
392dd2de | 238 | { |
758e79ec SP |
239 | mynumber u; |
240 | ||
241 | if (x <= 0) | |
242 | dx = -dx; | |
243 | u.x = big + fabs (x); | |
244 | x = fabs (x) - (u.x - big); | |
245 | ||
392dd2de SP |
246 | double xx, y, x1, x2, c1, c2, res, cor; |
247 | double s, sn, ssn, c, cs, ccs; | |
248 | xx = x * x; | |
249 | s = x * xx * (sn3 + xx * sn5); | |
250 | c = xx * (cs2 + xx * (cs4 + xx * cs6)); | |
251 | SINCOS_TABLE_LOOKUP (u, sn, ssn, cs, ccs); | |
252 | x1 = (x + t22) - t22; | |
253 | x2 = (x - x1) + dx; | |
254 | c1 = (cs + t22) - t22; | |
255 | c2 = (cs - c1) + ccs; | |
256 | cor = (ssn + s * ccs + cs * s + c2 * x + c1 * x2 - sn * x * dx) - sn * c; | |
257 | y = sn + c1 * x1; | |
258 | cor = cor + ((sn - y) + c1 * x1); | |
259 | res = y + cor; | |
260 | cor = (y - res) + cor; | |
3459931a | 261 | cor = 1.0005 * cor + __copysign (eps, cor); |
392dd2de SP |
262 | *corp = cor; |
263 | return res; | |
264 | } | |
265 | ||
36ee03e6 SP |
266 | /* Reduce range of X and compute sin of a + da. When SHIFT_QUADRANT is true, |
267 | the routine returns the cosine of a + da by rotating the quadrant once and | |
268 | computing the sine of the result. */ | |
6cce25f8 SP |
269 | static inline double |
270 | __always_inline | |
ead1ef37 | 271 | reduce_and_compute (double x, bool shift_quadrant) |
6cce25f8 | 272 | { |
975195e4 | 273 | double retval = 0, a, da; |
6cce25f8 | 274 | unsigned int n = __branred (x, &a, &da); |
ead1ef37 | 275 | int4 k = (n + shift_quadrant) % 4; |
6cce25f8 SP |
276 | switch (k) |
277 | { | |
32efd690 SP |
278 | case 2: |
279 | a = -a; | |
280 | da = -da; | |
1a822c61 | 281 | /* Fall through. */ |
32efd690 SP |
282 | case 0: |
283 | if (a * a < 0.01588) | |
284 | retval = bsloww (a, da, x, n); | |
285 | else | |
286 | retval = bsloww1 (a, da, x, n); | |
287 | break; | |
6cce25f8 | 288 | |
32efd690 SP |
289 | case 1: |
290 | case 3: | |
291 | retval = bsloww2 (a, da, x, n); | |
292 | break; | |
6cce25f8 SP |
293 | } |
294 | return retval; | |
295 | } | |
296 | ||
b3004556 SP |
297 | static inline int4 |
298 | __always_inline | |
299 | reduce_sincos_1 (double x, double *a, double *da) | |
300 | { | |
301 | mynumber v; | |
302 | ||
303 | double t = (x * hpinv + toint); | |
304 | double xn = t - toint; | |
305 | v.x = t; | |
306 | double y = (x - xn * mp1) - xn * mp2; | |
307 | int4 n = v.i[LOW_HALF] & 3; | |
308 | double db = xn * mp3; | |
309 | double b = y - db; | |
310 | db = (y - b) - db; | |
311 | ||
312 | *a = b; | |
313 | *da = db; | |
314 | ||
315 | return n; | |
316 | } | |
317 | ||
36ee03e6 SP |
318 | /* Compute sin (A + DA). cos can be computed by passing SHIFT_QUADRANT as |
319 | true, which results in shifting the quadrant N clockwise. */ | |
b3004556 SP |
320 | static double |
321 | __always_inline | |
b8b7e5e6 | 322 | do_sincos_1 (double a, double da, double x, int4 n, bool shift_quadrant) |
b3004556 | 323 | { |
758e79ec | 324 | double xx, retval, res, cor; |
b3004556 SP |
325 | double eps = fabs (x) * 1.2e-30; |
326 | ||
b8b7e5e6 | 327 | int k1 = (n + shift_quadrant) & 3; |
b3004556 SP |
328 | switch (k1) |
329 | { /* quarter of unit circle */ | |
330 | case 2: | |
331 | a = -a; | |
332 | da = -da; | |
1a822c61 | 333 | /* Fall through. */ |
b3004556 SP |
334 | case 0: |
335 | xx = a * a; | |
336 | if (xx < 0.01588) | |
337 | { | |
338 | /* Taylor series. */ | |
339 | res = TAYLOR_SIN (xx, a, da, cor); | |
3459931a | 340 | cor = 1.02 * cor + __copysign (eps, cor); |
b8b7e5e6 | 341 | retval = (res == res + cor) ? res : sloww (a, da, x, shift_quadrant); |
b3004556 SP |
342 | } |
343 | else | |
344 | { | |
758e79ec | 345 | res = do_sin (a, da, &cor); |
3459931a | 346 | cor = 1.035 * cor + __copysign (eps, cor); |
2bf49970 | 347 | retval = ((res == res + cor) ? __copysign (res, a) |
b8b7e5e6 | 348 | : sloww1 (a, da, x, shift_quadrant)); |
b3004556 SP |
349 | } |
350 | break; | |
351 | ||
352 | case 1: | |
353 | case 3: | |
758e79ec | 354 | res = do_cos (a, da, &cor); |
3459931a | 355 | cor = 1.025 * cor + __copysign (eps, cor); |
ba4e6884 | 356 | retval = ((res == res + cor) ? ((n & 2) ? -res : res) |
758e79ec SP |
357 | : sloww2 (a, da, x, n)); |
358 | break; | |
b3004556 SP |
359 | } |
360 | ||
361 | return retval; | |
362 | } | |
363 | ||
f7953c44 SP |
364 | static inline int4 |
365 | __always_inline | |
366 | reduce_sincos_2 (double x, double *a, double *da) | |
367 | { | |
368 | mynumber v; | |
369 | ||
370 | double t = (x * hpinv + toint); | |
371 | double xn = t - toint; | |
372 | v.x = t; | |
373 | double xn1 = (xn + 8.0e22) - 8.0e22; | |
374 | double xn2 = xn - xn1; | |
375 | double y = ((((x - xn1 * mp1) - xn1 * mp2) - xn2 * mp1) - xn2 * mp2); | |
376 | int4 n = v.i[LOW_HALF] & 3; | |
377 | double db = xn1 * pp3; | |
378 | t = y - db; | |
379 | db = (y - t) - db; | |
380 | db = (db - xn2 * pp3) - xn * pp4; | |
381 | double b = t + db; | |
382 | db = (t - b) + db; | |
383 | ||
384 | *a = b; | |
385 | *da = db; | |
386 | ||
387 | return n; | |
388 | } | |
389 | ||
36ee03e6 SP |
390 | /* Compute sin (A + DA). cos can be computed by passing SHIFT_QUADRANT as |
391 | true, which results in shifting the quadrant N clockwise. */ | |
f7953c44 SP |
392 | static double |
393 | __always_inline | |
b8b7e5e6 | 394 | do_sincos_2 (double a, double da, double x, int4 n, bool shift_quadrant) |
f7953c44 SP |
395 | { |
396 | double res, retval, cor, xx; | |
f7953c44 SP |
397 | |
398 | double eps = 1.0e-24; | |
399 | ||
b8b7e5e6 | 400 | int4 k = (n + shift_quadrant) & 3; |
f7953c44 SP |
401 | |
402 | switch (k) | |
403 | { | |
404 | case 2: | |
405 | a = -a; | |
406 | da = -da; | |
407 | /* Fall through. */ | |
408 | case 0: | |
409 | xx = a * a; | |
410 | if (xx < 0.01588) | |
411 | { | |
412 | /* Taylor series. */ | |
413 | res = TAYLOR_SIN (xx, a, da, cor); | |
3459931a | 414 | cor = 1.02 * cor + __copysign (eps, cor); |
f7953c44 SP |
415 | retval = (res == res + cor) ? res : bsloww (a, da, x, n); |
416 | } | |
417 | else | |
418 | { | |
758e79ec | 419 | res = do_sin (a, da, &cor); |
3459931a | 420 | cor = 1.035 * cor + __copysign (eps, cor); |
2bf49970 | 421 | retval = ((res == res + cor) ? __copysign (res, a) |
f7953c44 SP |
422 | : bsloww1 (a, da, x, n)); |
423 | } | |
424 | break; | |
425 | ||
426 | case 1: | |
427 | case 3: | |
758e79ec | 428 | res = do_cos (a, da, &cor); |
3459931a | 429 | cor = 1.025 * cor + __copysign (eps, cor); |
758e79ec SP |
430 | retval = ((res == res + cor) ? ((n & 2) ? -res : res) |
431 | : bsloww2 (a, da, x, n)); | |
432 | break; | |
f7953c44 SP |
433 | } |
434 | ||
435 | return retval; | |
436 | } | |
437 | ||
e4d82761 UD |
438 | /*******************************************************************/ |
439 | /* An ultimate sin routine. Given an IEEE double machine number x */ | |
440 | /* it computes the correctly rounded (to nearest) value of sin(x) */ | |
441 | /*******************************************************************/ | |
463ac90d SP |
442 | #ifdef IN_SINCOS |
443 | static double | |
444 | #else | |
31d3cc00 UD |
445 | double |
446 | SECTION | |
463ac90d | 447 | #endif |
6dbe713d SP |
448 | __sin (double x) |
449 | { | |
25e440c6 | 450 | double xx, res, t, cor; |
b3004556 SP |
451 | mynumber u; |
452 | int4 k, m; | |
6dbe713d SP |
453 | double retval = 0; |
454 | ||
463ac90d | 455 | #ifndef IN_SINCOS |
6dbe713d | 456 | SET_RESTORE_ROUND_53BIT (FE_TONEAREST); |
463ac90d | 457 | #endif |
6dbe713d SP |
458 | |
459 | u.x = x; | |
460 | m = u.i[HIGH_HALF]; | |
461 | k = 0x7fffffff & m; /* no sign */ | |
462 | if (k < 0x3e500000) /* if x->0 =>sin(x)=x */ | |
ad39cce0 | 463 | { |
d96164c3 | 464 | math_check_force_underflow (x); |
ad39cce0 JM |
465 | retval = x; |
466 | } | |
e4d82761 | 467 | /*---------------------------- 2^-26 < |x|< 0.25 ----------------------*/ |
6dbe713d SP |
468 | else if (k < 0x3fd00000) |
469 | { | |
470 | xx = x * x; | |
4aafb73c SP |
471 | /* Taylor series. */ |
472 | t = POLYNOMIAL (xx) * (xx * x); | |
6dbe713d SP |
473 | res = x + t; |
474 | cor = (x - res) + t; | |
475 | retval = (res == res + 1.07 * cor) ? res : slow (x); | |
6dbe713d | 476 | } /* else if (k < 0x3fd00000) */ |
e4d82761 | 477 | /*---------------------------- 0.25<|x|< 0.855469---------------------- */ |
6dbe713d SP |
478 | else if (k < 0x3feb6000) |
479 | { | |
25e440c6 | 480 | res = do_sin (x, 0, &cor); |
a87b5e95 | 481 | retval = (res == res + 1.096 * cor) ? res : slow1 (x); |
2bf49970 | 482 | retval = __copysign (retval, x); |
6dbe713d | 483 | } /* else if (k < 0x3feb6000) */ |
e4d82761 UD |
484 | |
485 | /*----------------------- 0.855469 <|x|<2.426265 ----------------------*/ | |
6dbe713d SP |
486 | else if (k < 0x400368fd) |
487 | { | |
488 | ||
9d84d0e5 | 489 | t = hp0 - fabs (x); |
758e79ec | 490 | res = do_cos (t, hp1, &cor); |
a87b5e95 | 491 | retval = (res == res + 1.020 * cor) ? res : slow2 (x); |
2bf49970 | 492 | retval = __copysign (retval, x); |
6dbe713d | 493 | } /* else if (k < 0x400368fd) */ |
e4d82761 | 494 | |
b3004556 | 495 | #ifndef IN_SINCOS |
e4d82761 | 496 | /*-------------------------- 2.426265<|x|< 105414350 ----------------------*/ |
6dbe713d SP |
497 | else if (k < 0x419921FB) |
498 | { | |
b3004556 SP |
499 | double a, da; |
500 | int4 n = reduce_sincos_1 (x, &a, &da); | |
b8b7e5e6 | 501 | retval = do_sincos_1 (a, da, x, n, false); |
6dbe713d | 502 | } /* else if (k < 0x419921FB ) */ |
e4d82761 UD |
503 | |
504 | /*---------------------105414350 <|x|< 281474976710656 --------------------*/ | |
6dbe713d SP |
505 | else if (k < 0x42F00000) |
506 | { | |
f7953c44 | 507 | double a, da; |
6dbe713d | 508 | |
f7953c44 | 509 | int4 n = reduce_sincos_2 (x, &a, &da); |
b8b7e5e6 | 510 | retval = do_sincos_2 (a, da, x, n, false); |
6dbe713d | 511 | } /* else if (k < 0x42F00000 ) */ |
e4d82761 UD |
512 | |
513 | /* -----------------281474976710656 <|x| <2^1024----------------------------*/ | |
6dbe713d | 514 | else if (k < 0x7ff00000) |
ead1ef37 | 515 | retval = reduce_and_compute (x, false); |
6dbe713d SP |
516 | |
517 | /*--------------------- |x| > 2^1024 ----------------------------------*/ | |
518 | else | |
519 | { | |
520 | if (k == 0x7ff00000 && u.i[LOW_HALF] == 0) | |
521 | __set_errno (EDOM); | |
522 | retval = x / x; | |
6dbe713d | 523 | } |
a045832d | 524 | #endif |
804360ed | 525 | |
6dbe713d | 526 | return retval; |
e4d82761 UD |
527 | } |
528 | ||
529 | ||
530 | /*******************************************************************/ | |
531 | /* An ultimate cos routine. Given an IEEE double machine number x */ | |
532 | /* it computes the correctly rounded (to nearest) value of cos(x) */ | |
533 | /*******************************************************************/ | |
534 | ||
463ac90d SP |
535 | #ifdef IN_SINCOS |
536 | static double | |
537 | #else | |
31d3cc00 UD |
538 | double |
539 | SECTION | |
463ac90d | 540 | #endif |
6dbe713d | 541 | __cos (double x) |
e4d82761 | 542 | { |
b3004556 SP |
543 | double y, xx, res, cor, a, da; |
544 | mynumber u; | |
545 | int4 k, m; | |
e4d82761 | 546 | |
804360ed JM |
547 | double retval = 0; |
548 | ||
463ac90d | 549 | #ifndef IN_SINCOS |
eb92c487 | 550 | SET_RESTORE_ROUND_53BIT (FE_TONEAREST); |
463ac90d | 551 | #endif |
804360ed | 552 | |
e4d82761 UD |
553 | u.x = x; |
554 | m = u.i[HIGH_HALF]; | |
6dbe713d | 555 | k = 0x7fffffff & m; |
e4d82761 | 556 | |
5eea0404 | 557 | /* |x|<2^-27 => cos(x)=1 */ |
6dbe713d | 558 | if (k < 0x3e400000) |
5eea0404 | 559 | retval = 1.0; |
6dbe713d SP |
560 | |
561 | else if (k < 0x3feb6000) | |
562 | { /* 2^-27 < |x| < 0.855469 */ | |
758e79ec | 563 | res = do_cos (x, 0, &cor); |
6dbe713d | 564 | retval = (res == res + 1.020 * cor) ? res : cslow2 (x); |
6dbe713d SP |
565 | } /* else if (k < 0x3feb6000) */ |
566 | ||
567 | else if (k < 0x400368fd) | |
568 | { /* 0.855469 <|x|<2.426265 */ ; | |
0e9be4db | 569 | y = hp0 - fabs (x); |
196f7f5d SP |
570 | a = y + hp1; |
571 | da = (y - a) + hp1; | |
6dbe713d SP |
572 | xx = a * a; |
573 | if (xx < 0.01588) | |
574 | { | |
8d561986 | 575 | res = TAYLOR_SIN (xx, a, da, cor); |
3459931a | 576 | cor = 1.02 * cor + __copysign (1.0e-31, cor); |
b8b7e5e6 | 577 | retval = (res == res + cor) ? res : sloww (a, da, x, true); |
6dbe713d SP |
578 | } |
579 | else | |
580 | { | |
758e79ec | 581 | res = do_sin (a, da, &cor); |
3459931a | 582 | cor = 1.035 * cor + __copysign (1.0e-31, cor); |
2bf49970 | 583 | retval = ((res == res + cor) ? __copysign (res, a) |
b8b7e5e6 | 584 | : sloww1 (a, da, x, true)); |
6dbe713d | 585 | } |
e4d82761 | 586 | |
6dbe713d | 587 | } /* else if (k < 0x400368fd) */ |
e4d82761 | 588 | |
e4d82761 | 589 | |
b3004556 | 590 | #ifndef IN_SINCOS |
6dbe713d SP |
591 | else if (k < 0x419921FB) |
592 | { /* 2.426265<|x|< 105414350 */ | |
b3004556 SP |
593 | double a, da; |
594 | int4 n = reduce_sincos_1 (x, &a, &da); | |
b8b7e5e6 | 595 | retval = do_sincos_1 (a, da, x, n, true); |
6dbe713d | 596 | } /* else if (k < 0x419921FB ) */ |
e4d82761 | 597 | |
6dbe713d SP |
598 | else if (k < 0x42F00000) |
599 | { | |
f7953c44 | 600 | double a, da; |
6dbe713d | 601 | |
f7953c44 | 602 | int4 n = reduce_sincos_2 (x, &a, &da); |
b8b7e5e6 | 603 | retval = do_sincos_2 (a, da, x, n, true); |
6dbe713d SP |
604 | } /* else if (k < 0x42F00000 ) */ |
605 | ||
6cce25f8 | 606 | /* 281474976710656 <|x| <2^1024 */ |
6dbe713d | 607 | else if (k < 0x7ff00000) |
ead1ef37 | 608 | retval = reduce_and_compute (x, true); |
e4d82761 | 609 | |
6dbe713d SP |
610 | else |
611 | { | |
612 | if (k == 0x7ff00000 && u.i[LOW_HALF] == 0) | |
613 | __set_errno (EDOM); | |
614 | retval = x / x; /* |x| > 2^1024 */ | |
e4d82761 | 615 | } |
a045832d | 616 | #endif |
e4d82761 | 617 | |
804360ed | 618 | return retval; |
e4d82761 UD |
619 | } |
620 | ||
621 | /************************************************************************/ | |
622 | /* Routine compute sin(x) for 2^-26 < |x|< 0.25 by Taylor with more */ | |
623 | /* precision and if still doesn't accurate enough by mpsin or dubsin */ | |
624 | /************************************************************************/ | |
625 | ||
54c86cca SP |
626 | static inline double |
627 | __always_inline | |
6dbe713d SP |
628 | slow (double x) |
629 | { | |
4aafb73c SP |
630 | double res, cor, w[2]; |
631 | res = TAYLOR_SLOW (x, 0, cor); | |
6dbe713d SP |
632 | if (res == res + 1.0007 * cor) |
633 | return res; | |
b7665e51 SP |
634 | |
635 | __dubsin (fabs (x), 0, w); | |
636 | if (w[0] == w[0] + 1.000000001 * w[1]) | |
2bf49970 | 637 | return __copysign (w[0], x); |
b7665e51 | 638 | |
2bf49970 | 639 | return __copysign (__mpsin (fabs (x), 0, false), x); |
e4d82761 | 640 | } |
6dbe713d | 641 | |
e4d82761 | 642 | /*******************************************************************************/ |
c5d5d574 | 643 | /* Routine compute sin(x) for 0.25<|x|< 0.855469 by __sincostab.tbl and Taylor */ |
e4d82761 UD |
644 | /* and if result still doesn't accurate enough by mpsin or dubsin */ |
645 | /*******************************************************************************/ | |
646 | ||
54c86cca SP |
647 | static inline double |
648 | __always_inline | |
6dbe713d SP |
649 | slow1 (double x) |
650 | { | |
758e79ec SP |
651 | double w[2], cor, res; |
652 | ||
653 | res = do_sin_slow (x, 0, 0, &cor); | |
392dd2de | 654 | if (res == res + cor) |
a87b5e95 | 655 | return res; |
b7665e51 SP |
656 | |
657 | __dubsin (fabs (x), 0, w); | |
658 | if (w[0] == w[0] + 1.000000005 * w[1]) | |
a87b5e95 | 659 | return w[0]; |
b7665e51 | 660 | |
a87b5e95 | 661 | return __mpsin (fabs (x), 0, false); |
e4d82761 | 662 | } |
6dbe713d | 663 | |
e4d82761 | 664 | /**************************************************************************/ |
af968f62 | 665 | /* Routine compute sin(x) for 0.855469 <|x|<2.426265 by __sincostab.tbl */ |
e4d82761 UD |
666 | /* and if result still doesn't accurate enough by mpsin or dubsin */ |
667 | /**************************************************************************/ | |
54c86cca SP |
668 | static inline double |
669 | __always_inline | |
6dbe713d SP |
670 | slow2 (double x) |
671 | { | |
758e79ec | 672 | double w[2], y, y1, y2, cor, res; |
7a74607f | 673 | |
9d84d0e5 | 674 | double t = hp0 - fabs (x); |
758e79ec | 675 | res = do_cos_slow (t, hp1, 0, &cor); |
392dd2de | 676 | if (res == res + cor) |
a87b5e95 | 677 | return res; |
b7665e51 SP |
678 | |
679 | y = fabs (x) - hp0; | |
680 | y1 = y - hp1; | |
681 | y2 = (y - y1) - hp1; | |
682 | __docos (y1, y2, w); | |
683 | if (w[0] == w[0] + 1.000000005 * w[1]) | |
a87b5e95 | 684 | return w[0]; |
b7665e51 | 685 | |
a87b5e95 | 686 | return __mpsin (fabs (x), 0, false); |
e4d82761 | 687 | } |
6dbe713d | 688 | |
36ee03e6 SP |
689 | /* Compute sin(x + dx) where X is small enough to use Taylor series around zero |
690 | and (x + dx) in the first or third quarter of the unit circle. ORIG is the | |
691 | original value of X for computing error of the result. If the result is not | |
692 | accurate enough, the routine calls mpsin or dubsin. SHIFT_QUADRANT rotates | |
693 | the unit circle by 1 to compute the cosine instead of sine. */ | |
54c86cca SP |
694 | static inline double |
695 | __always_inline | |
b8b7e5e6 | 696 | sloww (double x, double dx, double orig, bool shift_quadrant) |
6dbe713d | 697 | { |
4aafb73c | 698 | double y, t, res, cor, w[2], a, da, xn; |
8d561986 | 699 | mynumber v; |
e4d82761 | 700 | int4 n; |
4aafb73c | 701 | res = TAYLOR_SLOW (x, dx, cor); |
b3004556 | 702 | |
9d84d0e5 SP |
703 | double eps = fabs (orig) * 3.1e-30; |
704 | ||
3459931a | 705 | cor = 1.0005 * cor + __copysign (eps, cor); |
8d561986 | 706 | |
6dbe713d SP |
707 | if (res == res + cor) |
708 | return res; | |
b7665e51 | 709 | |
9d84d0e5 SP |
710 | a = fabs (x); |
711 | da = (x > 0) ? dx : -dx; | |
712 | __dubsin (a, da, w); | |
713 | eps = fabs (orig) * 1.1e-30; | |
3459931a | 714 | cor = 1.000000001 * w[1] + __copysign (eps, w[1]); |
b7665e51 SP |
715 | |
716 | if (w[0] == w[0] + cor) | |
2bf49970 | 717 | return __copysign (w[0], x); |
b7665e51 SP |
718 | |
719 | t = (orig * hpinv + toint); | |
720 | xn = t - toint; | |
721 | v.x = t; | |
722 | y = (orig - xn * mp1) - xn * mp2; | |
b8b7e5e6 | 723 | n = (v.i[LOW_HALF] + shift_quadrant) & 3; |
b7665e51 SP |
724 | da = xn * pp3; |
725 | t = y - da; | |
726 | da = (y - t) - da; | |
727 | y = xn * pp4; | |
728 | a = t - y; | |
729 | da = ((t - a) - y) + da; | |
b3004556 | 730 | |
cbf88869 | 731 | if (n & 2) |
6dbe713d | 732 | { |
b7665e51 SP |
733 | a = -a; |
734 | da = -da; | |
735 | } | |
9d84d0e5 SP |
736 | x = fabs (a); |
737 | dx = (a > 0) ? da : -da; | |
738 | __dubsin (x, dx, w); | |
739 | eps = fabs (orig) * 1.1e-40; | |
3459931a | 740 | cor = 1.000000001 * w[1] + __copysign (eps, w[1]); |
6dbe713d | 741 | |
b7665e51 | 742 | if (w[0] == w[0] + cor) |
2bf49970 | 743 | return __copysign (w[0], a); |
6dbe713d | 744 | |
b8b7e5e6 | 745 | return shift_quadrant ? __mpcos (orig, 0, true) : __mpsin (orig, 0, true); |
e4d82761 | 746 | } |
6dbe713d | 747 | |
36ee03e6 SP |
748 | /* Compute sin(x + dx) where X is in the first or third quarter of the unit |
749 | circle. ORIG is the original value of X for computing error of the result. | |
750 | If the result is not accurate enough, the routine calls mpsin or dubsin. | |
751 | SHIFT_QUADRANT rotates the unit circle by 1 to compute the cosine instead of | |
752 | sine. */ | |
54c86cca SP |
753 | static inline double |
754 | __always_inline | |
b8b7e5e6 | 755 | sloww1 (double x, double dx, double orig, bool shift_quadrant) |
6dbe713d | 756 | { |
758e79ec | 757 | double w[2], cor, res; |
6dbe713d | 758 | |
758e79ec | 759 | res = do_sin_slow (x, dx, 3.1e-30 * fabs (orig), &cor); |
6dbe713d SP |
760 | |
761 | if (res == res + cor) | |
2bf49970 | 762 | return __copysign (res, x); |
b7665e51 | 763 | |
758e79ec | 764 | dx = (x > 0 ? dx : -dx); |
9d84d0e5 | 765 | __dubsin (fabs (x), dx, w); |
b7665e51 | 766 | |
9d84d0e5 | 767 | double eps = 1.1e-30 * fabs (orig); |
3459931a | 768 | cor = 1.000000005 * w[1] + __copysign (eps, w[1]); |
6dbe713d | 769 | |
b7665e51 | 770 | if (w[0] == w[0] + cor) |
2bf49970 | 771 | return __copysign (w[0], x); |
6dbe713d | 772 | |
b8b7e5e6 | 773 | return shift_quadrant ? __mpcos (orig, 0, true) : __mpsin (orig, 0, true); |
e4d82761 | 774 | } |
6dbe713d | 775 | |
e4d82761 UD |
776 | /***************************************************************************/ |
777 | /* Routine compute sin(x+dx) (Double-Length number) where x in second or */ | |
778 | /* fourth quarter of unit circle.Routine receive also the original value */ | |
779 | /* and quarter(n= 1or 3)of x for computing error of result.And if result not*/ | |
780 | /* accurate enough routine calls mpsin1 or dubsin */ | |
781 | /***************************************************************************/ | |
782 | ||
54c86cca SP |
783 | static inline double |
784 | __always_inline | |
6dbe713d SP |
785 | sloww2 (double x, double dx, double orig, int n) |
786 | { | |
758e79ec | 787 | double w[2], cor, res; |
6dbe713d | 788 | |
758e79ec | 789 | res = do_cos_slow (x, dx, 3.1e-30 * fabs (orig), &cor); |
6dbe713d SP |
790 | |
791 | if (res == res + cor) | |
792 | return (n & 2) ? -res : res; | |
b7665e51 | 793 | |
758e79ec | 794 | dx = x > 0 ? dx : -dx; |
9d84d0e5 | 795 | __docos (fabs (x), dx, w); |
b7665e51 | 796 | |
9d84d0e5 | 797 | double eps = 1.1e-30 * fabs (orig); |
3459931a | 798 | cor = 1.000000005 * w[1] + __copysign (eps, w[1]); |
6dbe713d | 799 | |
b7665e51 SP |
800 | if (w[0] == w[0] + cor) |
801 | return (n & 2) ? -w[0] : w[0]; | |
6dbe713d | 802 | |
b3004556 | 803 | return (n & 1) ? __mpsin (orig, 0, true) : __mpcos (orig, 0, true); |
e4d82761 | 804 | } |
6dbe713d | 805 | |
e4d82761 UD |
806 | /***************************************************************************/ |
807 | /* Routine compute sin(x+dx) or cos(x+dx) (Double-Length number) where x */ | |
808 | /* is small enough to use Taylor series around zero and (x+dx) */ | |
809 | /* in first or third quarter of unit circle.Routine receive also */ | |
810 | /* (right argument) the original value of x for computing error of */ | |
811 | /* result.And if result not accurate enough routine calls other routines */ | |
812 | /***************************************************************************/ | |
813 | ||
54c86cca SP |
814 | static inline double |
815 | __always_inline | |
6dbe713d SP |
816 | bsloww (double x, double dx, double orig, int n) |
817 | { | |
9d84d0e5 | 818 | double res, cor, w[2], a, da; |
4aafb73c SP |
819 | |
820 | res = TAYLOR_SLOW (x, dx, cor); | |
3459931a | 821 | cor = 1.0005 * cor + __copysign (1.1e-24, cor); |
6dbe713d SP |
822 | if (res == res + cor) |
823 | return res; | |
b7665e51 | 824 | |
9d84d0e5 SP |
825 | a = fabs (x); |
826 | da = (x > 0) ? dx : -dx; | |
827 | __dubsin (a, da, w); | |
3459931a | 828 | cor = 1.000000001 * w[1] + __copysign (1.1e-24, w[1]); |
b7665e51 SP |
829 | |
830 | if (w[0] == w[0] + cor) | |
2bf49970 | 831 | return __copysign (w[0], x); |
b7665e51 SP |
832 | |
833 | return (n & 1) ? __mpcos (orig, 0, true) : __mpsin (orig, 0, true); | |
e4d82761 UD |
834 | } |
835 | ||
836 | /***************************************************************************/ | |
837 | /* Routine compute sin(x+dx) or cos(x+dx) (Double-Length number) where x */ | |
838 | /* in first or third quarter of unit circle.Routine receive also */ | |
839 | /* (right argument) the original value of x for computing error of result.*/ | |
840 | /* And if result not accurate enough routine calls other routines */ | |
841 | /***************************************************************************/ | |
842 | ||
54c86cca SP |
843 | static inline double |
844 | __always_inline | |
6dbe713d SP |
845 | bsloww1 (double x, double dx, double orig, int n) |
846 | { | |
758e79ec | 847 | double w[2], cor, res; |
6dbe713d | 848 | |
758e79ec | 849 | res = do_sin_slow (x, dx, 1.1e-24, &cor); |
6dbe713d SP |
850 | if (res == res + cor) |
851 | return (x > 0) ? res : -res; | |
b7665e51 | 852 | |
758e79ec | 853 | dx = (x > 0) ? dx : -dx; |
b7665e51 SP |
854 | __dubsin (fabs (x), dx, w); |
855 | ||
3459931a | 856 | cor = 1.000000005 * w[1] + __copysign (1.1e-24, w[1]); |
6dbe713d | 857 | |
b7665e51 | 858 | if (w[0] == w[0] + cor) |
2bf49970 | 859 | return __copysign (w[0], x); |
6dbe713d | 860 | |
b7665e51 | 861 | return (n & 1) ? __mpcos (orig, 0, true) : __mpsin (orig, 0, true); |
e4d82761 UD |
862 | } |
863 | ||
864 | /***************************************************************************/ | |
865 | /* Routine compute sin(x+dx) or cos(x+dx) (Double-Length number) where x */ | |
866 | /* in second or fourth quarter of unit circle.Routine receive also the */ | |
867 | /* original value and quarter(n= 1or 3)of x for computing error of result. */ | |
868 | /* And if result not accurate enough routine calls other routines */ | |
869 | /***************************************************************************/ | |
870 | ||
54c86cca SP |
871 | static inline double |
872 | __always_inline | |
6dbe713d SP |
873 | bsloww2 (double x, double dx, double orig, int n) |
874 | { | |
758e79ec | 875 | double w[2], cor, res; |
6dbe713d | 876 | |
758e79ec | 877 | res = do_cos_slow (x, dx, 1.1e-24, &cor); |
6dbe713d SP |
878 | if (res == res + cor) |
879 | return (n & 2) ? -res : res; | |
b7665e51 | 880 | |
758e79ec | 881 | dx = (x > 0) ? dx : -dx; |
b7665e51 SP |
882 | __docos (fabs (x), dx, w); |
883 | ||
3459931a | 884 | cor = 1.000000005 * w[1] + __copysign (1.1e-24, w[1]); |
6dbe713d | 885 | |
b7665e51 SP |
886 | if (w[0] == w[0] + cor) |
887 | return (n & 2) ? -w[0] : w[0]; | |
6dbe713d | 888 | |
b7665e51 | 889 | return (n & 1) ? __mpsin (orig, 0, true) : __mpcos (orig, 0, true); |
e4d82761 UD |
890 | } |
891 | ||
892 | /************************************************************************/ | |
893 | /* Routine compute cos(x) for 2^-27 < |x|< 0.25 by Taylor with more */ | |
894 | /* precision and if still doesn't accurate enough by mpcos or docos */ | |
895 | /************************************************************************/ | |
896 | ||
54c86cca SP |
897 | static inline double |
898 | __always_inline | |
6dbe713d SP |
899 | cslow2 (double x) |
900 | { | |
758e79ec | 901 | double w[2], cor, res; |
6dbe713d | 902 | |
758e79ec | 903 | res = do_cos_slow (x, 0, 0, &cor); |
392dd2de | 904 | if (res == res + cor) |
e4d82761 | 905 | return res; |
b7665e51 | 906 | |
758e79ec | 907 | __docos (fabs (x), 0, w); |
b7665e51 SP |
908 | if (w[0] == w[0] + 1.000000005 * w[1]) |
909 | return w[0]; | |
910 | ||
911 | return __mpcos (x, 0, false); | |
e4d82761 UD |
912 | } |
913 | ||
af968f62 | 914 | #ifndef __cos |
ca58f1db | 915 | weak_alias (__cos, cos) |
af968f62 UD |
916 | # ifdef NO_LONG_DOUBLE |
917 | strong_alias (__cos, __cosl) | |
918 | weak_alias (__cos, cosl) | |
919 | # endif | |
920 | #endif | |
921 | #ifndef __sin | |
ca58f1db | 922 | weak_alias (__sin, sin) |
af968f62 | 923 | # ifdef NO_LONG_DOUBLE |
ca58f1db UD |
924 | strong_alias (__sin, __sinl) |
925 | weak_alias (__sin, sinl) | |
af968f62 | 926 | # endif |
cccda09f | 927 | #endif |