diff options
Diffstat (limited to 'sysdeps/ieee754/dbl-64/dosincos.c')
| -rw-r--r-- | sysdeps/ieee754/dbl-64/dosincos.c | 84 |
1 files changed, 45 insertions, 39 deletions
diff --git a/sysdeps/ieee754/dbl-64/dosincos.c b/sysdeps/ieee754/dbl-64/dosincos.c index 2a5f8fc..692fc85 100644 --- a/sysdeps/ieee754/dbl-64/dosincos.c +++ b/sysdeps/ieee754/dbl-64/dosincos.c @@ -5,9 +5,9 @@ * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by - * the Free Software Foundation; either version 2 of the License, or + * the Free Software Foundation; either version 2 of the License, or * (at your option) any later version. - * + * * This program is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the @@ -15,10 +15,10 @@ * * You should have received a copy of the GNU Lesser General Public License * along with this program; if not, write to the Free Software - * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */ /********************************************************************/ -/* */ +/* */ /* MODULE_NAME: dosincos.c */ /* */ /* */ @@ -34,7 +34,7 @@ #include "endian.h" -#include "mydefs.h" +#include "mydefs.h" #include "sincos.tbl" #include "dla.h" #include "dosincos.h" @@ -42,15 +42,18 @@ /* Routine receive Double-Length number (x+dx) and computing sin(x+dx) */ /* as Double-Length number and store it at array v .It computes it by */ /* arithmetic action on Double-Length numbers */ -/*(x+dx) between 0 and PI/4 */ +/*(x+dx) between 0 and PI/4 */ /***********************************************************************/ void dubsin(double x, double dx, double v[]) { - double xx,y,yy,z,zz,r,s,p,hx,tx,hy,ty,q,c,cc,d,dd,d2,dd2,e,ee, + double r,s,p,hx,tx,hy,ty,q,c,cc,d,dd,d2,dd2,e,ee, sn,ssn,cs,ccs,ds,dss,dc,dcc; +#if 0 + double xx,y,yy,z,zz; +#endif mynumber u; int4 k; - + u.x=x+big.x; k = u.i[LOW_HALF]<<2; x=x-(u.x-big.x); @@ -63,26 +66,26 @@ void dubsin(double x, double dx, double v[]) { cs=sincos.x[k+2]; /* */ ccs=sincos.x[k+3]; /* */ MUL2(d2,dd2,s7.x,ss7.x,ds,dss,p,hx,tx,hy,ty,q,c,cc); /* Taylor */ - ADD2(ds,dss,s5.x,ss5.x,ds,dss,r,s); + ADD2(ds,dss,s5.x,ss5.x,ds,dss,r,s); MUL2(d2,dd2,ds,dss,ds,dss,p,hx,tx,hy,ty,q,c,cc); /* series */ - ADD2(ds,dss,s3.x,ss3.x,ds,dss,r,s); + ADD2(ds,dss,s3.x,ss3.x,ds,dss,r,s); MUL2(d2,dd2,ds,dss,ds,dss,p,hx,tx,hy,ty,q,c,cc); /* for sin */ MUL2(d,dd,ds,dss,ds,dss,p,hx,tx,hy,ty,q,c,cc); ADD2(ds,dss,d,dd,ds,dss,r,s); /* ds=sin(t) */ - + MUL2(d2,dd2,c8.x,cc8.x,dc,dcc,p,hx,tx,hy,ty,q,c,cc); ;/* Taylor */ - ADD2(dc,dcc,c6.x,cc6.x,dc,dcc,r,s); + ADD2(dc,dcc,c6.x,cc6.x,dc,dcc,r,s); MUL2(d2,dd2,dc,dcc,dc,dcc,p,hx,tx,hy,ty,q,c,cc); /* series */ - ADD2(dc,dcc,c4.x,cc4.x,dc,dcc,r,s); + ADD2(dc,dcc,c4.x,cc4.x,dc,dcc,r,s); MUL2(d2,dd2,dc,dcc,dc,dcc,p,hx,tx,hy,ty,q,c,cc); /* for cos */ - ADD2(dc,dcc,c2.x,cc2.x,dc,dcc,r,s); + ADD2(dc,dcc,c2.x,cc2.x,dc,dcc,r,s); MUL2(d2,dd2,dc,dcc,dc,dcc,p,hx,tx,hy,ty,q,c,cc); /* dc=cos(t) */ - + MUL2(cs,ccs,ds,dss,e,ee,p,hx,tx,hy,ty,q,c,cc); MUL2(dc,dcc,sn,ssn,dc,dcc,p,hx,tx,hy,ty,q,c,cc); - SUB2(e,ee,dc,dcc,e,ee,r,s); + SUB2(e,ee,dc,dcc,e,ee,r,s); ADD2(e,ee,sn,ssn,e,ee,r,s); /* e+ee=sin(x+dx) */ - + v[0]=e; v[1]=ee; } @@ -94,46 +97,49 @@ void dubsin(double x, double dx, double v[]) { /**********************************************************************/ void dubcos(double x, double dx, double v[]) { - double xx,y,yy,z,zz,r,s,p,hx,tx,hy,ty,q,c,cc,d,dd,d2,dd2,e,ee, + double r,s,p,hx,tx,hy,ty,q,c,cc,d,dd,d2,dd2,e,ee, sn,ssn,cs,ccs,ds,dss,dc,dcc; +#if 0 + double xx,y,yy,z,zz; +#endif mynumber u; - int4 k; - u.x=x+big.x; + int4 k; + u.x=x+big.x; k = u.i[LOW_HALF]<<2; x=x-(u.x-big.x); - d=x+dx; + d=x+dx; dd=(x-d)+dx; /* cos(x+dx)=cos(Xi+t)=cos(Xi)cos(t) - sin(Xi)sin(t) */ - MUL2(d,dd,d,dd,d2,dd2,p,hx,tx,hy,ty,q,c,cc); + MUL2(d,dd,d,dd,d2,dd2,p,hx,tx,hy,ty,q,c,cc); sn=sincos.x[k]; /* */ ssn=sincos.x[k+1]; /* sin(Xi) and cos(Xi) */ cs=sincos.x[k+2]; /* */ ccs=sincos.x[k+3]; /* */ MUL2(d2,dd2,s7.x,ss7.x,ds,dss,p,hx,tx,hy,ty,q,c,cc); - ADD2(ds,dss,s5.x,ss5.x,ds,dss,r,s); + ADD2(ds,dss,s5.x,ss5.x,ds,dss,r,s); MUL2(d2,dd2,ds,dss,ds,dss,p,hx,tx,hy,ty,q,c,cc); - ADD2(ds,dss,s3.x,ss3.x,ds,dss,r,s); + ADD2(ds,dss,s3.x,ss3.x,ds,dss,r,s); MUL2(d2,dd2,ds,dss,ds,dss,p,hx,tx,hy,ty,q,c,cc); MUL2(d,dd,ds,dss,ds,dss,p,hx,tx,hy,ty,q,c,cc); - ADD2(ds,dss,d,dd,ds,dss,r,s); - + ADD2(ds,dss,d,dd,ds,dss,r,s); + MUL2(d2,dd2,c8.x,cc8.x,dc,dcc,p,hx,tx,hy,ty,q,c,cc); - ADD2(dc,dcc,c6.x,cc6.x,dc,dcc,r,s); + ADD2(dc,dcc,c6.x,cc6.x,dc,dcc,r,s); MUL2(d2,dd2,dc,dcc,dc,dcc,p,hx,tx,hy,ty,q,c,cc); - ADD2(dc,dcc,c4.x,cc4.x,dc,dcc,r,s); + ADD2(dc,dcc,c4.x,cc4.x,dc,dcc,r,s); MUL2(d2,dd2,dc,dcc,dc,dcc,p,hx,tx,hy,ty,q,c,cc); - ADD2(dc,dcc,c2.x,cc2.x,dc,dcc,r,s); + ADD2(dc,dcc,c2.x,cc2.x,dc,dcc,r,s); MUL2(d2,dd2,dc,dcc,dc,dcc,p,hx,tx,hy,ty,q,c,cc); - + MUL2(cs,ccs,ds,dss,e,ee,p,hx,tx,hy,ty,q,c,cc); MUL2(dc,dcc,sn,ssn,dc,dcc,p,hx,tx,hy,ty,q,c,cc); MUL2(d2,dd2,s7.x,ss7.x,ds,dss,p,hx,tx,hy,ty,q,c,cc); ADD2(ds,dss,s5.x,ss5.x,ds,dss,r,s); MUL2(d2,dd2,ds,dss,ds,dss,p,hx,tx,hy,ty,q,c,cc); - ADD2(ds,dss,s3.x,ss3.x,ds,dss,r,s); + ADD2(ds,dss,s3.x,ss3.x,ds,dss,r,s); MUL2(d2,dd2,ds,dss,ds,dss,p,hx,tx,hy,ty,q,c,cc); - MUL2(d,dd,ds,dss,ds,dss,p,hx,tx,hy,ty,q,c,cc); - ADD2(ds,dss,d,dd,ds,dss,r,s); + MUL2(d,dd,ds,dss,ds,dss,p,hx,tx,hy,ty,q,c,cc); + ADD2(ds,dss,d,dd,ds,dss,r,s); MUL2(d2,dd2,c8.x,cc8.x,dc,dcc,p,hx,tx,hy,ty,q,c,cc); ADD2(dc,dcc,c6.x,cc6.x,dc,dcc,r,s); MUL2(d2,dd2,dc,dcc,dc,dcc,p,hx,tx,hy,ty,q,c,cc); @@ -143,9 +149,9 @@ void dubcos(double x, double dx, double v[]) { MUL2(d2,dd2,dc,dcc,dc,dcc,p,hx,tx,hy,ty,q,c,cc); MUL2(sn,ssn,ds,dss,e,ee,p,hx,tx,hy,ty,q,c,cc); MUL2(dc,dcc,cs,ccs,dc,dcc,p,hx,tx,hy,ty,q,c,cc); - ADD2(e,ee,dc,dcc,e,ee,r,s); - SUB2(cs,ccs,e,ee,e,ee,r,s); - + ADD2(e,ee,dc,dcc,e,ee,r,s); + SUB2(cs,ccs,e,ee,e,ee,r,s); + v[0]=e; v[1]=ee; } @@ -164,14 +170,14 @@ void docos(double x, double dx, double v[]) { yy=hp1.x-yy; y=p+yy; yy=(p-y)+yy; - if (y>0) {dubsin(y,yy,w); v[0]=w[0]; v[1]=w[1];} + if (y>0) {dubsin(y,yy,w); v[0]=w[0]; v[1]=w[1];} /* cos(x) = sin ( 90 - x ) */ - else {dubsin(-y,-yy,w); v[0]=-w[0]; v[1]=-w[1]; + else {dubsin(-y,-yy,w); v[0]=-w[0]; v[1]=-w[1]; } } else { /* y>= 3/4 * PI */ p=2.0*hp0.x-y; /* p = PI- y */ - yy=2.0*hp1.x-yy; + yy=2.0*hp1.x-yy; y=p+yy; yy=(p-y)+yy; dubcos(y,yy,w); |