2 * IBM Accurate Mathematical Library
3 * written by International Business Machines Corp.
4 * Copyright (C) 2001-2015 Free Software Foundation, Inc.
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.
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.
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/>.
19 /****************************************************************************/
21 /* MODULE_NAME:usncs.c */
35 /* FILES NEEDED: dla.h endian.h mpa.h mydefs.h usncs.h */
36 /* branred.c sincos32.c dosincos.c mpa.c */
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. */
44 /****************************************************************************/
54 #include <math_private.h>
57 /* Helper macros to compute sin of the input values. */
58 #define POLYNOMIAL2(xx) ((((s5 * (xx) + s4) * (xx) + s3) * (xx) + s2) * (xx))
60 #define POLYNOMIAL(xx) (POLYNOMIAL2 (xx) + s1)
62 /* The computed polynomial is a variation of the Taylor series expansion for
65 a - a^3/3! + a^5/5! - a^7/7! + a^9/9! + (1 - a^2) * da / 2
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) \
71 double t = ((POLYNOMIAL (xx) * (a) - 0.5 * (da)) * (xx) + (da)); \
72 double res = (a) + t; \
73 (cor) = ((a) - res) + t; \
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:
80 bb * x ^ 3 + 3 * aa * x * x1 * x2 + aa * x1 ^ 3 + aa * x2 ^ 3
82 The correction term is dx and bb + aa = -1/3!
84 #define TAYLOR_SLOW(x0, dx, cor) \
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; \
96 (cor) = (r - res) + t; \
100 #define SINCOS_TABLE_LOOKUP(u, sn, ssn, cs, ccs) \
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]; \
117 } __sincostab attribute_hidden
;
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;
126 static const double t22
= 0x1.8p22
;
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
, int n
);
136 static double sloww1 (double x
, double dx
, double orig
, int m
, int n
);
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
);
144 /* Given a number partitioned into U and X such that U is an index into the
145 sin/cos table, this macro computes the cosine of the number by combining
146 the sin and cos of X (as computed by a variation of the Taylor series) with
147 the values looked up from the sin/cos table to get the result in RES and a
148 correction value in COR. */
150 do_cos (mynumber u
, double x
, double *corp
)
152 double xx
, s
, sn
, ssn
, c
, cs
, ccs
, res
, cor
;
154 s
= x
+ x
* xx
* (sn3
+ xx
* sn5
);
155 c
= xx
* (cs2
+ xx
* (cs4
+ xx
* cs6
));
156 SINCOS_TABLE_LOOKUP (u
, sn
, ssn
, cs
, ccs
);
157 cor
= (ccs
- s
* ssn
- cs
* c
) - sn
* s
;
159 cor
= (cs
- res
) + cor
;
164 /* A more precise variant of DO_COS where the number is partitioned into U, X
165 and DX. EPS is the adjustment to the correction COR. */
167 do_cos_slow (mynumber u
, double x
, double dx
, double eps
, double *corp
)
169 double xx
, y
, x1
, x2
, e1
, e2
, res
, cor
;
170 double s
, sn
, ssn
, c
, cs
, ccs
;
172 s
= x
* xx
* (sn3
+ xx
* sn5
);
173 c
= x
* dx
+ xx
* (cs2
+ xx
* (cs4
+ xx
* cs6
));
174 SINCOS_TABLE_LOOKUP (u
, sn
, ssn
, cs
, ccs
);
175 x1
= (x
+ t22
) - t22
;
177 e1
= (sn
+ t22
) - t22
;
178 e2
= (sn
- e1
) + ssn
;
179 cor
= (ccs
- cs
* c
- e1
* x2
- e2
* x
) - sn
* s
;
181 cor
= cor
+ ((cs
- y
) - e1
* x1
);
183 cor
= (y
- res
) + cor
;
185 cor
= 1.0005 * cor
+ eps
;
187 cor
= 1.0005 * cor
- eps
;
192 /* Given a number partitioned into U and X and DX such that U is an index into
193 the sin/cos table, this macro computes the sine of the number by combining
194 the sin and cos of X (as computed by a variation of the Taylor series) with
195 the values looked up from the sin/cos table to get the result in RES and a
196 correction value in COR. */
198 do_sin (mynumber u
, double x
, double dx
, double *corp
)
200 double xx
, s
, sn
, ssn
, c
, cs
, ccs
, cor
, res
;
202 s
= x
+ (dx
+ x
* xx
* (sn3
+ xx
* sn5
));
203 c
= x
* dx
+ xx
* (cs2
+ xx
* (cs4
+ xx
* cs6
));
204 SINCOS_TABLE_LOOKUP (u
, sn
, ssn
, cs
, ccs
);
205 cor
= (ssn
+ s
* ccs
- sn
* c
) + cs
* s
;
207 cor
= (sn
- res
) + cor
;
212 /* A more precise variant of res = do_sin where the number is partitioned into U, X
213 and DX. EPS is the adjustment to the correction COR. */
215 do_sin_slow (mynumber u
, double x
, double dx
, double eps
, double *corp
)
217 double xx
, y
, x1
, x2
, c1
, c2
, res
, cor
;
218 double s
, sn
, ssn
, c
, cs
, ccs
;
220 s
= x
* xx
* (sn3
+ xx
* sn5
);
221 c
= xx
* (cs2
+ xx
* (cs4
+ xx
* cs6
));
222 SINCOS_TABLE_LOOKUP (u
, sn
, ssn
, cs
, ccs
);
223 x1
= (x
+ t22
) - t22
;
225 c1
= (cs
+ t22
) - t22
;
226 c2
= (cs
- c1
) + ccs
;
227 cor
= (ssn
+ s
* ccs
+ cs
* s
+ c2
* x
+ c1
* x2
- sn
* x
* dx
) - sn
* c
;
229 cor
= cor
+ ((sn
- y
) + c1
* x1
);
231 cor
= (y
- res
) + cor
;
233 cor
= 1.0005 * cor
+ eps
;
235 cor
= 1.0005 * cor
- eps
;
240 /* Reduce range of X and compute sin of a + da. K is the amount by which to
241 rotate the quadrants. This allows us to use the same routine to compute cos
242 by simply rotating the quadrants by 1. */
245 reduce_and_compute (double x
, unsigned int k
)
247 double retval
= 0, a
, da
;
248 unsigned int n
= __branred (x
, &a
, &da
);
254 retval
= bsloww (a
, da
, x
, n
);
256 retval
= bsloww1 (a
, da
, x
, n
);
260 retval
= bsloww (-a
, -da
, x
, n
);
262 retval
= bsloww1 (-a
, -da
, x
, n
);
267 retval
= bsloww2 (a
, da
, x
, n
);
275 reduce_sincos_1 (double x
, double *a
, double *da
)
279 double t
= (x
* hpinv
+ toint
);
280 double xn
= t
- toint
;
282 double y
= (x
- xn
* mp1
) - xn
* mp2
;
283 int4 n
= v
.i
[LOW_HALF
] & 3;
284 double db
= xn
* mp3
;
294 /* Compute sin (A + DA). cos can be computed by shifting the quadrant N
298 do_sincos_1 (double a
, double da
, double x
, int4 n
, int4 k
)
300 double xx
, retval
, res
, cor
, y
;
303 double eps
= fabs (x
) * 1.2e-30;
305 int k1
= (n
+ k
) & 3;
307 { /* quarter of unit circle */
316 res
= TAYLOR_SIN (xx
, a
, da
, cor
);
317 cor
= (cor
> 0) ? 1.02 * cor
+ eps
: 1.02 * cor
- eps
;
318 retval
= (res
== res
+ cor
) ? res
: sloww (a
, da
, x
, k
);
332 res
= do_sin (u
, y
, da
, &cor
);
333 cor
= (cor
> 0) ? 1.035 * cor
+ eps
: 1.035 * cor
- eps
;
334 retval
= ((res
== res
+ cor
) ? ((m
) ? res
: -res
)
335 : sloww1 (a
, da
, x
, m
, k
));
347 y
= a
- (u
.x
- big
) + da
;
348 res
= do_cos (u
, y
, &cor
);
349 cor
= (cor
> 0) ? 1.025 * cor
+ eps
: 1.025 * cor
- eps
;
350 retval
= ((res
== res
+ cor
) ? ((k1
& 2) ? -res
: res
)
351 : sloww2 (a
, da
, x
, n
));
360 reduce_sincos_2 (double x
, double *a
, double *da
)
364 double t
= (x
* hpinv
+ toint
);
365 double xn
= t
- toint
;
367 double xn1
= (xn
+ 8.0e22
) - 8.0e22
;
368 double xn2
= xn
- xn1
;
369 double y
= ((((x
- xn1
* mp1
) - xn1
* mp2
) - xn2
* mp1
) - xn2
* mp2
);
370 int4 n
= v
.i
[LOW_HALF
] & 3;
371 double db
= xn1
* pp3
;
374 db
= (db
- xn2
* pp3
) - xn
* pp4
;
384 /* Compute sin (A + DA). cos can be computed by shifting the quadrant N
388 do_sincos_2 (double a
, double da
, double x
, int4 n
, int4 k
)
390 double res
, retval
, cor
, xx
;
393 double eps
= 1.0e-24;
408 res
= TAYLOR_SIN (xx
, a
, da
, cor
);
409 cor
= (cor
> 0) ? 1.02 * cor
+ eps
: 1.02 * cor
- eps
;
410 retval
= (res
== res
+ cor
) ? res
: bsloww (a
, da
, x
, n
);
430 res
= do_sin (u
, y
, db
, &cor
);
431 cor
= (cor
> 0) ? 1.035 * cor
+ eps
: 1.035 * cor
- eps
;
432 retval
= ((res
== res
+ cor
) ? ((m
) ? res
: -res
)
433 : bsloww1 (a
, da
, x
, n
));
445 double y
= a
- (u
.x
- big
) + da
;
446 res
= do_cos (u
, y
, &cor
);
447 cor
= (cor
> 0) ? 1.025 * cor
+ eps
: 1.025 * cor
- eps
;
448 retval
= ((res
== res
+ cor
) ? ((n
& 2) ? -res
: res
)
449 : bsloww2 (a
, da
, x
, n
));
456 /*******************************************************************/
457 /* An ultimate sin routine. Given an IEEE double machine number x */
458 /* it computes the correctly rounded (to nearest) value of sin(x) */
459 /*******************************************************************/
468 double xx
, res
, t
, cor
, y
, s
, c
, sn
, ssn
, cs
, ccs
;
474 SET_RESTORE_ROUND_53BIT (FE_TONEAREST
);
479 k
= 0x7fffffff & m
; /* no sign */
480 if (k
< 0x3e500000) /* if x->0 =>sin(x)=x */
482 math_check_force_underflow (x
);
485 /*---------------------------- 2^-26 < |x|< 0.25 ----------------------*/
486 else if (k
< 0x3fd00000)
490 t
= POLYNOMIAL (xx
) * (xx
* x
);
493 retval
= (res
== res
+ 1.07 * cor
) ? res
: slow (x
);
494 } /* else if (k < 0x3fd00000) */
495 /*---------------------------- 0.25<|x|< 0.855469---------------------- */
496 else if (k
< 0x3feb6000)
498 u
.x
= (m
> 0) ? big
+ x
: big
- x
;
499 y
= (m
> 0) ? x
- (u
.x
- big
) : x
+ (u
.x
- big
);
501 s
= y
+ y
* xx
* (sn3
+ xx
* sn5
);
502 c
= xx
* (cs2
+ xx
* (cs4
+ xx
* cs6
));
503 SINCOS_TABLE_LOOKUP (u
, sn
, ssn
, cs
, ccs
);
509 cor
= (ssn
+ s
* ccs
- sn
* c
) + cs
* s
;
511 cor
= (sn
- res
) + cor
;
512 retval
= (res
== res
+ 1.096 * cor
) ? res
: slow1 (x
);
513 } /* else if (k < 0x3feb6000) */
515 /*----------------------- 0.855469 <|x|<2.426265 ----------------------*/
516 else if (k
< 0x400368fd)
519 y
= (m
> 0) ? hp0
- x
: hp0
+ x
;
523 y
= (y
- (u
.x
- big
)) + hp1
;
528 y
= (-hp1
) - (y
+ (u
.x
- big
));
530 res
= do_cos (u
, y
, &cor
);
531 retval
= (res
== res
+ 1.020 * cor
) ? ((m
> 0) ? res
: -res
) : slow2 (x
);
532 } /* else if (k < 0x400368fd) */
535 /*-------------------------- 2.426265<|x|< 105414350 ----------------------*/
536 else if (k
< 0x419921FB)
539 int4 n
= reduce_sincos_1 (x
, &a
, &da
);
540 retval
= do_sincos_1 (a
, da
, x
, n
, 0);
541 } /* else if (k < 0x419921FB ) */
543 /*---------------------105414350 <|x|< 281474976710656 --------------------*/
544 else if (k
< 0x42F00000)
548 int4 n
= reduce_sincos_2 (x
, &a
, &da
);
549 retval
= do_sincos_2 (a
, da
, x
, n
, 0);
550 } /* else if (k < 0x42F00000 ) */
552 /* -----------------281474976710656 <|x| <2^1024----------------------------*/
553 else if (k
< 0x7ff00000)
554 retval
= reduce_and_compute (x
, 0);
556 /*--------------------- |x| > 2^1024 ----------------------------------*/
559 if (k
== 0x7ff00000 && u
.i
[LOW_HALF
] == 0)
569 /*******************************************************************/
570 /* An ultimate cos routine. Given an IEEE double machine number x */
571 /* it computes the correctly rounded (to nearest) value of cos(x) */
572 /*******************************************************************/
582 double y
, xx
, res
, cor
, a
, da
;
589 SET_RESTORE_ROUND_53BIT (FE_TONEAREST
);
596 /* |x|<2^-27 => cos(x)=1 */
600 else if (k
< 0x3feb6000)
601 { /* 2^-27 < |x| < 0.855469 */
605 res
= do_cos (u
, y
, &cor
);
606 retval
= (res
== res
+ 1.020 * cor
) ? res
: cslow2 (x
);
607 } /* else if (k < 0x3feb6000) */
609 else if (k
< 0x400368fd)
610 { /* 0.855469 <|x|<2.426265 */ ;
617 res
= TAYLOR_SIN (xx
, a
, da
, cor
);
618 cor
= (cor
> 0) ? 1.02 * cor
+ 1.0e-31 : 1.02 * cor
- 1.0e-31;
619 retval
= (res
== res
+ cor
) ? res
: sloww (a
, da
, x
, 1);
635 res
= do_sin (u
, y
, da
, &cor
);
636 cor
= (cor
> 0) ? 1.035 * cor
+ 1.0e-31 : 1.035 * cor
- 1.0e-31;
637 retval
= ((res
== res
+ cor
) ? ((m
) ? res
: -res
)
638 : sloww1 (a
, da
, x
, m
, 1));
641 } /* else if (k < 0x400368fd) */
645 else if (k
< 0x419921FB)
646 { /* 2.426265<|x|< 105414350 */
648 int4 n
= reduce_sincos_1 (x
, &a
, &da
);
649 retval
= do_sincos_1 (a
, da
, x
, n
, 1);
650 } /* else if (k < 0x419921FB ) */
652 else if (k
< 0x42F00000)
656 int4 n
= reduce_sincos_2 (x
, &a
, &da
);
657 retval
= do_sincos_2 (a
, da
, x
, n
, 1);
658 } /* else if (k < 0x42F00000 ) */
660 /* 281474976710656 <|x| <2^1024 */
661 else if (k
< 0x7ff00000)
662 retval
= reduce_and_compute (x
, 1);
666 if (k
== 0x7ff00000 && u
.i
[LOW_HALF
] == 0)
668 retval
= x
/ x
; /* |x| > 2^1024 */
675 /************************************************************************/
676 /* Routine compute sin(x) for 2^-26 < |x|< 0.25 by Taylor with more */
677 /* precision and if still doesn't accurate enough by mpsin or dubsin */
678 /************************************************************************/
684 double res
, cor
, w
[2];
685 res
= TAYLOR_SLOW (x
, 0, cor
);
686 if (res
== res
+ 1.0007 * cor
)
689 __dubsin (fabs (x
), 0, w
);
690 if (w
[0] == w
[0] + 1.000000001 * w
[1])
691 return (x
> 0) ? w
[0] : -w
[0];
693 return (x
> 0) ? __mpsin (x
, 0, false) : -__mpsin (-x
, 0, false);
696 /*******************************************************************************/
697 /* Routine compute sin(x) for 0.25<|x|< 0.855469 by __sincostab.tbl and Taylor */
698 /* and if result still doesn't accurate enough by mpsin or dubsin */
699 /*******************************************************************************/
706 double w
[2], y
, cor
, res
;
710 res
= do_sin_slow (u
, y
, 0, 0, &cor
);
711 if (res
== res
+ cor
)
712 return (x
> 0) ? res
: -res
;
714 __dubsin (fabs (x
), 0, w
);
715 if (w
[0] == w
[0] + 1.000000005 * w
[1])
716 return (x
> 0) ? w
[0] : -w
[0];
718 return (x
> 0) ? __mpsin (x
, 0, false) : -__mpsin (-x
, 0, false);
721 /**************************************************************************/
722 /* Routine compute sin(x) for 0.855469 <|x|<2.426265 by __sincostab.tbl */
723 /* and if result still doesn't accurate enough by mpsin or dubsin */
724 /**************************************************************************/
730 double w
[2], y
, y1
, y2
, cor
, res
, del
;
743 y
= -(y
+ (u
.x
- big
));
746 res
= do_cos_slow (u
, y
, del
, 0, &cor
);
747 if (res
== res
+ cor
)
748 return (x
> 0) ? res
: -res
;
754 if (w
[0] == w
[0] + 1.000000005 * w
[1])
755 return (x
> 0) ? w
[0] : -w
[0];
757 return (x
> 0) ? __mpsin (x
, 0, false) : -__mpsin (-x
, 0, false);
760 /***************************************************************************/
761 /* Routine compute sin(x+dx) (Double-Length number) where x is small enough*/
762 /* to use Taylor series around zero and (x+dx) */
763 /* in first or third quarter of unit circle.Routine receive also */
764 /* (right argument) the original value of x for computing error of */
765 /* result.And if result not accurate enough routine calls mpsin1 or dubsin */
766 /***************************************************************************/
770 sloww (double x
, double dx
, double orig
, int k
)
772 double y
, t
, res
, cor
, w
[2], a
, da
, xn
;
775 res
= TAYLOR_SLOW (x
, dx
, cor
);
778 cor
= 1.0005 * cor
+ fabs (orig
) * 3.1e-30;
780 cor
= 1.0005 * cor
- fabs (orig
) * 3.1e-30;
782 if (res
== res
+ cor
)
785 (x
> 0) ? __dubsin (x
, dx
, w
) : __dubsin (-x
, -dx
, w
);
787 cor
= 1.000000001 * w
[1] + fabs (orig
) * 1.1e-30;
789 cor
= 1.000000001 * w
[1] - fabs (orig
) * 1.1e-30;
791 if (w
[0] == w
[0] + cor
)
792 return (x
> 0) ? w
[0] : -w
[0];
794 t
= (orig
* hpinv
+ toint
);
797 y
= (orig
- xn
* mp1
) - xn
* mp2
;
798 n
= (v
.i
[LOW_HALF
] + k
) & 3;
804 da
= ((t
- a
) - y
) + da
;
806 if (n
== 2 || n
== 1)
811 (a
> 0) ? __dubsin (a
, da
, w
) : __dubsin (-a
, -da
, w
);
813 cor
= 1.000000001 * w
[1] + fabs (orig
) * 1.1e-40;
815 cor
= 1.000000001 * w
[1] - fabs (orig
) * 1.1e-40;
817 if (w
[0] == w
[0] + cor
)
818 return (a
> 0) ? w
[0] : -w
[0];
820 return (n
& 1) ? __mpcos (orig
, 0, true) : __mpsin (orig
, 0, true);
823 /***************************************************************************/
824 /* Routine compute sin(x+dx) (Double-Length number) where x in first or */
825 /* third quarter of unit circle.Routine receive also (right argument) the */
826 /* original value of x for computing error of result.And if result not */
827 /* accurate enough routine calls mpsin1 or dubsin */
828 /***************************************************************************/
832 sloww1 (double x
, double dx
, double orig
, int m
, int k
)
835 double w
[2], y
, cor
, res
;
839 res
= do_sin_slow (u
, y
, dx
, 3.1e-30 * fabs (orig
), &cor
);
841 if (res
== res
+ cor
)
842 return (m
> 0) ? res
: -res
;
847 cor
= 1.000000005 * w
[1] + 1.1e-30 * fabs (orig
);
849 cor
= 1.000000005 * w
[1] - 1.1e-30 * fabs (orig
);
851 if (w
[0] == w
[0] + cor
)
852 return (m
> 0) ? w
[0] : -w
[0];
854 return (k
== 1) ? __mpcos (orig
, 0, true) : __mpsin (orig
, 0, true);
857 /***************************************************************************/
858 /* Routine compute sin(x+dx) (Double-Length number) where x in second or */
859 /* fourth quarter of unit circle.Routine receive also the original value */
860 /* and quarter(n= 1or 3)of x for computing error of result.And if result not*/
861 /* accurate enough routine calls mpsin1 or dubsin */
862 /***************************************************************************/
866 sloww2 (double x
, double dx
, double orig
, int n
)
869 double w
[2], y
, cor
, res
;
873 res
= do_cos_slow (u
, y
, dx
, 3.1e-30 * fabs (orig
), &cor
);
875 if (res
== res
+ cor
)
876 return (n
& 2) ? -res
: res
;
881 cor
= 1.000000005 * w
[1] + 1.1e-30 * fabs (orig
);
883 cor
= 1.000000005 * w
[1] - 1.1e-30 * fabs (orig
);
885 if (w
[0] == w
[0] + cor
)
886 return (n
& 2) ? -w
[0] : w
[0];
888 return (n
& 1) ? __mpsin (orig
, 0, true) : __mpcos (orig
, 0, true);
891 /***************************************************************************/
892 /* Routine compute sin(x+dx) or cos(x+dx) (Double-Length number) where x */
893 /* is small enough to use Taylor series around zero and (x+dx) */
894 /* in first or third quarter of unit circle.Routine receive also */
895 /* (right argument) the original value of x for computing error of */
896 /* result.And if result not accurate enough routine calls other routines */
897 /***************************************************************************/
901 bsloww (double x
, double dx
, double orig
, int n
)
903 double res
, cor
, w
[2];
905 res
= TAYLOR_SLOW (x
, dx
, cor
);
906 cor
= (cor
> 0) ? 1.0005 * cor
+ 1.1e-24 : 1.0005 * cor
- 1.1e-24;
907 if (res
== res
+ cor
)
910 (x
> 0) ? __dubsin (x
, dx
, w
) : __dubsin (-x
, -dx
, w
);
912 cor
= 1.000000001 * w
[1] + 1.1e-24;
914 cor
= 1.000000001 * w
[1] - 1.1e-24;
916 if (w
[0] == w
[0] + cor
)
917 return (x
> 0) ? w
[0] : -w
[0];
919 return (n
& 1) ? __mpcos (orig
, 0, true) : __mpsin (orig
, 0, true);
922 /***************************************************************************/
923 /* Routine compute sin(x+dx) or cos(x+dx) (Double-Length number) where x */
924 /* in first or third quarter of unit circle.Routine receive also */
925 /* (right argument) the original value of x for computing error of result.*/
926 /* And if result not accurate enough routine calls other routines */
927 /***************************************************************************/
931 bsloww1 (double x
, double dx
, double orig
, int n
)
934 double w
[2], y
, cor
, res
;
939 dx
= (x
> 0) ? dx
: -dx
;
940 res
= do_sin_slow (u
, y
, dx
, 1.1e-24, &cor
);
941 if (res
== res
+ cor
)
942 return (x
> 0) ? res
: -res
;
944 __dubsin (fabs (x
), dx
, w
);
947 cor
= 1.000000005 * w
[1] + 1.1e-24;
949 cor
= 1.000000005 * w
[1] - 1.1e-24;
951 if (w
[0] == w
[0] + cor
)
952 return (x
> 0) ? w
[0] : -w
[0];
954 return (n
& 1) ? __mpcos (orig
, 0, true) : __mpsin (orig
, 0, true);
957 /***************************************************************************/
958 /* Routine compute sin(x+dx) or cos(x+dx) (Double-Length number) where x */
959 /* in second or fourth quarter of unit circle.Routine receive also the */
960 /* original value and quarter(n= 1or 3)of x for computing error of result. */
961 /* And if result not accurate enough routine calls other routines */
962 /***************************************************************************/
966 bsloww2 (double x
, double dx
, double orig
, int n
)
969 double w
[2], y
, cor
, res
;
974 dx
= (x
> 0) ? dx
: -dx
;
975 res
= do_cos_slow (u
, y
, dx
, 1.1e-24, &cor
);
976 if (res
== res
+ cor
)
977 return (n
& 2) ? -res
: res
;
979 __docos (fabs (x
), dx
, w
);
982 cor
= 1.000000005 * w
[1] + 1.1e-24;
984 cor
= 1.000000005 * w
[1] - 1.1e-24;
986 if (w
[0] == w
[0] + cor
)
987 return (n
& 2) ? -w
[0] : w
[0];
989 return (n
& 1) ? __mpsin (orig
, 0, true) : __mpcos (orig
, 0, true);
992 /************************************************************************/
993 /* Routine compute cos(x) for 2^-27 < |x|< 0.25 by Taylor with more */
994 /* precision and if still doesn't accurate enough by mpcos or docos */
995 /************************************************************************/
1002 double w
[2], y
, cor
, res
;
1006 y
= y
- (u
.x
- big
);
1007 res
= do_cos_slow (u
, y
, 0, 0, &cor
);
1008 if (res
== res
+ cor
)
1013 if (w
[0] == w
[0] + 1.000000005 * w
[1])
1016 return __mpcos (x
, 0, false);
1020 weak_alias (__cos
, cos
)
1021 # ifdef NO_LONG_DOUBLE
1022 strong_alias (__cos
, __cosl
)
1023 weak_alias (__cos
, cosl
)
1027 weak_alias (__sin
, sin
)
1028 # ifdef NO_LONG_DOUBLE
1029 strong_alias (__sin
, __sinl
)
1030 weak_alias (__sin
, sinl
)