1 | /* |
2 | * IBM Accurate Mathematical Library |
3 | * written by International Business Machines Corp. |
4 | * Copyright (C) 2001-2017 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 <fenv.h> |
56 | |
57 | /* Helper macros to compute sin of the input values. */ |
58 | #define POLYNOMIAL2(xx) ((((s5 * (xx) + s4) * (xx) + s3) * (xx) + s2) * (xx)) |
59 | |
60 | #define POLYNOMIAL(xx) (POLYNOMIAL2 (xx) + s1) |
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. */ |
69 | #define TAYLOR_SIN(xx, a, da, cor) \ |
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; \ |
89 | double y = aa * x1 * x1 * x1; \ |
90 | double r = (x0) + y; \ |
91 | double x2 = ((x0) - x1) + (dx); \ |
92 | double t = (((POLYNOMIAL2 (xx) + bb) * xx + 3.0 * aa * x1 * x2) \ |
93 | * (x0) + aa * x2 * x2 * x2 + (dx)); \ |
94 | t = (((x0) - r) + y) + t; \ |
95 | double res = r + t; \ |
96 | (cor) = (r - res) + t; \ |
97 | res; \ |
98 | }) |
99 | |
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 | |
109 | #ifndef SECTION |
110 | # define SECTION |
111 | #endif |
112 | |
113 | extern const union |
114 | { |
115 | int4 i[880]; |
116 | double x[440]; |
117 | } __sincostab attribute_hidden; |
118 | |
119 | static const double |
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 | |
126 | static const double t22 = 0x1.8p22; |
127 | |
128 | void __dubsin (double x, double dx, double w[]); |
129 | void __docos (double x, double dx, double w[]); |
130 | double __mpsin (double x, double dx, bool reduce_range); |
131 | double __mpcos (double x, double dx, bool reduce_range); |
132 | static double slow (double x); |
133 | static double slow1 (double x); |
134 | static double slow2 (double x); |
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); |
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); |
143 | |
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. */ |
148 | static inline double |
149 | __always_inline |
150 | do_cos (double x, double dx, double *corp) |
151 | { |
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 | |
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 | |
172 | /* A more precise variant of DO_COS. EPS is the adjustment to the correction |
173 | COR. */ |
174 | static inline double |
175 | __always_inline |
176 | do_cos_slow (double x, double dx, double eps, double *corp) |
177 | { |
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 | |
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; |
201 | cor = 1.0005 * cor + __copysign (eps, cor); |
202 | *corp = cor; |
203 | return res; |
204 | } |
205 | |
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. */ |
210 | static inline double |
211 | __always_inline |
212 | do_sin (double x, double dx, double *corp) |
213 | { |
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 | |
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 | |
233 | /* A more precise variant of DO_SIN. EPS is the adjustment to the correction |
234 | COR. */ |
235 | static inline double |
236 | __always_inline |
237 | do_sin_slow (double x, double dx, double eps, double *corp) |
238 | { |
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 | |
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; |
261 | cor = 1.0005 * cor + __copysign (eps, cor); |
262 | *corp = cor; |
263 | return res; |
264 | } |
265 | |
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. */ |
269 | static inline double |
270 | __always_inline |
271 | reduce_and_compute (double x, bool shift_quadrant) |
272 | { |
273 | double retval = 0, a, da; |
274 | unsigned int n = __branred (x, &a, &da); |
275 | int4 k = (n + shift_quadrant) % 4; |
276 | switch (k) |
277 | { |
278 | case 2: |
279 | a = -a; |
280 | da = -da; |
281 | /* Fall through. */ |
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; |
288 | |
289 | case 1: |
290 | case 3: |
291 | retval = bsloww2 (a, da, x, n); |
292 | break; |
293 | } |
294 | return retval; |
295 | } |
296 | |
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 | |
318 | /* Compute sin (A + DA). cos can be computed by passing SHIFT_QUADRANT as |
319 | true, which results in shifting the quadrant N clockwise. */ |
320 | static double |
321 | __always_inline |
322 | do_sincos_1 (double a, double da, double x, int4 n, bool shift_quadrant) |
323 | { |
324 | double xx, retval, res, cor; |
325 | double eps = fabs (x) * 1.2e-30; |
326 | |
327 | int k1 = (n + shift_quadrant) & 3; |
328 | switch (k1) |
329 | { /* quarter of unit circle */ |
330 | case 2: |
331 | a = -a; |
332 | da = -da; |
333 | /* Fall through. */ |
334 | case 0: |
335 | xx = a * a; |
336 | if (xx < 0.01588) |
337 | { |
338 | /* Taylor series. */ |
339 | res = TAYLOR_SIN (xx, a, da, cor); |
340 | cor = 1.02 * cor + __copysign (eps, cor); |
341 | retval = (res == res + cor) ? res : sloww (a, da, x, shift_quadrant); |
342 | } |
343 | else |
344 | { |
345 | res = do_sin (a, da, &cor); |
346 | cor = 1.035 * cor + __copysign (eps, cor); |
347 | retval = ((res == res + cor) ? __copysign (res, a) |
348 | : sloww1 (a, da, x, shift_quadrant)); |
349 | } |
350 | break; |
351 | |
352 | case 1: |
353 | case 3: |
354 | res = do_cos (a, da, &cor); |
355 | cor = 1.025 * cor + __copysign (eps, cor); |
356 | retval = ((res == res + cor) ? ((n & 2) ? -res : res) |
357 | : sloww2 (a, da, x, n)); |
358 | break; |
359 | } |
360 | |
361 | return retval; |
362 | } |
363 | |
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 | |
390 | /* Compute sin (A + DA). cos can be computed by passing SHIFT_QUADRANT as |
391 | true, which results in shifting the quadrant N clockwise. */ |
392 | static double |
393 | __always_inline |
394 | do_sincos_2 (double a, double da, double x, int4 n, bool shift_quadrant) |
395 | { |
396 | double res, retval, cor, xx; |
397 | |
398 | double eps = 1.0e-24; |
399 | |
400 | int4 k = (n + shift_quadrant) & 3; |
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); |
414 | cor = 1.02 * cor + __copysign (eps, cor); |
415 | retval = (res == res + cor) ? res : bsloww (a, da, x, n); |
416 | } |
417 | else |
418 | { |
419 | res = do_sin (a, da, &cor); |
420 | cor = 1.035 * cor + __copysign (eps, cor); |
421 | retval = ((res == res + cor) ? __copysign (res, a) |
422 | : bsloww1 (a, da, x, n)); |
423 | } |
424 | break; |
425 | |
426 | case 1: |
427 | case 3: |
428 | res = do_cos (a, da, &cor); |
429 | cor = 1.025 * cor + __copysign (eps, cor); |
430 | retval = ((res == res + cor) ? ((n & 2) ? -res : res) |
431 | : bsloww2 (a, da, x, n)); |
432 | break; |
433 | } |
434 | |
435 | return retval; |
436 | } |
437 | |
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 | /*******************************************************************/ |
442 | #ifdef IN_SINCOS |
443 | static double |
444 | #else |
445 | double |
446 | SECTION |
447 | #endif |
448 | __sin (double x) |
449 | { |
450 | double xx, res, t, cor; |
451 | mynumber u; |
452 | int4 k, m; |
453 | double retval = 0; |
454 | |
455 | #ifndef IN_SINCOS |
456 | SET_RESTORE_ROUND_53BIT (FE_TONEAREST); |
457 | #endif |
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 */ |
463 | { |
464 | math_check_force_underflow (x); |
465 | retval = x; |
466 | } |
467 | /*---------------------------- 2^-26 < |x|< 0.25 ----------------------*/ |
468 | else if (k < 0x3fd00000) |
469 | { |
470 | xx = x * x; |
471 | /* Taylor series. */ |
472 | t = POLYNOMIAL (xx) * (xx * x); |
473 | res = x + t; |
474 | cor = (x - res) + t; |
475 | retval = (res == res + 1.07 * cor) ? res : slow (x); |
476 | } /* else if (k < 0x3fd00000) */ |
477 | /*---------------------------- 0.25<|x|< 0.855469---------------------- */ |
478 | else if (k < 0x3feb6000) |
479 | { |
480 | res = do_sin (x, 0, &cor); |
481 | retval = (res == res + 1.096 * cor) ? res : slow1 (x); |
482 | retval = __copysign (retval, x); |
483 | } /* else if (k < 0x3feb6000) */ |
484 | |
485 | /*----------------------- 0.855469 <|x|<2.426265 ----------------------*/ |
486 | else if (k < 0x400368fd) |
487 | { |
488 | |
489 | t = hp0 - fabs (x); |
490 | res = do_cos (t, hp1, &cor); |
491 | retval = (res == res + 1.020 * cor) ? res : slow2 (x); |
492 | retval = __copysign (retval, x); |
493 | } /* else if (k < 0x400368fd) */ |
494 | |
495 | #ifndef IN_SINCOS |
496 | /*-------------------------- 2.426265<|x|< 105414350 ----------------------*/ |
497 | else if (k < 0x419921FB) |
498 | { |
499 | double a, da; |
500 | int4 n = reduce_sincos_1 (x, &a, &da); |
501 | retval = do_sincos_1 (a, da, x, n, false); |
502 | } /* else if (k < 0x419921FB ) */ |
503 | |
504 | /*---------------------105414350 <|x|< 281474976710656 --------------------*/ |
505 | else if (k < 0x42F00000) |
506 | { |
507 | double a, da; |
508 | |
509 | int4 n = reduce_sincos_2 (x, &a, &da); |
510 | retval = do_sincos_2 (a, da, x, n, false); |
511 | } /* else if (k < 0x42F00000 ) */ |
512 | |
513 | /* -----------------281474976710656 <|x| <2^1024----------------------------*/ |
514 | else if (k < 0x7ff00000) |
515 | retval = reduce_and_compute (x, false); |
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; |
523 | } |
524 | #endif |
525 | |
526 | return retval; |
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 | |
535 | #ifdef IN_SINCOS |
536 | static double |
537 | #else |
538 | double |
539 | SECTION |
540 | #endif |
541 | __cos (double x) |
542 | { |
543 | double y, xx, res, cor, a, da; |
544 | mynumber u; |
545 | int4 k, m; |
546 | |
547 | double retval = 0; |
548 | |
549 | #ifndef IN_SINCOS |
550 | SET_RESTORE_ROUND_53BIT (FE_TONEAREST); |
551 | #endif |
552 | |
553 | u.x = x; |
554 | m = u.i[HIGH_HALF]; |
555 | k = 0x7fffffff & m; |
556 | |
557 | /* |x|<2^-27 => cos(x)=1 */ |
558 | if (k < 0x3e400000) |
559 | retval = 1.0; |
560 | |
561 | else if (k < 0x3feb6000) |
562 | { /* 2^-27 < |x| < 0.855469 */ |
563 | res = do_cos (x, 0, &cor); |
564 | retval = (res == res + 1.020 * cor) ? res : cslow2 (x); |
565 | } /* else if (k < 0x3feb6000) */ |
566 | |
567 | else if (k < 0x400368fd) |
568 | { /* 0.855469 <|x|<2.426265 */ ; |
569 | y = hp0 - fabs (x); |
570 | a = y + hp1; |
571 | da = (y - a) + hp1; |
572 | xx = a * a; |
573 | if (xx < 0.01588) |
574 | { |
575 | res = TAYLOR_SIN (xx, a, da, cor); |
576 | cor = 1.02 * cor + __copysign (1.0e-31, cor); |
577 | retval = (res == res + cor) ? res : sloww (a, da, x, true); |
578 | } |
579 | else |
580 | { |
581 | res = do_sin (a, da, &cor); |
582 | cor = 1.035 * cor + __copysign (1.0e-31, cor); |
583 | retval = ((res == res + cor) ? __copysign (res, a) |
584 | : sloww1 (a, da, x, true)); |
585 | } |
586 | |
587 | } /* else if (k < 0x400368fd) */ |
588 | |
589 | |
590 | #ifndef IN_SINCOS |
591 | else if (k < 0x419921FB) |
592 | { /* 2.426265<|x|< 105414350 */ |
593 | double a, da; |
594 | int4 n = reduce_sincos_1 (x, &a, &da); |
595 | retval = do_sincos_1 (a, da, x, n, true); |
596 | } /* else if (k < 0x419921FB ) */ |
597 | |
598 | else if (k < 0x42F00000) |
599 | { |
600 | double a, da; |
601 | |
602 | int4 n = reduce_sincos_2 (x, &a, &da); |
603 | retval = do_sincos_2 (a, da, x, n, true); |
604 | } /* else if (k < 0x42F00000 ) */ |
605 | |
606 | /* 281474976710656 <|x| <2^1024 */ |
607 | else if (k < 0x7ff00000) |
608 | retval = reduce_and_compute (x, true); |
609 | |
610 | else |
611 | { |
612 | if (k == 0x7ff00000 && u.i[LOW_HALF] == 0) |
613 | __set_errno (EDOM); |
614 | retval = x / x; /* |x| > 2^1024 */ |
615 | } |
616 | #endif |
617 | |
618 | return retval; |
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 | |
626 | static inline double |
627 | __always_inline |
628 | slow (double x) |
629 | { |
630 | double res, cor, w[2]; |
631 | res = TAYLOR_SLOW (x, 0, cor); |
632 | if (res == res + 1.0007 * cor) |
633 | return res; |
634 | |
635 | __dubsin (fabs (x), 0, w); |
636 | if (w[0] == w[0] + 1.000000001 * w[1]) |
637 | return __copysign (w[0], x); |
638 | |
639 | return __copysign (__mpsin (fabs (x), 0, false), x); |
640 | } |
641 | |
642 | /*******************************************************************************/ |
643 | /* Routine compute sin(x) for 0.25<|x|< 0.855469 by __sincostab.tbl and Taylor */ |
644 | /* and if result still doesn't accurate enough by mpsin or dubsin */ |
645 | /*******************************************************************************/ |
646 | |
647 | static inline double |
648 | __always_inline |
649 | slow1 (double x) |
650 | { |
651 | double w[2], cor, res; |
652 | |
653 | res = do_sin_slow (x, 0, 0, &cor); |
654 | if (res == res + cor) |
655 | return res; |
656 | |
657 | __dubsin (fabs (x), 0, w); |
658 | if (w[0] == w[0] + 1.000000005 * w[1]) |
659 | return w[0]; |
660 | |
661 | return __mpsin (fabs (x), 0, false); |
662 | } |
663 | |
664 | /**************************************************************************/ |
665 | /* Routine compute sin(x) for 0.855469 <|x|<2.426265 by __sincostab.tbl */ |
666 | /* and if result still doesn't accurate enough by mpsin or dubsin */ |
667 | /**************************************************************************/ |
668 | static inline double |
669 | __always_inline |
670 | slow2 (double x) |
671 | { |
672 | double w[2], y, y1, y2, cor, res; |
673 | |
674 | double t = hp0 - fabs (x); |
675 | res = do_cos_slow (t, hp1, 0, &cor); |
676 | if (res == res + cor) |
677 | return res; |
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]) |
684 | return w[0]; |
685 | |
686 | return __mpsin (fabs (x), 0, false); |
687 | } |
688 | |
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. */ |
694 | static inline double |
695 | __always_inline |
696 | sloww (double x, double dx, double orig, bool shift_quadrant) |
697 | { |
698 | double y, t, res, cor, w[2], a, da, xn; |
699 | mynumber v; |
700 | int4 n; |
701 | res = TAYLOR_SLOW (x, dx, cor); |
702 | |
703 | double eps = fabs (orig) * 3.1e-30; |
704 | |
705 | cor = 1.0005 * cor + __copysign (eps, cor); |
706 | |
707 | if (res == res + cor) |
708 | return res; |
709 | |
710 | a = fabs (x); |
711 | da = (x > 0) ? dx : -dx; |
712 | __dubsin (a, da, w); |
713 | eps = fabs (orig) * 1.1e-30; |
714 | cor = 1.000000001 * w[1] + __copysign (eps, w[1]); |
715 | |
716 | if (w[0] == w[0] + cor) |
717 | return __copysign (w[0], x); |
718 | |
719 | t = (orig * hpinv + toint); |
720 | xn = t - toint; |
721 | v.x = t; |
722 | y = (orig - xn * mp1) - xn * mp2; |
723 | n = (v.i[LOW_HALF] + shift_quadrant) & 3; |
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; |
730 | |
731 | if (n & 2) |
732 | { |
733 | a = -a; |
734 | da = -da; |
735 | } |
736 | x = fabs (a); |
737 | dx = (a > 0) ? da : -da; |
738 | __dubsin (x, dx, w); |
739 | eps = fabs (orig) * 1.1e-40; |
740 | cor = 1.000000001 * w[1] + __copysign (eps, w[1]); |
741 | |
742 | if (w[0] == w[0] + cor) |
743 | return __copysign (w[0], a); |
744 | |
745 | return shift_quadrant ? __mpcos (orig, 0, true) : __mpsin (orig, 0, true); |
746 | } |
747 | |
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. */ |
753 | static inline double |
754 | __always_inline |
755 | sloww1 (double x, double dx, double orig, bool shift_quadrant) |
756 | { |
757 | double w[2], cor, res; |
758 | |
759 | res = do_sin_slow (x, dx, 3.1e-30 * fabs (orig), &cor); |
760 | |
761 | if (res == res + cor) |
762 | return __copysign (res, x); |
763 | |
764 | dx = (x > 0 ? dx : -dx); |
765 | __dubsin (fabs (x), dx, w); |
766 | |
767 | double eps = 1.1e-30 * fabs (orig); |
768 | cor = 1.000000005 * w[1] + __copysign (eps, w[1]); |
769 | |
770 | if (w[0] == w[0] + cor) |
771 | return __copysign (w[0], x); |
772 | |
773 | return shift_quadrant ? __mpcos (orig, 0, true) : __mpsin (orig, 0, true); |
774 | } |
775 | |
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 | |
783 | static inline double |
784 | __always_inline |
785 | sloww2 (double x, double dx, double orig, int n) |
786 | { |
787 | double w[2], cor, res; |
788 | |
789 | res = do_cos_slow (x, dx, 3.1e-30 * fabs (orig), &cor); |
790 | |
791 | if (res == res + cor) |
792 | return (n & 2) ? -res : res; |
793 | |
794 | dx = x > 0 ? dx : -dx; |
795 | __docos (fabs (x), dx, w); |
796 | |
797 | double eps = 1.1e-30 * fabs (orig); |
798 | cor = 1.000000005 * w[1] + __copysign (eps, w[1]); |
799 | |
800 | if (w[0] == w[0] + cor) |
801 | return (n & 2) ? -w[0] : w[0]; |
802 | |
803 | return (n & 1) ? __mpsin (orig, 0, true) : __mpcos (orig, 0, true); |
804 | } |
805 | |
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 | |
814 | static inline double |
815 | __always_inline |
816 | bsloww (double x, double dx, double orig, int n) |
817 | { |
818 | double res, cor, w[2], a, da; |
819 | |
820 | res = TAYLOR_SLOW (x, dx, cor); |
821 | cor = 1.0005 * cor + __copysign (1.1e-24, cor); |
822 | if (res == res + cor) |
823 | return res; |
824 | |
825 | a = fabs (x); |
826 | da = (x > 0) ? dx : -dx; |
827 | __dubsin (a, da, w); |
828 | cor = 1.000000001 * w[1] + __copysign (1.1e-24, w[1]); |
829 | |
830 | if (w[0] == w[0] + cor) |
831 | return __copysign (w[0], x); |
832 | |
833 | return (n & 1) ? __mpcos (orig, 0, true) : __mpsin (orig, 0, true); |
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 | |
843 | static inline double |
844 | __always_inline |
845 | bsloww1 (double x, double dx, double orig, int n) |
846 | { |
847 | double w[2], cor, res; |
848 | |
849 | res = do_sin_slow (x, dx, 1.1e-24, &cor); |
850 | if (res == res + cor) |
851 | return (x > 0) ? res : -res; |
852 | |
853 | dx = (x > 0) ? dx : -dx; |
854 | __dubsin (fabs (x), dx, w); |
855 | |
856 | cor = 1.000000005 * w[1] + __copysign (1.1e-24, w[1]); |
857 | |
858 | if (w[0] == w[0] + cor) |
859 | return __copysign (w[0], x); |
860 | |
861 | return (n & 1) ? __mpcos (orig, 0, true) : __mpsin (orig, 0, true); |
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 | |
871 | static inline double |
872 | __always_inline |
873 | bsloww2 (double x, double dx, double orig, int n) |
874 | { |
875 | double w[2], cor, res; |
876 | |
877 | res = do_cos_slow (x, dx, 1.1e-24, &cor); |
878 | if (res == res + cor) |
879 | return (n & 2) ? -res : res; |
880 | |
881 | dx = (x > 0) ? dx : -dx; |
882 | __docos (fabs (x), dx, w); |
883 | |
884 | cor = 1.000000005 * w[1] + __copysign (1.1e-24, w[1]); |
885 | |
886 | if (w[0] == w[0] + cor) |
887 | return (n & 2) ? -w[0] : w[0]; |
888 | |
889 | return (n & 1) ? __mpsin (orig, 0, true) : __mpcos (orig, 0, true); |
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 | |
897 | static inline double |
898 | __always_inline |
899 | cslow2 (double x) |
900 | { |
901 | double w[2], cor, res; |
902 | |
903 | res = do_cos_slow (x, 0, 0, &cor); |
904 | if (res == res + cor) |
905 | return res; |
906 | |
907 | __docos (fabs (x), 0, w); |
908 | if (w[0] == w[0] + 1.000000005 * w[1]) |
909 | return w[0]; |
910 | |
911 | return __mpcos (x, 0, false); |
912 | } |
913 | |
914 | #ifndef __cos |
915 | weak_alias (__cos, cos) |
916 | # ifdef NO_LONG_DOUBLE |
917 | strong_alias (__cos, __cosl) |
918 | weak_alias (__cos, cosl) |
919 | # endif |
920 | #endif |
921 | #ifndef __sin |
922 | weak_alias (__sin, sin) |
923 | # ifdef NO_LONG_DOUBLE |
924 | strong_alias (__sin, __sinl) |
925 | weak_alias (__sin, sinl) |
926 | # endif |
927 | #endif |
928 | |