sysdeps/ieee754/dbl-64/dosincos.c

   1 /*
   2  * IBM Accurate Mathematical Library
   3  * Copyright (c) International Business Machines Corp., 2001
   4  *
   5  * This program is free software; you can redistribute it and/or modify
   6  * it under the terms of the GNU Lesser General Public License as published by
   7  * the Free Software Foundation; either version 2 of the License, or
   8  * (at your option) any later version.
   9  *
  10  * This program is distributed in the hope that it will be useful,
  11  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  13  * GNU General Public License for more details.
  14  *
  15  * You should have received a copy of the GNU Lesser General Public License
  16  * along with this program; if not, write to the Free Software
  17  * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
  18  */
  19 /********************************************************************/
  20 /*                                                                  */
  21 /* MODULE_NAME: dosincos.c                                          */
  22 /*                                                                  */
  23 /*                                                                  */
  24 /* FUNCTIONS:   dubsin                                              */
  25 /*              dubcos                                              */
  26 /*              docos                                               */
  27 /* FILES NEEDED: endian.h mydefs.h dla.h dosincos.h                 */
  28 /*               sincos.tbl                                         */
  29 /*                                                                  */
  30 /* Routines compute sin() and cos() as Double-Length numbers         */
  31 /********************************************************************/
  32
  33
  34
  35 #include "endian.h"
  36 #include "mydefs.h"
  37 #include "sincos.tbl"
  38 #include "dla.h"
  39 #include "dosincos.h"
  40 #include "math_private.h"
  41
  42 /***********************************************************************/
  43 /* Routine receive Double-Length number (x+dx) and computing sin(x+dx) */
  44 /* as Double-Length number and store it at array v .It computes it by  */
  45 /* arithmetic action on Double-Length numbers                          */
  46 /*(x+dx) between 0 and PI/4                                            */
  47 /***********************************************************************/
  48
  49 void __dubsin(double x, double dx, double v[]) {
  50   double r,s,p,hx,tx,hy,ty,q,c,cc,d,dd,d2,dd2,e,ee,
  51     sn,ssn,cs,ccs,ds,dss,dc,dcc;
  52 #if 0
  53   double xx,y,yy,z,zz;
  54 #endif
  55   mynumber u;
  56   int4 k;
  57
  58   u.x=x+big.x;
  59   k = u.i[LOW_HALF]<<2;
  60   x=x-(u.x-big.x);
  61   d=x+dx;
  62   dd=(x-d)+dx;
  63          /* sin(x+dx)=sin(Xi+t)=sin(Xi)*cos(t) + cos(Xi)sin(t) where t ->0 */
  64   MUL2(d,dd,d,dd,d2,dd2,p,hx,tx,hy,ty,q,c,cc);
  65   sn=sincos.x[k];     /*                                  */
  66   ssn=sincos.x[k+1];  /*      sin(Xi) and cos(Xi)         */
  67   cs=sincos.x[k+2];   /*                                  */
  68   ccs=sincos.x[k+3];  /*                                  */
  69   MUL2(d2,dd2,s7.x,ss7.x,ds,dss,p,hx,tx,hy,ty,q,c,cc);  /* Taylor    */
  70   ADD2(ds,dss,s5.x,ss5.x,ds,dss,r,s);
  71   MUL2(d2,dd2,ds,dss,ds,dss,p,hx,tx,hy,ty,q,c,cc);      /* series    */
  72   ADD2(ds,dss,s3.x,ss3.x,ds,dss,r,s);
  73   MUL2(d2,dd2,ds,dss,ds,dss,p,hx,tx,hy,ty,q,c,cc);      /* for sin   */
  74   MUL2(d,dd,ds,dss,ds,dss,p,hx,tx,hy,ty,q,c,cc);
  75   ADD2(ds,dss,d,dd,ds,dss,r,s);                         /* ds=sin(t) */
  76
  77   MUL2(d2,dd2,c8.x,cc8.x,dc,dcc,p,hx,tx,hy,ty,q,c,cc); ;/* Taylor    */
  78   ADD2(dc,dcc,c6.x,cc6.x,dc,dcc,r,s);
  79   MUL2(d2,dd2,dc,dcc,dc,dcc,p,hx,tx,hy,ty,q,c,cc);      /* series    */
  80   ADD2(dc,dcc,c4.x,cc4.x,dc,dcc,r,s);
  81   MUL2(d2,dd2,dc,dcc,dc,dcc,p,hx,tx,hy,ty,q,c,cc);      /* for cos   */
  82   ADD2(dc,dcc,c2.x,cc2.x,dc,dcc,r,s);
  83   MUL2(d2,dd2,dc,dcc,dc,dcc,p,hx,tx,hy,ty,q,c,cc);      /* dc=cos(t) */
  84
  85   MUL2(cs,ccs,ds,dss,e,ee,p,hx,tx,hy,ty,q,c,cc);
  86   MUL2(dc,dcc,sn,ssn,dc,dcc,p,hx,tx,hy,ty,q,c,cc);
  87   SUB2(e,ee,dc,dcc,e,ee,r,s);
  88   ADD2(e,ee,sn,ssn,e,ee,r,s);                    /* e+ee=sin(x+dx) */
  89
  90   v[0]=e;
  91   v[1]=ee;
  92 }
  93 /**********************************************************************/
  94 /* Routine receive Double-Length number (x+dx) and computes cos(x+dx) */
  95 /* as Double-Length number and store it in array v .It computes it by */
  96 /* arithmetic action on Double-Length numbers                         */
  97 /*(x+dx) between 0 and PI/4                                           */
  98 /**********************************************************************/
  99
 100 void __dubcos(double x, double dx, double v[]) {
 101   double r,s,p,hx,tx,hy,ty,q,c,cc,d,dd,d2,dd2,e,ee,
 102     sn,ssn,cs,ccs,ds,dss,dc,dcc;
 103 #if 0
 104   double xx,y,yy,z,zz;
 105 #endif
 106   mynumber u;
 107   int4 k;
 108   u.x=x+big.x;
 109   k = u.i[LOW_HALF]<<2;
 110   x=x-(u.x-big.x);
 111   d=x+dx;
 112   dd=(x-d)+dx;  /* cos(x+dx)=cos(Xi+t)=cos(Xi)cos(t) - sin(Xi)sin(t) */
 113   MUL2(d,dd,d,dd,d2,dd2,p,hx,tx,hy,ty,q,c,cc);
 114   sn=sincos.x[k];     /*                                  */
 115   ssn=sincos.x[k+1];  /*      sin(Xi) and cos(Xi)         */
 116   cs=sincos.x[k+2];   /*                                  */
 117   ccs=sincos.x[k+3];  /*                                  */
 118   MUL2(d2,dd2,s7.x,ss7.x,ds,dss,p,hx,tx,hy,ty,q,c,cc);
 119   ADD2(ds,dss,s5.x,ss5.x,ds,dss,r,s);
 120   MUL2(d2,dd2,ds,dss,ds,dss,p,hx,tx,hy,ty,q,c,cc);
 121   ADD2(ds,dss,s3.x,ss3.x,ds,dss,r,s);
 122   MUL2(d2,dd2,ds,dss,ds,dss,p,hx,tx,hy,ty,q,c,cc);
 123   MUL2(d,dd,ds,dss,ds,dss,p,hx,tx,hy,ty,q,c,cc);
 124   ADD2(ds,dss,d,dd,ds,dss,r,s);
 125
 126   MUL2(d2,dd2,c8.x,cc8.x,dc,dcc,p,hx,tx,hy,ty,q,c,cc);
 127   ADD2(dc,dcc,c6.x,cc6.x,dc,dcc,r,s);
 128   MUL2(d2,dd2,dc,dcc,dc,dcc,p,hx,tx,hy,ty,q,c,cc);
 129   ADD2(dc,dcc,c4.x,cc4.x,dc,dcc,r,s);
 130   MUL2(d2,dd2,dc,dcc,dc,dcc,p,hx,tx,hy,ty,q,c,cc);
 131   ADD2(dc,dcc,c2.x,cc2.x,dc,dcc,r,s);
 132   MUL2(d2,dd2,dc,dcc,dc,dcc,p,hx,tx,hy,ty,q,c,cc);
 133
 134   MUL2(cs,ccs,ds,dss,e,ee,p,hx,tx,hy,ty,q,c,cc);
 135   MUL2(dc,dcc,sn,ssn,dc,dcc,p,hx,tx,hy,ty,q,c,cc);
 136
 137   MUL2(d2,dd2,s7.x,ss7.x,ds,dss,p,hx,tx,hy,ty,q,c,cc);
 138   ADD2(ds,dss,s5.x,ss5.x,ds,dss,r,s);
 139   MUL2(d2,dd2,ds,dss,ds,dss,p,hx,tx,hy,ty,q,c,cc);
 140   ADD2(ds,dss,s3.x,ss3.x,ds,dss,r,s);
 141   MUL2(d2,dd2,ds,dss,ds,dss,p,hx,tx,hy,ty,q,c,cc);
 142   MUL2(d,dd,ds,dss,ds,dss,p,hx,tx,hy,ty,q,c,cc);
 143   ADD2(ds,dss,d,dd,ds,dss,r,s);
 144   MUL2(d2,dd2,c8.x,cc8.x,dc,dcc,p,hx,tx,hy,ty,q,c,cc);
 145   ADD2(dc,dcc,c6.x,cc6.x,dc,dcc,r,s);
 146   MUL2(d2,dd2,dc,dcc,dc,dcc,p,hx,tx,hy,ty,q,c,cc);
 147   ADD2(dc,dcc,c4.x,cc4.x,dc,dcc,r,s);
 148   MUL2(d2,dd2,dc,dcc,dc,dcc,p,hx,tx,hy,ty,q,c,cc);
 149   ADD2(dc,dcc,c2.x,cc2.x,dc,dcc,r,s);
 150   MUL2(d2,dd2,dc,dcc,dc,dcc,p,hx,tx,hy,ty,q,c,cc);
 151   MUL2(sn,ssn,ds,dss,e,ee,p,hx,tx,hy,ty,q,c,cc);
 152   MUL2(dc,dcc,cs,ccs,dc,dcc,p,hx,tx,hy,ty,q,c,cc);
 153   ADD2(e,ee,dc,dcc,e,ee,r,s);
 154   SUB2(cs,ccs,e,ee,e,ee,r,s);
 155
 156   v[0]=e;
 157   v[1]=ee;
 158 }
 159 /**********************************************************************/
 160 /* Routine receive Double-Length number (x+dx) and computes cos(x+dx) */
 161 /* as Double-Length number and store it in array v                    */
 162 /**********************************************************************/
 163 void __docos(double x, double dx, double v[]) {
 164   double y,yy,p,w[2];
 165   if (x>0) {y=x; yy=dx;}
 166      else {y=-x; yy=-dx;}
 167   if (y<0.5*hp0.x)                                 /*  y< PI/4    */
 168            {__dubcos(y,yy,w); v[0]=w[0]; v[1]=w[1];}
 169      else if (y<1.5*hp0.x) {                       /* y< 3/4 * PI */
 170        p=hp0.x-y;  /* p = PI/2 - y */
 171        yy=hp1.x-yy;
 172        y=p+yy;
 173        yy=(p-y)+yy;
 174        if (y>0) {__dubsin(y,yy,w); v[0]=w[0]; v[1]=w[1];}
 175                                        /* cos(x) = sin ( 90 -  x ) */
 176          else {__dubsin(-y,-yy,w); v[0]=-w[0]; v[1]=-w[1];
 177          }
 178      }
 179   else { /* y>= 3/4 * PI */
 180     p=2.0*hp0.x-y;    /* p = PI- y */
 181     yy=2.0*hp1.x-yy;
 182     y=p+yy;
 183     yy=(p-y)+yy;
 184     __dubcos(y,yy,w);
 185     v[0]=-w[0];
 186     v[1]=-w[1];
 187   }
 188 }