]>
git.ipfire.org Git - thirdparty/gcc.git/blob - libgcc-math/dbl-64/s_sin.c
2 * IBM Accurate Mathematical Library
3 * written by International Business Machines Corp.
4 * Copyright (C) 2001 Free Software Foundation
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, write to the Free Software
18 * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
20 /****************************************************************************/
22 /* MODULE_NAME:usncs.c */
39 /* FILES NEEDED: dla.h endian.h mpa.h mydefs.h usncs.h */
40 /* branred.c sincos32.c dosincos.c mpa.c */
43 /* An ultimate sin and routine. Given an IEEE double machine number x */
44 /* it computes the correctly rounded (to nearest) value of sin(x) or cos(x) */
45 /* Assumption: Machine arithmetic operations are performed in */
46 /* round to nearest mode of IEEE 754 standard. */
48 /****************************************************************************/
56 #include "math_private.h"
59 sn3
= -1.66666666666664880952546298448555E-01,
60 sn5
= 8.33333214285722277379541354343671E-03,
61 cs2
= 4.99999999999999999999950396842453E-01,
62 cs4
= -4.16666666666664434524222570944589E-02,
63 cs6
= 1.38888874007937613028114285595617E-03;
65 void __dubsin(double x
, double dx
, double w
[]);
66 void __docos(double x
, double dx
, double w
[]);
67 double __mpsin(double x
, double dx
);
68 double __mpcos(double x
, double dx
);
69 double __mpsin1(double x
);
70 double __mpcos1(double x
);
71 static double slow(double x
);
72 static double slow1(double x
);
73 static double slow2(double x
);
74 static double sloww(double x
, double dx
, double orig
);
75 static double sloww1(double x
, double dx
, double orig
);
76 static double sloww2(double x
, double dx
, double orig
, int n
);
77 static double bsloww(double x
, double dx
, double orig
, int n
);
78 static double bsloww1(double x
, double dx
, double orig
, int n
);
79 static double bsloww2(double x
, double dx
, double orig
, int n
);
80 int __branred(double x
, double *a
, double *aa
);
81 static double cslow2(double x
);
82 static double csloww(double x
, double dx
, double orig
);
83 static double csloww1(double x
, double dx
, double orig
);
84 static double csloww2(double x
, double dx
, double orig
, int n
);
85 /*******************************************************************/
86 /* An ultimate sin routine. Given an IEEE double machine number x */
87 /* it computes the correctly rounded (to nearest) value of sin(x) */
88 /*******************************************************************/
89 double __sin(double x
){
90 double xx
,res
,t
,cor
,y
,s
,c
,sn
,ssn
,cs
,ccs
,xn
,a
,da
,db
,eps
,xn1
,xn2
;
102 k
= 0x7fffffff&m
; /* no sign */
103 if (k
< 0x3e500000) /* if x->0 =>sin(x)=x */
105 /*---------------------------- 2^-26 < |x|< 0.25 ----------------------*/
106 else if (k
< 0x3fd00000){
109 t
= ((((s5
.x
*xx
+ s4
.x
)*xx
+ s3
.x
)*xx
+ s2
.x
)*xx
+ s1
.x
)*(xx
*x
);
112 return (res
== res
+ 1.07*cor
)? res
: slow(x
);
113 } /* else if (k < 0x3fd00000) */
114 /*---------------------------- 0.25<|x|< 0.855469---------------------- */
115 else if (k
< 0x3feb6000) {
116 u
.x
=(m
>0)?big
.x
+x
:big
.x
-x
;
117 y
=(m
>0)?x
-(u
.x
-big
.x
):x
+(u
.x
-big
.x
);
119 s
= y
+ y
*xx
*(sn3
+xx
*sn5
);
120 c
= xx
*(cs2
+xx
*(cs4
+ xx
*cs6
));
122 sn
=(m
>0)?sincos
.x
[k
]:-sincos
.x
[k
];
123 ssn
=(m
>0)?sincos
.x
[k
+1]:-sincos
.x
[k
+1];
126 cor
=(ssn
+s
*ccs
-sn
*c
)+cs
*s
;
129 return (res
==res
+1.025*cor
)? res
: slow1(x
);
130 } /* else if (k < 0x3feb6000) */
132 /*----------------------- 0.855469 <|x|<2.426265 ----------------------*/
133 else if (k
< 0x400368fd ) {
135 y
= (m
>0)? hp0
.x
-x
:hp0
.x
+x
;
138 y
= (y
-(u
.x
-big
.x
))+hp1
.x
;
142 y
= (-hp1
.x
) - (y
+(u
.x
-big
.x
));
145 s
= y
+ y
*xx
*(sn3
+xx
*sn5
);
146 c
= xx
*(cs2
+xx
*(cs4
+ xx
*cs6
));
152 cor
=(ccs
-s
*ssn
-cs
*c
)-sn
*s
;
155 return (res
==res
+1.020*cor
)? ((m
>0)?res
:-res
) : slow2(x
);
156 } /* else if (k < 0x400368fd) */
158 /*-------------------------- 2.426265<|x|< 105414350 ----------------------*/
159 else if (k
< 0x419921FB ) {
160 t
= (x
*hpinv
.x
+ toint
.x
);
163 y
= (x
- xn
*mp1
.x
) - xn
*mp2
.x
;
168 eps
= ABS(x
)*1.2e-30;
170 switch (n
) { /* quarter of unit circle */
174 if (n
) {a
=-a
;da
=-da
;}
177 t
= (((((s5
.x
*xx
+ s4
.x
)*xx
+ s3
.x
)*xx
+ s2
.x
)*xx
+ s1
.x
)*a
- 0.5*da
)*xx
+da
;
180 cor
= (cor
>0)? 1.02*cor
+eps
: 1.02*cor
-eps
;
181 return (res
== res
+ cor
)? res
: sloww(a
,da
,x
);
191 s
= y
+ (db
+y
*xx
*(sn3
+xx
*sn5
));
192 c
= y
*db
+xx
*(cs2
+xx
*(cs4
+ xx
*cs6
));
198 cor
=(ssn
+s
*ccs
-sn
*c
)+cs
*s
;
201 cor
= (cor
>0)? 1.035*cor
+eps
: 1.035*cor
-eps
;
202 return (res
==res
+cor
)? ((m
)?res
:-res
) : sloww1(a
,da
,x
);
218 s
= y
+ y
*xx
*(sn3
+xx
*sn5
);
219 c
= xx
*(cs2
+xx
*(cs4
+ xx
*cs6
));
220 cor
=(ccs
-s
*ssn
-cs
*c
)-sn
*s
;
223 cor
= (cor
>0)? 1.025*cor
+eps
: 1.025*cor
-eps
;
224 return (res
==res
+cor
)? ((n
&2)?-res
:res
) : sloww2(a
,da
,x
,n
);
230 } /* else if (k < 0x419921FB ) */
232 /*---------------------105414350 <|x|< 281474976710656 --------------------*/
233 else if (k
< 0x42F00000 ) {
234 t
= (x
*hpinv
.x
+ toint
.x
);
237 xn1
= (xn
+8.0e22
)-8.0e22
;
239 y
= ((((x
- xn1
*mp1
.x
) - xn1
*mp2
.x
)-xn2
*mp1
.x
)-xn2
*mp2
.x
);
244 da
= (da
- xn2
*pp3
.x
) -xn
*pp4
.x
;
253 if (n
) {a
=-a
;da
=-da
;}
256 t
= (((((s5
.x
*xx
+ s4
.x
)*xx
+ s3
.x
)*xx
+ s2
.x
)*xx
+ s1
.x
)*a
- 0.5*da
)*xx
+da
;
259 cor
= (cor
>0)? 1.02*cor
+eps
: 1.02*cor
-eps
;
260 return (res
== res
+ cor
)? res
: bsloww(a
,da
,x
,n
);
263 if (a
>0) {m
=1;t
=a
;db
=da
;}
264 else {m
=0;t
=-a
;db
=-da
;}
268 s
= y
+ (db
+y
*xx
*(sn3
+xx
*sn5
));
269 c
= y
*db
+xx
*(cs2
+xx
*(cs4
+ xx
*cs6
));
275 cor
=(ssn
+s
*ccs
-sn
*c
)+cs
*s
;
278 cor
= (cor
>0)? 1.035*cor
+eps
: 1.035*cor
-eps
;
279 return (res
==res
+cor
)? ((m
)?res
:-res
) : bsloww1(a
,da
,x
,n
);
295 s
= y
+ y
*xx
*(sn3
+xx
*sn5
);
296 c
= xx
*(cs2
+xx
*(cs4
+ xx
*cs6
));
297 cor
=(ccs
-s
*ssn
-cs
*c
)-sn
*s
;
300 cor
= (cor
>0)? 1.025*cor
+eps
: 1.025*cor
-eps
;
301 return (res
==res
+cor
)? ((n
&2)?-res
:res
) : bsloww2(a
,da
,x
,n
);
307 } /* else if (k < 0x42F00000 ) */
309 /* -----------------281474976710656 <|x| <2^1024----------------------------*/
310 else if (k
< 0x7ff00000) {
312 n
= __branred(x
,&a
,&da
);
315 if (a
*a
< 0.01588) return bsloww(a
,da
,x
,n
);
316 else return bsloww1(a
,da
,x
,n
);
319 if (a
*a
< 0.01588) return bsloww(-a
,-da
,x
,n
);
320 else return bsloww1(-a
,-da
,x
,n
);
325 return bsloww2(a
,da
,x
,n
);
329 } /* else if (k < 0x7ff00000 ) */
331 /*--------------------- |x| > 2^1024 ----------------------------------*/
333 return 0; /* unreachable */
337 /*******************************************************************/
338 /* An ultimate cos routine. Given an IEEE double machine number x */
339 /* it computes the correctly rounded (to nearest) value of cos(x) */
340 /*******************************************************************/
342 double __cos(double x
)
344 double y
,xx
,res
,t
,cor
,s
,c
,sn
,ssn
,cs
,ccs
,xn
,a
,da
,db
,eps
,xn1
,xn2
;
352 if (k
< 0x3e400000 ) return 1.0; /* |x|<2^-27 => cos(x)=1 */
354 else if (k
< 0x3feb6000 ) {/* 2^-27 < |x| < 0.855469 */
359 s
= y
+ y
*xx
*(sn3
+xx
*sn5
);
360 c
= xx
*(cs2
+xx
*(cs4
+ xx
*cs6
));
366 cor
=(ccs
-s
*ssn
-cs
*c
)-sn
*s
;
369 return (res
==res
+1.020*cor
)? res
: cslow2(x
);
371 } /* else if (k < 0x3feb6000) */
373 else if (k
< 0x400368fd ) {/* 0.855469 <|x|<2.426265 */;
379 t
= (((((s5
.x
*xx
+ s4
.x
)*xx
+ s3
.x
)*xx
+ s2
.x
)*xx
+ s1
.x
)*a
- 0.5*da
)*xx
+da
;
382 cor
= (cor
>0)? 1.02*cor
+1.0e-31 : 1.02*cor
-1.0e-31;
383 return (res
== res
+ cor
)? res
: csloww(a
,da
,x
);
386 if (a
>0) {m
=1;t
=a
;db
=da
;}
387 else {m
=0;t
=-a
;db
=-da
;}
391 s
= y
+ (db
+y
*xx
*(sn3
+xx
*sn5
));
392 c
= y
*db
+xx
*(cs2
+xx
*(cs4
+ xx
*cs6
));
398 cor
=(ssn
+s
*ccs
-sn
*c
)+cs
*s
;
401 cor
= (cor
>0)? 1.035*cor
+1.0e-31 : 1.035*cor
-1.0e-31;
402 return (res
==res
+cor
)? ((m
)?res
:-res
) : csloww1(a
,da
,x
);
405 } /* else if (k < 0x400368fd) */
408 else if (k
< 0x419921FB ) {/* 2.426265<|x|< 105414350 */
409 t
= (x
*hpinv
.x
+ toint
.x
);
412 y
= (x
- xn
*mp1
.x
) - xn
*mp2
.x
;
417 eps
= ABS(x
)*1.2e-30;
423 if (n
== 1) {a
=-a
;da
=-da
;}
425 t
= (((((s5
.x
*xx
+ s4
.x
)*xx
+ s3
.x
)*xx
+ s2
.x
)*xx
+ s1
.x
)*a
- 0.5*da
)*xx
+da
;
428 cor
= (cor
>0)? 1.02*cor
+eps
: 1.02*cor
-eps
;
429 return (res
== res
+ cor
)? res
: csloww(a
,da
,x
);
432 if (a
>0) {m
=1;t
=a
;db
=da
;}
433 else {m
=0;t
=-a
;db
=-da
;}
437 s
= y
+ (db
+y
*xx
*(sn3
+xx
*sn5
));
438 c
= y
*db
+xx
*(cs2
+xx
*(cs4
+ xx
*cs6
));
444 cor
=(ssn
+s
*ccs
-sn
*c
)+cs
*s
;
447 cor
= (cor
>0)? 1.035*cor
+eps
: 1.035*cor
-eps
;
448 return (res
==res
+cor
)? ((m
)?res
:-res
) : csloww1(a
,da
,x
);
454 if (a
<0) {a
=-a
;da
=-da
;}
463 s
= y
+ y
*xx
*(sn3
+xx
*sn5
);
464 c
= xx
*(cs2
+xx
*(cs4
+ xx
*cs6
));
465 cor
=(ccs
-s
*ssn
-cs
*c
)-sn
*s
;
468 cor
= (cor
>0)? 1.025*cor
+eps
: 1.025*cor
-eps
;
469 return (res
==res
+cor
)? ((n
)?-res
:res
) : csloww2(a
,da
,x
,n
);
475 } /* else if (k < 0x419921FB ) */
478 else if (k
< 0x42F00000 ) {
479 t
= (x
*hpinv
.x
+ toint
.x
);
482 xn1
= (xn
+8.0e22
)-8.0e22
;
484 y
= ((((x
- xn1
*mp1
.x
) - xn1
*mp2
.x
)-xn2
*mp1
.x
)-xn2
*mp2
.x
);
489 da
= (da
- xn2
*pp3
.x
) -xn
*pp4
.x
;
498 if (n
==1) {a
=-a
;da
=-da
;}
500 t
= (((((s5
.x
*xx
+ s4
.x
)*xx
+ s3
.x
)*xx
+ s2
.x
)*xx
+ s1
.x
)*a
- 0.5*da
)*xx
+da
;
503 cor
= (cor
>0)? 1.02*cor
+eps
: 1.02*cor
-eps
;
504 return (res
== res
+ cor
)? res
: bsloww(a
,da
,x
,n
);
507 if (a
>0) {m
=1;t
=a
;db
=da
;}
508 else {m
=0;t
=-a
;db
=-da
;}
512 s
= y
+ (db
+y
*xx
*(sn3
+xx
*sn5
));
513 c
= y
*db
+xx
*(cs2
+xx
*(cs4
+ xx
*cs6
));
519 cor
=(ssn
+s
*ccs
-sn
*c
)+cs
*s
;
522 cor
= (cor
>0)? 1.035*cor
+eps
: 1.035*cor
-eps
;
523 return (res
==res
+cor
)? ((m
)?res
:-res
) : bsloww1(a
,da
,x
,n
);
529 if (a
<0) {a
=-a
;da
=-da
;}
538 s
= y
+ y
*xx
*(sn3
+xx
*sn5
);
539 c
= xx
*(cs2
+xx
*(cs4
+ xx
*cs6
));
540 cor
=(ccs
-s
*ssn
-cs
*c
)-sn
*s
;
543 cor
= (cor
>0)? 1.025*cor
+eps
: 1.025*cor
-eps
;
544 return (res
==res
+cor
)? ((n
)?-res
:res
) : bsloww2(a
,da
,x
,n
);
549 } /* else if (k < 0x42F00000 ) */
551 else if (k
< 0x7ff00000) {/* 281474976710656 <|x| <2^1024 */
553 n
= __branred(x
,&a
,&da
);
556 if (a
*a
< 0.01588) return bsloww(-a
,-da
,x
,n
);
557 else return bsloww1(-a
,-da
,x
,n
);
560 if (a
*a
< 0.01588) return bsloww(a
,da
,x
,n
);
561 else return bsloww1(a
,da
,x
,n
);
566 return bsloww2(a
,da
,x
,n
);
570 } /* else if (k < 0x7ff00000 ) */
575 else return x
/ x
; /* |x| > 2^1024 */
580 /************************************************************************/
581 /* Routine compute sin(x) for 2^-26 < |x|< 0.25 by Taylor with more */
582 /* precision and if still doesn't accurate enough by mpsin or dubsin */
583 /************************************************************************/
585 static double slow(double x
) {
586 static const double th2_36
= 206158430208.0; /* 1.5*2**37 */
587 double y
,x1
,x2
,xx
,r
,t
,res
,cor
,w
[2];
588 x1
=(x
+th2_36
)-th2_36
;
593 t
= (((((s5
.x
*xx
+ s4
.x
)*xx
+ s3
.x
)*xx
+ s2
.x
)*xx
+ bb
.x
)*xx
+ 3.0*aa
.x
*x1
*x2
)*x
+aa
.x
*x2
*x2
*x2
;
597 if (res
== res
+ 1.0007*cor
) return res
;
599 __dubsin(ABS(x
),0,w
);
600 if (w
[0] == w
[0]+1.000000001*w
[1]) return (x
>0)?w
[0]:-w
[0];
601 else return (x
>0)?__mpsin(x
,0):-__mpsin(-x
,0);
604 /*******************************************************************************/
605 /* Routine compute sin(x) for 0.25<|x|< 0.855469 by sincos.tbl and Taylor */
606 /* and if result still doesn't accurate enough by mpsin or dubsin */
607 /*******************************************************************************/
609 static double slow1(double x
) {
611 double sn
,ssn
,cs
,ccs
,s
,c
,w
[2],y
,y1
,y2
,c1
,c2
,xx
,cor
,res
;
612 static const double t22
= 6291456.0;
618 s
= y
*xx
*(sn3
+xx
*sn5
);
619 c
= xx
*(cs2
+xx
*(cs4
+ xx
*cs6
));
621 sn
=sincos
.x
[k
]; /* Data */
622 ssn
=sincos
.x
[k
+1]; /* from */
623 cs
=sincos
.x
[k
+2]; /* tables */
624 ccs
=sincos
.x
[k
+3]; /* sincos.tbl */
629 cor
=(ssn
+s
*ccs
+cs
*s
+c2
*y
+c1
*y2
)-sn
*c
;
631 cor
= cor
+((sn
-y
)+c1
*y1
);
634 if (res
== res
+1.0005*cor
) return (x
>0)?res
:-res
;
636 __dubsin(ABS(x
),0,w
);
637 if (w
[0] == w
[0]+1.000000005*w
[1]) return (x
>0)?w
[0]:-w
[0];
638 else return (x
>0)?__mpsin(x
,0):-__mpsin(-x
,0);
641 /**************************************************************************/
642 /* Routine compute sin(x) for 0.855469 <|x|<2.426265 by sincos.tbl */
643 /* and if result still doesn't accurate enough by mpsin or dubsin */
644 /**************************************************************************/
645 static double slow2(double x
) {
647 double sn
,ssn
,cs
,ccs
,s
,c
,w
[2],y
,y1
,y2
,e1
,e2
,xx
,cor
,res
,del
;
648 static const double t22
= 6291456.0;
659 y
= -(y
+(u
.x
-big
.x
));
663 s
= y
*xx
*(sn3
+xx
*sn5
);
664 c
= y
*del
+xx
*(cs2
+xx
*(cs4
+ xx
*cs6
));
674 cor
=(ccs
-cs
*c
-e1
*y2
-e2
*y
)-sn
*s
;
676 cor
= cor
+((cs
-y
)-e1
*y1
);
679 if (res
== res
+1.0005*cor
) return (x
>0)?res
:-res
;
685 if (w
[0] == w
[0]+1.000000005*w
[1]) return (x
>0)?w
[0]:-w
[0];
686 else return (x
>0)?__mpsin(x
,0):-__mpsin(-x
,0);
689 /***************************************************************************/
690 /* Routine compute sin(x+dx) (Double-Length number) where x is small enough*/
691 /* to use Taylor series around zero and (x+dx) */
692 /* in first or third quarter of unit circle.Routine receive also */
693 /* (right argument) the original value of x for computing error of */
694 /* result.And if result not accurate enough routine calls mpsin1 or dubsin */
695 /***************************************************************************/
697 static double sloww(double x
,double dx
, double orig
) {
698 static const double th2_36
= 206158430208.0; /* 1.5*2**37 */
699 double y
,x1
,x2
,xx
,r
,t
,res
,cor
,w
[2],a
,da
,xn
;
700 union {int4 i
[2]; double x
;} v
;
702 x1
=(x
+th2_36
)-th2_36
;
707 t
= (((((s5
.x
*xx
+ s4
.x
)*xx
+ s3
.x
)*xx
+ s2
.x
)*xx
+ bb
.x
)*xx
+ 3.0*aa
.x
*x1
*x2
)*x
+aa
.x
*x2
*x2
*x2
+dx
;
711 cor
= (cor
>0)? 1.0005*cor
+ABS(orig
)*3.1e-30 : 1.0005*cor
-ABS(orig
)*3.1e-30;
712 if (res
== res
+ cor
) return res
;
714 (x
>0)? __dubsin(x
,dx
,w
) : __dubsin(-x
,-dx
,w
);
715 cor
= (w
[1]>0)? 1.000000001*w
[1] + ABS(orig
)*1.1e-30 : 1.000000001*w
[1] - ABS(orig
)*1.1e-30;
716 if (w
[0] == w
[0]+cor
) return (x
>0)?w
[0]:-w
[0];
718 t
= (orig
*hpinv
.x
+ toint
.x
);
721 y
= (orig
- xn
*mp1
.x
) - xn
*mp2
.x
;
729 if (n
&2) {a
=-a
; da
=-da
;}
730 (a
>0)? __dubsin(a
,da
,w
) : __dubsin(-a
,-da
,w
);
731 cor
= (w
[1]>0)? 1.000000001*w
[1] + ABS(orig
)*1.1e-40 : 1.000000001*w
[1] - ABS(orig
)*1.1e-40;
732 if (w
[0] == w
[0]+cor
) return (a
>0)?w
[0]:-w
[0];
733 else return __mpsin1(orig
);
737 /***************************************************************************/
738 /* Routine compute sin(x+dx) (Double-Length number) where x in first or */
739 /* third quarter of unit circle.Routine receive also (right argument) the */
740 /* original value of x for computing error of result.And if result not */
741 /* accurate enough routine calls mpsin1 or dubsin */
742 /***************************************************************************/
744 static double sloww1(double x
, double dx
, double orig
) {
746 double sn
,ssn
,cs
,ccs
,s
,c
,w
[2],y
,y1
,y2
,c1
,c2
,xx
,cor
,res
;
747 static const double t22
= 6291456.0;
754 s
= y
*xx
*(sn3
+xx
*sn5
);
755 c
= xx
*(cs2
+xx
*(cs4
+ xx
*cs6
));
765 cor
=(ssn
+s
*ccs
+cs
*s
+c2
*y
+c1
*y2
-sn
*y
*dx
)-sn
*c
;
767 cor
= cor
+((sn
-y
)+c1
*y1
);
770 cor
= (cor
>0)? 1.0005*cor
+3.1e-30*ABS(orig
) : 1.0005*cor
-3.1e-30*ABS(orig
);
771 if (res
== res
+ cor
) return (x
>0)?res
:-res
;
773 __dubsin(ABS(x
),dx
,w
);
774 cor
= (w
[1]>0)? 1.000000005*w
[1]+1.1e-30*ABS(orig
) : 1.000000005*w
[1]-1.1e-30*ABS(orig
);
775 if (w
[0] == w
[0]+cor
) return (x
>0)?w
[0]:-w
[0];
776 else return __mpsin1(orig
);
779 /***************************************************************************/
780 /* Routine compute sin(x+dx) (Double-Length number) where x in second or */
781 /* fourth quarter of unit circle.Routine receive also the original value */
782 /* and quarter(n= 1or 3)of x for computing error of result.And if result not*/
783 /* accurate enough routine calls mpsin1 or dubsin */
784 /***************************************************************************/
786 static double sloww2(double x
, double dx
, double orig
, int n
) {
788 double sn
,ssn
,cs
,ccs
,s
,c
,w
[2],y
,y1
,y2
,e1
,e2
,xx
,cor
,res
;
789 static const double t22
= 6291456.0;
796 s
= y
*xx
*(sn3
+xx
*sn5
);
797 c
= y
*dx
+xx
*(cs2
+xx
*(cs4
+ xx
*cs6
));
808 cor
=(ccs
-cs
*c
-e1
*y2
-e2
*y
)-sn
*s
;
810 cor
= cor
+((cs
-y
)-e1
*y1
);
813 cor
= (cor
>0)? 1.0005*cor
+3.1e-30*ABS(orig
) : 1.0005*cor
-3.1e-30*ABS(orig
);
814 if (res
== res
+ cor
) return (n
&2)?-res
:res
;
816 __docos(ABS(x
),dx
,w
);
817 cor
= (w
[1]>0)? 1.000000005*w
[1]+1.1e-30*ABS(orig
) : 1.000000005*w
[1]-1.1e-30*ABS(orig
);
818 if (w
[0] == w
[0]+cor
) return (n
&2)?-w
[0]:w
[0];
819 else return __mpsin1(orig
);
822 /***************************************************************************/
823 /* Routine compute sin(x+dx) or cos(x+dx) (Double-Length number) where x */
824 /* is small enough to use Taylor series around zero and (x+dx) */
825 /* in first or third quarter of unit circle.Routine receive also */
826 /* (right argument) the original value of x for computing error of */
827 /* result.And if result not accurate enough routine calls other routines */
828 /***************************************************************************/
830 static double bsloww(double x
,double dx
, double orig
,int n
) {
831 static const double th2_36
= 206158430208.0; /* 1.5*2**37 */
832 double y
,x1
,x2
,xx
,r
,t
,res
,cor
,w
[2];
835 union {int4 i
[2]; double x
;} v
;
837 x1
=(x
+th2_36
)-th2_36
;
842 t
= (((((s5
.x
*xx
+ s4
.x
)*xx
+ s3
.x
)*xx
+ s2
.x
)*xx
+ bb
.x
)*xx
+ 3.0*aa
.x
*x1
*x2
)*x
+aa
.x
*x2
*x2
*x2
+dx
;
846 cor
= (cor
>0)? 1.0005*cor
+1.1e-24 : 1.0005*cor
-1.1e-24;
847 if (res
== res
+ cor
) return res
;
849 (x
>0)? __dubsin(x
,dx
,w
) : __dubsin(-x
,-dx
,w
);
850 cor
= (w
[1]>0)? 1.000000001*w
[1] + 1.1e-24 : 1.000000001*w
[1] - 1.1e-24;
851 if (w
[0] == w
[0]+cor
) return (x
>0)?w
[0]:-w
[0];
852 else return (n
&1)?__mpcos1(orig
):__mpsin1(orig
);
856 /***************************************************************************/
857 /* Routine compute sin(x+dx) or cos(x+dx) (Double-Length number) where x */
858 /* in first or third quarter of unit circle.Routine receive also */
859 /* (right argument) the original value of x for computing error of result.*/
860 /* And if result not accurate enough routine calls other routines */
861 /***************************************************************************/
863 static double bsloww1(double x
, double dx
, double orig
,int n
) {
865 double sn
,ssn
,cs
,ccs
,s
,c
,w
[2],y
,y1
,y2
,c1
,c2
,xx
,cor
,res
;
866 static const double t22
= 6291456.0;
873 s
= y
*xx
*(sn3
+xx
*sn5
);
874 c
= xx
*(cs2
+xx
*(cs4
+ xx
*cs6
));
884 cor
=(ssn
+s
*ccs
+cs
*s
+c2
*y
+c1
*y2
-sn
*y
*dx
)-sn
*c
;
886 cor
= cor
+((sn
-y
)+c1
*y1
);
889 cor
= (cor
>0)? 1.0005*cor
+1.1e-24 : 1.0005*cor
-1.1e-24;
890 if (res
== res
+ cor
) return (x
>0)?res
:-res
;
892 __dubsin(ABS(x
),dx
,w
);
893 cor
= (w
[1]>0)? 1.000000005*w
[1]+1.1e-24: 1.000000005*w
[1]-1.1e-24;
894 if (w
[0] == w
[0]+cor
) return (x
>0)?w
[0]:-w
[0];
895 else return (n
&1)?__mpcos1(orig
):__mpsin1(orig
);
899 /***************************************************************************/
900 /* Routine compute sin(x+dx) or cos(x+dx) (Double-Length number) where x */
901 /* in second or fourth quarter of unit circle.Routine receive also the */
902 /* original value and quarter(n= 1or 3)of x for computing error of result. */
903 /* And if result not accurate enough routine calls other routines */
904 /***************************************************************************/
906 static double bsloww2(double x
, double dx
, double orig
, int n
) {
908 double sn
,ssn
,cs
,ccs
,s
,c
,w
[2],y
,y1
,y2
,e1
,e2
,xx
,cor
,res
;
909 static const double t22
= 6291456.0;
916 s
= y
*xx
*(sn3
+xx
*sn5
);
917 c
= y
*dx
+xx
*(cs2
+xx
*(cs4
+ xx
*cs6
));
928 cor
=(ccs
-cs
*c
-e1
*y2
-e2
*y
)-sn
*s
;
930 cor
= cor
+((cs
-y
)-e1
*y1
);
933 cor
= (cor
>0)? 1.0005*cor
+1.1e-24 : 1.0005*cor
-1.1e-24;
934 if (res
== res
+ cor
) return (n
&2)?-res
:res
;
936 __docos(ABS(x
),dx
,w
);
937 cor
= (w
[1]>0)? 1.000000005*w
[1]+1.1e-24 : 1.000000005*w
[1]-1.1e-24;
938 if (w
[0] == w
[0]+cor
) return (n
&2)?-w
[0]:w
[0];
939 else return (n
&1)?__mpsin1(orig
):__mpcos1(orig
);
943 /************************************************************************/
944 /* Routine compute cos(x) for 2^-27 < |x|< 0.25 by Taylor with more */
945 /* precision and if still doesn't accurate enough by mpcos or docos */
946 /************************************************************************/
948 static double cslow2(double x
) {
950 double sn
,ssn
,cs
,ccs
,s
,c
,w
[2],y
,y1
,y2
,e1
,e2
,xx
,cor
,res
;
951 static const double t22
= 6291456.0;
957 s
= y
*xx
*(sn3
+xx
*sn5
);
958 c
= xx
*(cs2
+xx
*(cs4
+ xx
*cs6
));
968 cor
=(ccs
-cs
*c
-e1
*y2
-e2
*y
)-sn
*s
;
970 cor
= cor
+((cs
-y
)-e1
*y1
);
973 if (res
== res
+1.0005*cor
)
978 if (w
[0] == w
[0]+1.000000005*w
[1]) return w
[0];
979 else return __mpcos(x
,0);
983 /***************************************************************************/
984 /* Routine compute cos(x+dx) (Double-Length number) where x is small enough*/
985 /* to use Taylor series around zero and (x+dx) .Routine receive also */
986 /* (right argument) the original value of x for computing error of */
987 /* result.And if result not accurate enough routine calls other routines */
988 /***************************************************************************/
991 static double csloww(double x
,double dx
, double orig
) {
992 static const double th2_36
= 206158430208.0; /* 1.5*2**37 */
993 double y
,x1
,x2
,xx
,r
,t
,res
,cor
,w
[2],a
,da
,xn
;
994 union {int4 i
[2]; double x
;} v
;
996 x1
=(x
+th2_36
)-th2_36
;
1002 t
= (((((s5
.x
*xx
+ s4
.x
)*xx
+ s3
.x
)*xx
+ s2
.x
)*xx
+ bb
.x
)*xx
+ 3.0*aa
.x
*x1
*x2
)*x
+aa
.x
*x2
*x2
*x2
+dx
;
1006 cor
= (cor
>0)? 1.0005*cor
+ABS(orig
)*3.1e-30 : 1.0005*cor
-ABS(orig
)*3.1e-30;
1007 if (res
== res
+ cor
) return res
;
1009 (x
>0)? __dubsin(x
,dx
,w
) : __dubsin(-x
,-dx
,w
);
1010 cor
= (w
[1]>0)? 1.000000001*w
[1] + ABS(orig
)*1.1e-30 : 1.000000001*w
[1] - ABS(orig
)*1.1e-30;
1011 if (w
[0] == w
[0]+cor
) return (x
>0)?w
[0]:-w
[0];
1013 t
= (orig
*hpinv
.x
+ toint
.x
);
1016 y
= (orig
- xn
*mp1
.x
) - xn
*mp2
.x
;
1024 if (n
==1) {a
=-a
; da
=-da
;}
1025 (a
>0)? __dubsin(a
,da
,w
) : __dubsin(-a
,-da
,w
);
1026 cor
= (w
[1]>0)? 1.000000001*w
[1] + ABS(orig
)*1.1e-40 : 1.000000001*w
[1] - ABS(orig
)*1.1e-40;
1027 if (w
[0] == w
[0]+cor
) return (a
>0)?w
[0]:-w
[0];
1028 else return __mpcos1(orig
);
1033 /***************************************************************************/
1034 /* Routine compute sin(x+dx) (Double-Length number) where x in first or */
1035 /* third quarter of unit circle.Routine receive also (right argument) the */
1036 /* original value of x for computing error of result.And if result not */
1037 /* accurate enough routine calls other routines */
1038 /***************************************************************************/
1040 static double csloww1(double x
, double dx
, double orig
) {
1042 double sn
,ssn
,cs
,ccs
,s
,c
,w
[2],y
,y1
,y2
,c1
,c2
,xx
,cor
,res
;
1043 static const double t22
= 6291456.0;
1050 s
= y
*xx
*(sn3
+xx
*sn5
);
1051 c
= xx
*(cs2
+xx
*(cs4
+ xx
*cs6
));
1061 cor
=(ssn
+s
*ccs
+cs
*s
+c2
*y
+c1
*y2
-sn
*y
*dx
)-sn
*c
;
1063 cor
= cor
+((sn
-y
)+c1
*y1
);
1066 cor
= (cor
>0)? 1.0005*cor
+3.1e-30*ABS(orig
) : 1.0005*cor
-3.1e-30*ABS(orig
);
1067 if (res
== res
+ cor
) return (x
>0)?res
:-res
;
1069 __dubsin(ABS(x
),dx
,w
);
1070 cor
= (w
[1]>0)? 1.000000005*w
[1]+1.1e-30*ABS(orig
) : 1.000000005*w
[1]-1.1e-30*ABS(orig
);
1071 if (w
[0] == w
[0]+cor
) return (x
>0)?w
[0]:-w
[0];
1072 else return __mpcos1(orig
);
1077 /***************************************************************************/
1078 /* Routine compute sin(x+dx) (Double-Length number) where x in second or */
1079 /* fourth quarter of unit circle.Routine receive also the original value */
1080 /* and quarter(n= 1or 3)of x for computing error of result.And if result not*/
1081 /* accurate enough routine calls other routines */
1082 /***************************************************************************/
1084 static double csloww2(double x
, double dx
, double orig
, int n
) {
1086 double sn
,ssn
,cs
,ccs
,s
,c
,w
[2],y
,y1
,y2
,e1
,e2
,xx
,cor
,res
;
1087 static const double t22
= 6291456.0;
1094 s
= y
*xx
*(sn3
+xx
*sn5
);
1095 c
= y
*dx
+xx
*(cs2
+xx
*(cs4
+ xx
*cs6
));
1106 cor
=(ccs
-cs
*c
-e1
*y2
-e2
*y
)-sn
*s
;
1108 cor
= cor
+((cs
-y
)-e1
*y1
);
1111 cor
= (cor
>0)? 1.0005*cor
+3.1e-30*ABS(orig
) : 1.0005*cor
-3.1e-30*ABS(orig
);
1112 if (res
== res
+ cor
) return (n
)?-res
:res
;
1114 __docos(ABS(x
),dx
,w
);
1115 cor
= (w
[1]>0)? 1.000000005*w
[1]+1.1e-30*ABS(orig
) : 1.000000005*w
[1]-1.1e-30*ABS(orig
);
1116 if (w
[0] == w
[0]+cor
) return (n
)?-w
[0]:w
[0];
1117 else return __mpcos1(orig
);
1121 weak_alias (__cos
, cos
)
1122 weak_alias (__sin
, sin
)
1124 #ifdef NO_LONG_DOUBLE
1125 strong_alias (__sin
, __sinl
)
1126 weak_alias (__sin
, sinl
)
1127 strong_alias (__cos
, __cosl
)
1128 weak_alias (__cos
, cosl
)