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
113extern const union
114{
115 int4 i[880];
116 double x[440];
117} __sincostab attribute_hidden;
118
119static 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
126static const double t22 = 0x1.8p22;
127
128void __dubsin (double x, double dx, double w[]);
129void __docos (double x, double dx, double w[]);
130double __mpsin (double x, double dx, bool reduce_range);
131double __mpcos (double x, double dx, bool reduce_range);
132static double slow (double x);
133static double slow1 (double x);
134static double slow2 (double x);
135static double sloww (double x, double dx, double orig, bool shift_quadrant);
136static double sloww1 (double x, double dx, double orig, bool shift_quadrant);
137static double sloww2 (double x, double dx, double orig, int n);
138static double bsloww (double x, double dx, double orig, int n);
139static double bsloww1 (double x, double dx, double orig, int n);
140static double bsloww2 (double x, double dx, double orig, int n);
141int __branred (double x, double *a, double *aa);
142static 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. */
148static inline double
149__always_inline
150do_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. */
174static inline double
175__always_inline
176do_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. */
210static inline double
211__always_inline
212do_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. */
235static inline double
236__always_inline
237do_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. */
269static inline double
270__always_inline
271reduce_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
297static inline int4
298__always_inline
299reduce_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. */
320static double
321__always_inline
322do_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
364static inline int4
365__always_inline
366reduce_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. */
392static double
393__always_inline
394do_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
443static double
444#else
445double
446SECTION
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
536static double
537#else
538double
539SECTION
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
626static inline double
627__always_inline
628slow (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
647static inline double
648__always_inline
649slow1 (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/**************************************************************************/
668static inline double
669__always_inline
670slow2 (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. */
694static inline double
695__always_inline
696sloww (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. */
753static inline double
754__always_inline
755sloww1 (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
783static inline double
784__always_inline
785sloww2 (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
814static inline double
815__always_inline
816bsloww (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
843static inline double
844__always_inline
845bsloww1 (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
871static inline double
872__always_inline
873bsloww2 (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
897static inline double
898__always_inline
899cslow2 (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
915weak_alias (__cos, cos)
916# ifdef NO_LONG_DOUBLE
917strong_alias (__cos, __cosl)
918weak_alias (__cos, cosl)
919# endif
920#endif
921#ifndef __sin
922weak_alias (__sin, sin)
923# ifdef NO_LONG_DOUBLE
924strong_alias (__sin, __sinl)
925weak_alias (__sin, sinl)
926# endif
927#endif
928