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