sysdeps/ieee754/dbl-64/dosincos.c

   1 /*
   2  * IBM Accurate Mathematical Library
   3  * written by International Business Machines Corp.
   4  * Copyright (C) 2001, 2011 Free Software Foundation
   5  *
   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.1 of the License, or
   9  * (at your option) any later version.
  10  *
  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 Lesser General Public License for more details.
  15  *
  16  * You should have received a copy of the GNU Lesser General Public License
  17  * along with this program; if not, see <http://www.gnu.org/licenses/>.
  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 <dla.h>
  38 #include "dosincos.h"
  39 #include "math_private.h"
  40
  41 #ifndef SECTION
  42 # define SECTION
  43 #endif
  44
  45 extern const union
  46 {
  47   int4 i[880];
  48   double x[440];
  49 } __sincostab attribute_hidden;
  50
  51 /***********************************************************************/
  52 /* Routine receive Double-Length number (x+dx) and computing sin(x+dx) */
  53 /* as Double-Length number and store it at array v .It computes it by  */
  54 /* arithmetic action on Double-Length numbers                          */
  55 /*(x+dx) between 0 and PI/4                                            */
  56 /***********************************************************************/
  57
  58 void
  59 SECTION
  60 __dubsin(double x, double dx, double v[]) {
  61   double r,s,c,cc,d,dd,d2,dd2,e,ee,
  62     sn,ssn,cs,ccs,ds,dss,dc,dcc;
  63 #ifndef DLA_FMS
  64   double p,hx,tx,hy,ty,q;
  65 #endif
  66 #if 0
  67   double xx,y,yy,z,zz;
  68 #endif
  69   mynumber u;
  70   int4 k;
  71
  72   u.x=x+big.x;
  73   k = u.i[LOW_HALF]<<2;
  74   x=x-(u.x-big.x);
  75   d=x+dx;
  76   dd=(x-d)+dx;
  77          /* sin(x+dx)=sin(Xi+t)=sin(Xi)*cos(t) + cos(Xi)sin(t) where t ->0 */
  78   MUL2(d,dd,d,dd,d2,dd2,p,hx,tx,hy,ty,q,c,cc);
  79   sn=__sincostab.x[k];     /*                                  */
  80   ssn=__sincostab.x[k+1];  /*      sin(Xi) and cos(Xi)         */
  81   cs=__sincostab.x[k+2];   /*                                  */
  82   ccs=__sincostab.x[k+3];  /*                                  */
  83   MUL2(d2,dd2,s7.x,ss7.x,ds,dss,p,hx,tx,hy,ty,q,c,cc);  /* Taylor    */
  84   ADD2(ds,dss,s5.x,ss5.x,ds,dss,r,s);
  85   MUL2(d2,dd2,ds,dss,ds,dss,p,hx,tx,hy,ty,q,c,cc);      /* series    */
  86   ADD2(ds,dss,s3.x,ss3.x,ds,dss,r,s);
  87   MUL2(d2,dd2,ds,dss,ds,dss,p,hx,tx,hy,ty,q,c,cc);      /* for sin   */
  88   MUL2(d,dd,ds,dss,ds,dss,p,hx,tx,hy,ty,q,c,cc);
  89   ADD2(ds,dss,d,dd,ds,dss,r,s);                         /* ds=sin(t) */
  90
  91   MUL2(d2,dd2,c8.x,cc8.x,dc,dcc,p,hx,tx,hy,ty,q,c,cc); ;/* Taylor    */
  92   ADD2(dc,dcc,c6.x,cc6.x,dc,dcc,r,s);
  93   MUL2(d2,dd2,dc,dcc,dc,dcc,p,hx,tx,hy,ty,q,c,cc);      /* series    */
  94   ADD2(dc,dcc,c4.x,cc4.x,dc,dcc,r,s);
  95   MUL2(d2,dd2,dc,dcc,dc,dcc,p,hx,tx,hy,ty,q,c,cc);      /* for cos   */
  96   ADD2(dc,dcc,c2.x,cc2.x,dc,dcc,r,s);
  97   MUL2(d2,dd2,dc,dcc,dc,dcc,p,hx,tx,hy,ty,q,c,cc);      /* dc=cos(t) */
  98
  99   MUL2(cs,ccs,ds,dss,e,ee,p,hx,tx,hy,ty,q,c,cc);
 100   MUL2(dc,dcc,sn,ssn,dc,dcc,p,hx,tx,hy,ty,q,c,cc);
 101   SUB2(e,ee,dc,dcc,e,ee,r,s);
 102   ADD2(e,ee,sn,ssn,e,ee,r,s);                    /* e+ee=sin(x+dx) */
 103
 104   v[0]=e;
 105   v[1]=ee;
 106 }
 107 /**********************************************************************/
 108 /* Routine receive Double-Length number (x+dx) and computes cos(x+dx) */
 109 /* as Double-Length number and store it in array v .It computes it by */
 110 /* arithmetic action on Double-Length numbers                         */
 111 /*(x+dx) between 0 and PI/4                                           */
 112 /**********************************************************************/
 113
 114 void
 115 SECTION
 116 __dubcos(double x, double dx, double v[]) {
 117   double r,s,c,cc,d,dd,d2,dd2,e,ee,
 118     sn,ssn,cs,ccs,ds,dss,dc,dcc;
 119 #ifndef DLA_FMS
 120   double p,hx,tx,hy,ty,q;
 121 #endif
 122 #if 0
 123   double xx,y,yy,z,zz;
 124 #endif
 125   mynumber u;
 126   int4 k;
 127   u.x=x+big.x;
 128   k = u.i[LOW_HALF]<<2;
 129   x=x-(u.x-big.x);
 130   d=x+dx;
 131   dd=(x-d)+dx;  /* cos(x+dx)=cos(Xi+t)=cos(Xi)cos(t) - sin(Xi)sin(t) */
 132   MUL2(d,dd,d,dd,d2,dd2,p,hx,tx,hy,ty,q,c,cc);
 133   sn=__sincostab.x[k];     /*                                  */
 134   ssn=__sincostab.x[k+1];  /*      sin(Xi) and cos(Xi)         */
 135   cs=__sincostab.x[k+2];   /*                                  */
 136   ccs=__sincostab.x[k+3];  /*                                  */
 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
 145   MUL2(d2,dd2,c8.x,cc8.x,dc,dcc,p,hx,tx,hy,ty,q,c,cc);
 146   ADD2(dc,dcc,c6.x,cc6.x,dc,dcc,r,s);
 147   MUL2(d2,dd2,dc,dcc,dc,dcc,p,hx,tx,hy,ty,q,c,cc);
 148   ADD2(dc,dcc,c4.x,cc4.x,dc,dcc,r,s);
 149   MUL2(d2,dd2,dc,dcc,dc,dcc,p,hx,tx,hy,ty,q,c,cc);
 150   ADD2(dc,dcc,c2.x,cc2.x,dc,dcc,r,s);
 151   MUL2(d2,dd2,dc,dcc,dc,dcc,p,hx,tx,hy,ty,q,c,cc);
 152
 153   MUL2(cs,ccs,ds,dss,e,ee,p,hx,tx,hy,ty,q,c,cc);
 154   MUL2(dc,dcc,sn,ssn,dc,dcc,p,hx,tx,hy,ty,q,c,cc);
 155
 156   MUL2(d2,dd2,s7.x,ss7.x,ds,dss,p,hx,tx,hy,ty,q,c,cc);
 157   ADD2(ds,dss,s5.x,ss5.x,ds,dss,r,s);
 158   MUL2(d2,dd2,ds,dss,ds,dss,p,hx,tx,hy,ty,q,c,cc);
 159   ADD2(ds,dss,s3.x,ss3.x,ds,dss,r,s);
 160   MUL2(d2,dd2,ds,dss,ds,dss,p,hx,tx,hy,ty,q,c,cc);
 161   MUL2(d,dd,ds,dss,ds,dss,p,hx,tx,hy,ty,q,c,cc);
 162   ADD2(ds,dss,d,dd,ds,dss,r,s);
 163   MUL2(d2,dd2,c8.x,cc8.x,dc,dcc,p,hx,tx,hy,ty,q,c,cc);
 164   ADD2(dc,dcc,c6.x,cc6.x,dc,dcc,r,s);
 165   MUL2(d2,dd2,dc,dcc,dc,dcc,p,hx,tx,hy,ty,q,c,cc);
 166   ADD2(dc,dcc,c4.x,cc4.x,dc,dcc,r,s);
 167   MUL2(d2,dd2,dc,dcc,dc,dcc,p,hx,tx,hy,ty,q,c,cc);
 168   ADD2(dc,dcc,c2.x,cc2.x,dc,dcc,r,s);
 169   MUL2(d2,dd2,dc,dcc,dc,dcc,p,hx,tx,hy,ty,q,c,cc);
 170   MUL2(sn,ssn,ds,dss,e,ee,p,hx,tx,hy,ty,q,c,cc);
 171   MUL2(dc,dcc,cs,ccs,dc,dcc,p,hx,tx,hy,ty,q,c,cc);
 172   ADD2(e,ee,dc,dcc,e,ee,r,s);
 173   SUB2(cs,ccs,e,ee,e,ee,r,s);
 174
 175   v[0]=e;
 176   v[1]=ee;
 177 }
 178 /**********************************************************************/
 179 /* Routine receive Double-Length number (x+dx) and computes cos(x+dx) */
 180 /* as Double-Length number and store it in array v                    */
 181 /**********************************************************************/
 182 void
 183 SECTION
 184 __docos(double x, double dx, double v[]) {
 185   double y,yy,p,w[2];
 186   if (x>0) {y=x; yy=dx;}
 187      else {y=-x; yy=-dx;}
 188   if (y<0.5*hp0.x)                                 /*  y< PI/4    */
 189            {__dubcos(y,yy,w); v[0]=w[0]; v[1]=w[1];}
 190      else if (y<1.5*hp0.x) {                       /* y< 3/4 * PI */
 191        p=hp0.x-y;  /* p = PI/2 - y */
 192        yy=hp1.x-yy;
 193        y=p+yy;
 194        yy=(p-y)+yy;
 195        if (y>0) {__dubsin(y,yy,w); v[0]=w[0]; v[1]=w[1];}
 196                                        /* cos(x) = sin ( 90 -  x ) */
 197          else {__dubsin(-y,-yy,w); v[0]=-w[0]; v[1]=-w[1];
 198          }
 199      }
 200   else { /* y>= 3/4 * PI */
 201     p=2.0*hp0.x-y;    /* p = PI- y */
 202     yy=2.0*hp1.x-yy;
 203     y=p+yy;
 204     yy=(p-y)+yy;
 205     __dubcos(y,yy,w);
 206     v[0]=-w[0];
 207     v[1]=-w[1];
 208   }
 209 }