aboutsummaryrefslogtreecommitdiff
path: root/sysdeps/ieee754/dbl-64/dosincos.c
diff options
context:
space:
mode:
Diffstat (limited to 'sysdeps/ieee754/dbl-64/dosincos.c')
-rw-r--r--sysdeps/ieee754/dbl-64/dosincos.c84
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);
generated by cgit v1.1 at 2024-11-26 00:18:23 +0000