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 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 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)