Back to index

glibc  2.9
dosincos.c
Go to the documentation of this file.
00001 /*
00002  * IBM Accurate Mathematical Library
00003  * written by International Business Machines Corp.
00004  * Copyright (C) 2001 Free Software Foundation
00005  *
00006  * This program is free software; you can redistribute it and/or modify
00007  * it under the terms of the GNU Lesser General Public License as published by
00008  * the Free Software Foundation; either version 2.1 of the License, or
00009  * (at your option) any later version.
00010  *
00011  * This program is distributed in the hope that it will be useful,
00012  * but WITHOUT ANY WARRANTY; without even the implied warranty of
00013  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00014  * GNU Lesser General Public License for more details.
00015  *
00016  * You should have received a copy of the GNU Lesser General Public License
00017  * along with this program; if not, write to the Free Software
00018  * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
00019  */
00020 /********************************************************************/
00021 /*                                                                  */
00022 /* MODULE_NAME: dosincos.c                                          */
00023 /*                                                                  */
00024 /*                                                                  */
00025 /* FUNCTIONS:   dubsin                                              */
00026 /*              dubcos                                              */
00027 /*              docos                                               */
00028 /* FILES NEEDED: endian.h mydefs.h dla.h dosincos.h                 */
00029 /*               sincos.tbl                                         */
00030 /*                                                                  */
00031 /* Routines compute sin() and cos() as Double-Length numbers         */
00032 /********************************************************************/
00033 
00034 
00035 
00036 #include "endian.h"
00037 #include "mydefs.h"
00038 #include "sincos.tbl"
00039 #include "dla.h"
00040 #include "dosincos.h"
00041 #include "math_private.h"
00042 
00043 /***********************************************************************/
00044 /* Routine receive Double-Length number (x+dx) and computing sin(x+dx) */
00045 /* as Double-Length number and store it at array v .It computes it by  */
00046 /* arithmetic action on Double-Length numbers                          */
00047 /*(x+dx) between 0 and PI/4                                            */
00048 /***********************************************************************/
00049 
00050 void __dubsin(double x, double dx, double v[]) {
00051   double r,s,p,hx,tx,hy,ty,q,c,cc,d,dd,d2,dd2,e,ee,
00052     sn,ssn,cs,ccs,ds,dss,dc,dcc;
00053 #if 0
00054   double xx,y,yy,z,zz;
00055 #endif
00056   mynumber u;
00057   int4 k;
00058 
00059   u.x=x+big.x;
00060   k = u.i[LOW_HALF]<<2;
00061   x=x-(u.x-big.x);
00062   d=x+dx;
00063   dd=(x-d)+dx;
00064          /* sin(x+dx)=sin(Xi+t)=sin(Xi)*cos(t) + cos(Xi)sin(t) where t ->0 */
00065   MUL2(d,dd,d,dd,d2,dd2,p,hx,tx,hy,ty,q,c,cc);
00066   sn=sincos.x[k];     /*                                  */
00067   ssn=sincos.x[k+1];  /*      sin(Xi) and cos(Xi)         */
00068   cs=sincos.x[k+2];   /*                                  */
00069   ccs=sincos.x[k+3];  /*                                  */
00070   MUL2(d2,dd2,s7.x,ss7.x,ds,dss,p,hx,tx,hy,ty,q,c,cc);  /* Taylor    */
00071   ADD2(ds,dss,s5.x,ss5.x,ds,dss,r,s);
00072   MUL2(d2,dd2,ds,dss,ds,dss,p,hx,tx,hy,ty,q,c,cc);      /* series    */
00073   ADD2(ds,dss,s3.x,ss3.x,ds,dss,r,s);
00074   MUL2(d2,dd2,ds,dss,ds,dss,p,hx,tx,hy,ty,q,c,cc);      /* for sin   */
00075   MUL2(d,dd,ds,dss,ds,dss,p,hx,tx,hy,ty,q,c,cc);
00076   ADD2(ds,dss,d,dd,ds,dss,r,s);                         /* ds=sin(t) */
00077 
00078   MUL2(d2,dd2,c8.x,cc8.x,dc,dcc,p,hx,tx,hy,ty,q,c,cc); ;/* Taylor    */
00079   ADD2(dc,dcc,c6.x,cc6.x,dc,dcc,r,s);
00080   MUL2(d2,dd2,dc,dcc,dc,dcc,p,hx,tx,hy,ty,q,c,cc);      /* series    */
00081   ADD2(dc,dcc,c4.x,cc4.x,dc,dcc,r,s);
00082   MUL2(d2,dd2,dc,dcc,dc,dcc,p,hx,tx,hy,ty,q,c,cc);      /* for cos   */
00083   ADD2(dc,dcc,c2.x,cc2.x,dc,dcc,r,s);
00084   MUL2(d2,dd2,dc,dcc,dc,dcc,p,hx,tx,hy,ty,q,c,cc);      /* dc=cos(t) */
00085 
00086   MUL2(cs,ccs,ds,dss,e,ee,p,hx,tx,hy,ty,q,c,cc);
00087   MUL2(dc,dcc,sn,ssn,dc,dcc,p,hx,tx,hy,ty,q,c,cc);
00088   SUB2(e,ee,dc,dcc,e,ee,r,s);
00089   ADD2(e,ee,sn,ssn,e,ee,r,s);                    /* e+ee=sin(x+dx) */
00090 
00091   v[0]=e;
00092   v[1]=ee;
00093 }
00094 /**********************************************************************/
00095 /* Routine receive Double-Length number (x+dx) and computes cos(x+dx) */
00096 /* as Double-Length number and store it in array v .It computes it by */
00097 /* arithmetic action on Double-Length numbers                         */
00098 /*(x+dx) between 0 and PI/4                                           */
00099 /**********************************************************************/
00100 
00101 void __dubcos(double x, double dx, double v[]) {
00102   double r,s,p,hx,tx,hy,ty,q,c,cc,d,dd,d2,dd2,e,ee,
00103     sn,ssn,cs,ccs,ds,dss,dc,dcc;
00104 #if 0
00105   double xx,y,yy,z,zz;
00106 #endif
00107   mynumber u;
00108   int4 k;
00109   u.x=x+big.x;
00110   k = u.i[LOW_HALF]<<2;
00111   x=x-(u.x-big.x);
00112   d=x+dx;
00113   dd=(x-d)+dx;  /* cos(x+dx)=cos(Xi+t)=cos(Xi)cos(t) - sin(Xi)sin(t) */
00114   MUL2(d,dd,d,dd,d2,dd2,p,hx,tx,hy,ty,q,c,cc);
00115   sn=sincos.x[k];     /*                                  */
00116   ssn=sincos.x[k+1];  /*      sin(Xi) and cos(Xi)         */
00117   cs=sincos.x[k+2];   /*                                  */
00118   ccs=sincos.x[k+3];  /*                                  */
00119   MUL2(d2,dd2,s7.x,ss7.x,ds,dss,p,hx,tx,hy,ty,q,c,cc);
00120   ADD2(ds,dss,s5.x,ss5.x,ds,dss,r,s);
00121   MUL2(d2,dd2,ds,dss,ds,dss,p,hx,tx,hy,ty,q,c,cc);
00122   ADD2(ds,dss,s3.x,ss3.x,ds,dss,r,s);
00123   MUL2(d2,dd2,ds,dss,ds,dss,p,hx,tx,hy,ty,q,c,cc);
00124   MUL2(d,dd,ds,dss,ds,dss,p,hx,tx,hy,ty,q,c,cc);
00125   ADD2(ds,dss,d,dd,ds,dss,r,s);
00126 
00127   MUL2(d2,dd2,c8.x,cc8.x,dc,dcc,p,hx,tx,hy,ty,q,c,cc);
00128   ADD2(dc,dcc,c6.x,cc6.x,dc,dcc,r,s);
00129   MUL2(d2,dd2,dc,dcc,dc,dcc,p,hx,tx,hy,ty,q,c,cc);
00130   ADD2(dc,dcc,c4.x,cc4.x,dc,dcc,r,s);
00131   MUL2(d2,dd2,dc,dcc,dc,dcc,p,hx,tx,hy,ty,q,c,cc);
00132   ADD2(dc,dcc,c2.x,cc2.x,dc,dcc,r,s);
00133   MUL2(d2,dd2,dc,dcc,dc,dcc,p,hx,tx,hy,ty,q,c,cc);
00134 
00135   MUL2(cs,ccs,ds,dss,e,ee,p,hx,tx,hy,ty,q,c,cc);
00136   MUL2(dc,dcc,sn,ssn,dc,dcc,p,hx,tx,hy,ty,q,c,cc);
00137 
00138   MUL2(d2,dd2,s7.x,ss7.x,ds,dss,p,hx,tx,hy,ty,q,c,cc);
00139   ADD2(ds,dss,s5.x,ss5.x,ds,dss,r,s);
00140   MUL2(d2,dd2,ds,dss,ds,dss,p,hx,tx,hy,ty,q,c,cc);
00141   ADD2(ds,dss,s3.x,ss3.x,ds,dss,r,s);
00142   MUL2(d2,dd2,ds,dss,ds,dss,p,hx,tx,hy,ty,q,c,cc);
00143   MUL2(d,dd,ds,dss,ds,dss,p,hx,tx,hy,ty,q,c,cc);
00144   ADD2(ds,dss,d,dd,ds,dss,r,s);
00145   MUL2(d2,dd2,c8.x,cc8.x,dc,dcc,p,hx,tx,hy,ty,q,c,cc);
00146   ADD2(dc,dcc,c6.x,cc6.x,dc,dcc,r,s);
00147   MUL2(d2,dd2,dc,dcc,dc,dcc,p,hx,tx,hy,ty,q,c,cc);
00148   ADD2(dc,dcc,c4.x,cc4.x,dc,dcc,r,s);
00149   MUL2(d2,dd2,dc,dcc,dc,dcc,p,hx,tx,hy,ty,q,c,cc);
00150   ADD2(dc,dcc,c2.x,cc2.x,dc,dcc,r,s);
00151   MUL2(d2,dd2,dc,dcc,dc,dcc,p,hx,tx,hy,ty,q,c,cc);
00152   MUL2(sn,ssn,ds,dss,e,ee,p,hx,tx,hy,ty,q,c,cc);
00153   MUL2(dc,dcc,cs,ccs,dc,dcc,p,hx,tx,hy,ty,q,c,cc);
00154   ADD2(e,ee,dc,dcc,e,ee,r,s);
00155   SUB2(cs,ccs,e,ee,e,ee,r,s);
00156 
00157   v[0]=e;
00158   v[1]=ee;
00159 }
00160 /**********************************************************************/
00161 /* Routine receive Double-Length number (x+dx) and computes cos(x+dx) */
00162 /* as Double-Length number and store it in array v                    */
00163 /**********************************************************************/
00164 void __docos(double x, double dx, double v[]) {
00165   double y,yy,p,w[2];
00166   if (x>0) {y=x; yy=dx;}
00167      else {y=-x; yy=-dx;}
00168   if (y<0.5*hp0.x)                                 /*  y< PI/4    */
00169            {__dubcos(y,yy,w); v[0]=w[0]; v[1]=w[1];}
00170      else if (y<1.5*hp0.x) {                       /* y< 3/4 * PI */
00171        p=hp0.x-y;  /* p = PI/2 - y */
00172        yy=hp1.x-yy;
00173        y=p+yy;
00174        yy=(p-y)+yy;
00175        if (y>0) {__dubsin(y,yy,w); v[0]=w[0]; v[1]=w[1];}
00176                                        /* cos(x) = sin ( 90 -  x ) */
00177          else {__dubsin(-y,-yy,w); v[0]=-w[0]; v[1]=-w[1];
00178         }
00179      }
00180   else { /* y>= 3/4 * PI */
00181     p=2.0*hp0.x-y;    /* p = PI- y */
00182     yy=2.0*hp1.x-yy;
00183     y=p+yy;
00184     yy=(p-y)+yy;
00185     __dubcos(y,yy,w);
00186     v[0]=-w[0];
00187     v[1]=-w[1];
00188   }
00189 }