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
114extern const union
115{
116 int4 i[880];
117 double x[440];
118} __sincostab attribute_hidden;
119
120static 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
127static const double t22 = 0x1.8p22;
128
129void __dubsin (double x, double dx, double w[]);
130void __docos (double x, double dx, double w[]);
131double __mpsin (double x, double dx, bool reduce_range);
132double __mpcos (double x, double dx, bool reduce_range);
133static double slow (double x);
134static double slow1 (double x);
135static double slow2 (double x);
136static double sloww (double x, double dx, double orig, bool shift_quadrant);
137static double sloww1 (double x, double dx, double orig, bool shift_quadrant);
138static double sloww2 (double x, double dx, double orig, int n);
139static double bsloww (double x, double dx, double orig, int n);
140static double bsloww1 (double x, double dx, double orig, int n);
141static double bsloww2 (double x, double dx, double orig, int n);
142int __branred (double x, double *a, double *aa);
143static 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. */
149static inline double
150__always_inline
151do_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. */
175static inline double
176__always_inline
177do_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. */
211static inline double
212__always_inline
213do_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. */
236static inline double
237__always_inline
238do_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. */
270static inline double
271__always_inline
272reduce_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
298static inline int4
299__always_inline
300reduce_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. */
321static double
322__always_inline
323do_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
365static inline int4
366__always_inline
367reduce_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. */
393static double
394__always_inline
395do_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
444static double
445#else
446double
447SECTION
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
537static double
538#else
539double
540SECTION
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
627static inline double
628__always_inline
629slow (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
648static inline double
649__always_inline
650slow1 (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/**************************************************************************/
669static inline double
670__always_inline
671slow2 (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. */
695static inline double
696__always_inline
697sloww (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. */
754static inline double
755__always_inline
756sloww1 (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
784static inline double
785__always_inline
786sloww2 (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
815static inline double
816__always_inline
817bsloww (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
844static inline double
845__always_inline
846bsloww1 (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
872static inline double
873__always_inline
874bsloww2 (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
898static inline double
899__always_inline
900cslow2 (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
916libm_alias_double (__cos, cos)
917#endif
918#ifndef __sin
919libm_alias_double (__sin, sin)
920#endif
921