libgcc-math/dbl-64/sincos32.c

   1 /*
   2  * IBM Accurate Mathematical Library
   3  * written by International Business Machines Corp.
   4  * Copyright (C) 2001 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, write to the Free Software
  18  * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
  19  */
  20 /****************************************************************/
  21 /*  MODULE_NAME: sincos32.c                                     */
  22 /*                                                              */
  23 /*  FUNCTIONS: ss32                                             */
  24 /*             cc32                                             */
  25 /*             c32                                              */
  26 /*             sin32                                            */
  27 /*             cos32                                            */
  28 /*             mpsin                                            */
  29 /*             mpcos                                            */
  30 /*             mpranred                                         */
  31 /*             mpsin1                                           */
  32 /*             mpcos1                                           */
  33 /*                                                              */
  34 /* FILES NEEDED: endian.h mpa.h sincos32.h                      */
  35 /*               mpa.c                                          */
  36 /*                                                              */
  37 /* Multi Precision sin() and cos() function with p=32  for sin()*/
  38 /* cos() arcsin() and arccos() routines                         */
  39 /* In addition mpranred() routine  performs range  reduction of */
  40 /* a double number x into multi precision number   y,           */
  41 /* such that y=x-n*pi/2, abs(y)<pi/4,  n=0,+-1,+-2,....         */
  42 /****************************************************************/
  43 #include "endian.h"
  44 #include "mpa.h"
  45 #include "sincos32.h"
  46 #include "math_private.h"
  47
  48 /****************************************************************/
  49 /* Compute Multi-Precision sin() function for given p.  Receive */
  50 /* Multi  Precision number x and result stored at y             */
  51 /****************************************************************/
  52 static void ss32(mp_no *x, mp_no *y, int p) {
  53   int i;
  54   double a;
  55 #if 0
  56   double b;
  57   static const mp_no mpone = {1,{1.0,1.0}};
  58 #endif
  59   mp_no mpt1,x2,gor,sum ,mpk={1,{1.0}};
  60 #if 0
  61   mp_no mpt2;
  62 #endif
  63   for (i=1;i<=p;i++) mpk.d[i]=0;
  64
  65   __mul(x,x,&x2,p);
  66   __cpy(&oofac27,&gor,p);
  67   __cpy(&gor,&sum,p);
  68   for (a=27.0;a>1.0;a-=2.0) {
  69     mpk.d[1]=a*(a-1.0);
  70     __mul(&gor,&mpk,&mpt1,p);
  71     __cpy(&mpt1,&gor,p);
  72     __mul(&x2,&sum,&mpt1,p);
  73     __sub(&gor,&mpt1,&sum,p);
  74   }
  75   __mul(x,&sum,y,p);
  76 }
  77
  78 /**********************************************************************/
  79 /* Compute Multi-Precision cos() function for given p. Receive Multi  */
  80 /* Precision number x and result stored at y                          */
  81 /**********************************************************************/
  82 static void cc32(mp_no *x, mp_no *y, int p) {
  83   int i;
  84   double a;
  85 #if 0
  86   double b;
  87   static const mp_no mpone = {1,{1.0,1.0}};
  88 #endif
  89   mp_no mpt1,x2,gor,sum ,mpk={1,{1.0}};
  90 #if 0
  91   mp_no mpt2;
  92 #endif
  93   for (i=1;i<=p;i++) mpk.d[i]=0;
  94
  95   __mul(x,x,&x2,p);
  96   mpk.d[1]=27.0;
  97   __mul(&oofac27,&mpk,&gor,p);
  98   __cpy(&gor,&sum,p);
  99   for (a=26.0;a>2.0;a-=2.0) {
 100     mpk.d[1]=a*(a-1.0);
 101     __mul(&gor,&mpk,&mpt1,p);
 102     __cpy(&mpt1,&gor,p);
 103     __mul(&x2,&sum,&mpt1,p);
 104     __sub(&gor,&mpt1,&sum,p);
 105   }
 106   __mul(&x2,&sum,y,p);
 107 }
 108
 109 /***************************************************************************/
 110 /* c32()   computes both sin(x), cos(x) as Multi precision numbers         */
 111 /***************************************************************************/
 112 void __c32(mp_no *x, mp_no *y, mp_no *z, int p) {
 113   static const mp_no mpt={1,{1.0,2.0}}, one={1,{1.0,1.0}};
 114   mp_no u,t,t1,t2,c,s;
 115   int i;
 116   __cpy(x,&u,p);
 117   u.e=u.e-1;
 118   cc32(&u,&c,p);
 119   ss32(&u,&s,p);
 120   for (i=0;i<24;i++) {
 121     __mul(&c,&s,&t,p);
 122     __sub(&s,&t,&t1,p);
 123     __add(&t1,&t1,&s,p);
 124     __sub(&mpt,&c,&t1,p);
 125     __mul(&t1,&c,&t2,p);
 126     __add(&t2,&t2,&c,p);
 127   }
 128   __sub(&one,&c,y,p);
 129   __cpy(&s,z,p);
 130 }
 131
 132 /************************************************************************/
 133 /*Routine receive double x and two double results of sin(x) and return  */
 134 /*result which is more accurate                                         */
 135 /*Computing sin(x) with multi precision routine c32                     */
 136 /************************************************************************/
 137 double __sin32(double x, double res, double res1) {
 138   int p;
 139   mp_no a,b,c;
 140   p=32;
 141   __dbl_mp(res,&a,p);
 142   __dbl_mp(0.5*(res1-res),&b,p);
 143   __add(&a,&b,&c,p);
 144   if (x>0.8)
 145   { __sub(&hp,&c,&a,p);
 146     __c32(&a,&b,&c,p);
 147   }
 148   else __c32(&c,&a,&b,p);     /* b=sin(0.5*(res+res1))  */
 149   __dbl_mp(x,&c,p);           /* c = x                  */
 150   __sub(&b,&c,&a,p);
 151   /* if a>0 return min(res,res1), otherwise return max(res,res1) */
 152   if (a.d[0]>0)  return (res<res1)?res:res1;
 153   else  return (res>res1)?res:res1;
 154 }
 155
 156 /************************************************************************/
 157 /*Routine receive double x and two double results of cos(x) and return  */
 158 /*result which is more accurate                                         */
 159 /*Computing cos(x) with multi precision routine c32                     */
 160 /************************************************************************/
 161 double __cos32(double x, double res, double res1) {
 162   int p;
 163   mp_no a,b,c;
 164   p=32;
 165   __dbl_mp(res,&a,p);
 166   __dbl_mp(0.5*(res1-res),&b,p);
 167   __add(&a,&b,&c,p);
 168   if (x>2.4)
 169   { __sub(&pi,&c,&a,p);
 170     __c32(&a,&b,&c,p);
 171     b.d[0]=-b.d[0];
 172   }
 173   else if (x>0.8)
 174        { __sub(&hp,&c,&a,p);
 175          __c32(&a,&c,&b,p);
 176        }
 177   else __c32(&c,&b,&a,p);     /* b=cos(0.5*(res+res1))  */
 178   __dbl_mp(x,&c,p);    /* c = x                  */
 179   __sub(&b,&c,&a,p);
 180              /* if a>0 return max(res,res1), otherwise return min(res,res1) */
 181   if (a.d[0]>0)  return (res>res1)?res:res1;
 182   else  return (res<res1)?res:res1;
 183 }
 184
 185 /*******************************************************************/
 186 /*Compute sin(x+dx) as Multi Precision number and return result as */
 187 /* double                                                          */
 188 /*******************************************************************/
 189 double __mpsin(double x, double dx) {
 190   int p;
 191   double y;
 192   mp_no a,b,c;
 193   p=32;
 194   __dbl_mp(x,&a,p);
 195   __dbl_mp(dx,&b,p);
 196   __add(&a,&b,&c,p);
 197   if (x>0.8) { __sub(&hp,&c,&a,p); __c32(&a,&b,&c,p); }
 198   else __c32(&c,&a,&b,p);     /* b = sin(x+dx)     */
 199   __mp_dbl(&b,&y,p);
 200   return y;
 201 }
 202
 203 /*******************************************************************/
 204 /* Compute cos()of double-length number (x+dx) as Multi Precision  */
 205 /* number and return result as double                              */
 206 /*******************************************************************/
 207 double __mpcos(double x, double dx) {
 208   int p;
 209   double y;
 210   mp_no a,b,c;
 211   p=32;
 212   __dbl_mp(x,&a,p);
 213   __dbl_mp(dx,&b,p);
 214   __add(&a,&b,&c,p);
 215   if (x>0.8)
 216   { __sub(&hp,&c,&b,p);
 217     __c32(&b,&c,&a,p);
 218   }
 219   else __c32(&c,&a,&b,p);     /* a = cos(x+dx)     */
 220   __mp_dbl(&a,&y,p);
 221   return y;
 222 }
 223
 224 /******************************************************************/
 225 /* mpranred() performs range reduction of a double number x into  */
 226 /* multi precision number y, such that y=x-n*pi/2, abs(y)<pi/4,   */
 227 /* n=0,+-1,+-2,....                                               */
 228 /* Return int which indicates in which quarter of circle x is     */
 229 /******************************************************************/
 230 int __mpranred(double x, mp_no *y, int p)
 231 {
 232   number v;
 233   double t,xn;
 234   int i,k,n;
 235   static const mp_no one = {1,{1.0,1.0}};
 236   mp_no a,b,c;
 237
 238   if (ABS(x) < 2.8e14) {
 239     t = (x*hpinv.d + toint.d);
 240     xn = t - toint.d;
 241     v.d = t;
 242     n =v.i[LOW_HALF]&3;
 243     __dbl_mp(xn,&a,p);
 244     __mul(&a,&hp,&b,p);
 245     __dbl_mp(x,&c,p);
 246     __sub(&c,&b,y,p);
 247     return n;
 248   }
 249   else {                      /* if x is very big more precision required */
 250     __dbl_mp(x,&a,p);
 251     a.d[0]=1.0;
 252     k = a.e-5;
 253     if (k < 0) k=0;
 254     b.e = -k;
 255     b.d[0] = 1.0;
 256     for (i=0;i<p;i++) b.d[i+1] = toverp[i+k];
 257     __mul(&a,&b,&c,p);
 258     t = c.d[c.e];
 259     for (i=1;i<=p-c.e;i++) c.d[i]=c.d[i+c.e];
 260     for (i=p+1-c.e;i<=p;i++) c.d[i]=0;
 261     c.e=0;
 262     if (c.d[1] >=  8388608.0)
 263     { t +=1.0;
 264       __sub(&c,&one,&b,p);
 265       __mul(&b,&hp,y,p);
 266     }
 267     else __mul(&c,&hp,y,p);
 268     n = (int) t;
 269     if (x < 0) { y->d[0] = - y->d[0]; n = -n; }
 270     return (n&3);
 271   }
 272 }
 273
 274 /*******************************************************************/
 275 /* Multi-Precision sin() function subroutine, for p=32.  It is     */
 276 /* based on the routines mpranred() and c32().                     */
 277 /*******************************************************************/
 278 double __mpsin1(double x)
 279 {
 280   int p;
 281   int n;
 282   mp_no u,s,c;
 283   double y;
 284   p=32;
 285   n=__mpranred(x,&u,p);               /* n is 0, 1, 2 or 3 */
 286   __c32(&u,&c,&s,p);
 287   switch (n) {                      /* in which quarter of unit circle y is*/
 288   case 0:
 289     __mp_dbl(&s,&y,p);
 290     return y;
 291     break;
 292
 293   case 2:
 294     __mp_dbl(&s,&y,p);
 295     return -y;
 296     break;
 297
 298   case 1:
 299     __mp_dbl(&c,&y,p);
 300     return y;
 301     break;
 302
 303   case 3:
 304     __mp_dbl(&c,&y,p);
 305     return -y;
 306     break;
 307
 308   }
 309   return 0;                     /* unreachable, to make the compiler happy */
 310 }
 311
 312 /*****************************************************************/
 313 /* Multi-Precision cos() function subroutine, for p=32.  It is   */
 314 /* based  on the routines mpranred() and c32().                  */
 315 /*****************************************************************/
 316
 317 double __mpcos1(double x)
 318 {
 319   int p;
 320   int n;
 321   mp_no u,s,c;
 322   double y;
 323
 324   p=32;
 325   n=__mpranred(x,&u,p);              /* n is 0, 1, 2 or 3 */
 326   __c32(&u,&c,&s,p);
 327   switch (n) {                     /* in what quarter of unit circle y is*/
 328
 329   case 0:
 330     __mp_dbl(&c,&y,p);
 331     return y;
 332     break;
 333
 334   case 2:
 335     __mp_dbl(&c,&y,p);
 336     return -y;
 337     break;
 338
 339   case 1:
 340     __mp_dbl(&s,&y,p);
 341     return -y;
 342     break;
 343
 344   case 3:
 345     __mp_dbl(&s,&y,p);
 346     return y;
 347     break;
 348
 349   }
 350   return 0;                     /* unreachable, to make the compiler happy */
 351 }
 352 /******************************************************************/