sysdeps/ieee754/dbl-64/dosincos.c

   1 /*
   2  * IBM Accurate Mathematical Library
   3  * written by International Business Machines Corp.
   4  * Copyright (C) 2001-2014 Free Software Foundation, Inc.
   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 {
  62   double r, s, c, cc, d, dd, d2, dd2, e, ee,
  63          sn, ssn, cs, ccs, ds, dss, dc, dcc;
  64 #ifndef DLA_FMS
  65   double p, hx, tx, hy, ty, q;
  66 #endif
  67   mynumber u;
  68   int4 k;
  69
  70   u.x = x + big.x;
  71   k = u.i[LOW_HALF] << 2;
  72   x = x - (u.x - big.x);
  73   d = x + dx;
  74   dd = (x - d) + dx;
  75   /* sin(x+dx)=sin(Xi+t)=sin(Xi)*cos(t) + cos(Xi)sin(t) where t ->0 */
  76   MUL2 (d, dd, d, dd, d2, dd2, p, hx, tx, hy, ty, q, c, cc);
  77   sn = __sincostab.x[k];       /*                                  */
  78   ssn = __sincostab.x[k + 1];  /*      sin(Xi) and cos(Xi)         */
  79   cs = __sincostab.x[k + 2];   /*                                  */
  80   ccs = __sincostab.x[k + 3];  /*                                  */
  81   /* Taylor series for sin ds=sin(t) */
  82   MUL2 (d2, dd2, s7.x, ss7.x, ds, dss, p, hx, tx, hy, ty, q, c, cc);
  83   ADD2 (ds, dss, s5.x, ss5.x, ds, dss, r, s);
  84   MUL2 (d2, dd2, ds, dss, ds, dss, p, hx, tx, hy, ty, q, c, cc);
  85   ADD2 (ds, dss, s3.x, ss3.x, ds, dss, r, s);
  86   MUL2 (d2, dd2, ds, dss, ds, dss, p, hx, tx, hy, ty, q, c, cc);
  87   MUL2 (d, dd, ds, dss, ds, dss, p, hx, tx, hy, ty, q, c, cc);
  88   ADD2 (ds, dss, d, dd, ds, dss, r, s);
  89
  90   /* Taylor series for cos dc=cos(t) */
  91   MUL2 (d2, dd2, c8.x, cc8.x, dc, dcc, p, hx, tx, hy, ty, q, c, cc);
  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);
  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);
  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);
  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 {
 118   double r, s, c, cc, d, dd, d2, dd2, e, ee,
 119          sn, ssn, cs, ccs, ds, dss, dc, dcc;
 120 #ifndef DLA_FMS
 121   double p, hx, tx, hy, ty, q;
 122 #endif
 123   mynumber u;
 124   int4 k;
 125   u.x = x + big.x;
 126   k = u.i[LOW_HALF] << 2;
 127   x = x - (u.x - big.x);
 128   d = x + dx;
 129   dd = (x - d) + dx;  /* cos(x+dx)=cos(Xi+t)=cos(Xi)cos(t) - sin(Xi)sin(t) */
 130   MUL2 (d, dd, d, dd, d2, dd2, p, hx, tx, hy, ty, q, c, cc);
 131   sn = __sincostab.x[k];     /*                                  */
 132   ssn = __sincostab.x[k + 1];  /*      sin(Xi) and cos(Xi)         */
 133   cs = __sincostab.x[k + 2];   /*                                  */
 134   ccs = __sincostab.x[k + 3];  /*                                  */
 135   MUL2 (d2, dd2, s7.x, ss7.x, ds, dss, p, hx, tx, hy, ty, q, c, cc);
 136   ADD2 (ds, dss, s5.x, ss5.x, ds, dss, r, s);
 137   MUL2 (d2, dd2, ds, dss, ds, dss, p, hx, tx, hy, ty, q, c, cc);
 138   ADD2 (ds, dss, s3.x, ss3.x, ds, dss, r, s);
 139   MUL2 (d2, dd2, ds, dss, ds, dss, p, hx, tx, hy, ty, q, c, cc);
 140   MUL2 (d, dd, ds, dss, ds, dss, p, hx, tx, hy, ty, q, c, cc);
 141   ADD2 (ds, dss, d, dd, ds, dss, r, s);
 142
 143   MUL2 (d2, dd2, c8.x, cc8.x, dc, dcc, p, hx, tx, hy, ty, q, c, cc);
 144   ADD2 (dc, dcc, c6.x, cc6.x, dc, dcc, r, s);
 145   MUL2 (d2, dd2, dc, dcc, dc, dcc, p, hx, tx, hy, ty, q, c, cc);
 146   ADD2 (dc, dcc, c4.x, cc4.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, c2.x, cc2.x, dc, dcc, r, s);
 149   MUL2 (d2, dd2, dc, dcc, dc, dcc, p, hx, tx, hy, ty, q, c, cc);
 150
 151   MUL2 (cs, ccs, ds, dss, e, ee, p, hx, tx, hy, ty, q, c, cc);
 152   MUL2 (dc, dcc, sn, ssn, dc, dcc, p, hx, tx, hy, ty, q, c, cc);
 153
 154   MUL2 (d2, dd2, s7.x, ss7.x, ds, dss, p, hx, tx, hy, ty, q, c, cc);
 155   ADD2 (ds, dss, s5.x, ss5.x, ds, dss, r, s);
 156   MUL2 (d2, dd2, ds, dss, ds, dss, p, hx, tx, hy, ty, q, c, cc);
 157   ADD2 (ds, dss, s3.x, ss3.x, ds, dss, r, s);
 158   MUL2 (d2, dd2, ds, dss, ds, dss, p, hx, tx, hy, ty, q, c, cc);
 159   MUL2 (d, dd, ds, dss, ds, dss, p, hx, tx, hy, ty, q, c, cc);
 160   ADD2 (ds, dss, d, dd, ds, dss, r, s);
 161   MUL2 (d2, dd2, c8.x, cc8.x, dc, dcc, p, hx, tx, hy, ty, q, c, cc);
 162   ADD2 (dc, dcc, c6.x, cc6.x, dc, dcc, r, s);
 163   MUL2 (d2, dd2, dc, dcc, dc, dcc, p, hx, tx, hy, ty, q, c, cc);
 164   ADD2 (dc, dcc, c4.x, cc4.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, c2.x, cc2.x, dc, dcc, r, s);
 167   MUL2 (d2, dd2, dc, dcc, dc, dcc, p, hx, tx, hy, ty, q, c, cc);
 168   MUL2 (sn, ssn, ds, dss, e, ee, p, hx, tx, hy, ty, q, c, cc);
 169   MUL2 (dc, dcc, cs, ccs, dc, dcc, p, hx, tx, hy, ty, q, c, cc);
 170   ADD2 (e, ee, dc, dcc, e, ee, r, s);
 171   SUB2 (cs, ccs, e, ee, e, ee, r, s);
 172
 173   v[0] = e;
 174   v[1] = ee;
 175 }
 176 /**********************************************************************/
 177 /* Routine receive Double-Length number (x+dx) and computes cos(x+dx) */
 178 /* as Double-Length number and store it in array v                    */
 179 /**********************************************************************/
 180 void
 181 SECTION
 182 __docos (double x, double dx, double v[])
 183 {
 184   double y, yy, p, w[2];
 185   if (x > 0)
 186     {
 187       y = x; yy = dx;
 188     }
 189   else
 190     {
 191       y = -x; yy = -dx;
 192     }
 193   if (y < 0.5 * hp0.x)                                 /*  y< PI/4    */
 194     {
 195       __dubcos (y, yy, w); v[0] = w[0]; v[1] = w[1];
 196     }
 197   else if (y < 1.5 * hp0.x)                        /* y< 3/4 * PI */
 198     {
 199       p = hp0.x - y; /* p = PI/2 - y */
 200       yy = hp1.x - yy;
 201       y = p + yy;
 202       yy = (p - y) + yy;
 203       if (y > 0)
 204         {
 205           __dubsin (y, yy, w); v[0] = w[0]; v[1] = w[1];
 206         }
 207       /* cos(x) = sin ( 90 -  x ) */
 208       else
 209         {
 210           __dubsin (-y, -yy, w); v[0] = -w[0]; v[1] = -w[1];
 211         }
 212     }
 213   else   /* y>= 3/4 * PI */
 214     {
 215       p = 2.0 * hp0.x - y; /* p = PI- y */
 216       yy = 2.0 * hp1.x - yy;
 217       y = p + yy;
 218       yy = (p - y) + yy;
 219       __dubcos (y, yy, w);
 220       v[0] = -w[0];
 221       v[1] = -w[1];
 222     }
 223 }