libgcc-math/dbl-64/e_pow.c

   1 /*
   2  * IBM Accurate Mathematical Library
   3  * written by International Business Machines Corp.
   4  * Copyright (C) 2001, 2002, 2004 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: upow.c                                                    */
  22 /*                                                                         */
  23 /*  FUNCTIONS: upow                                                        */
  24 /*             power1                                                      */
  25 /*             my_log2                                                        */
  26 /*             log1                                                        */
  27 /*             checkint                                                    */
  28 /* FILES NEEDED: dla.h endian.h mpa.h mydefs.h                             */
  29 /*               halfulp.c mpexp.c mplog.c slowexp.c slowpow.c mpa.c       */
  30 /*                          uexp.c  upow.c                                 */
  31 /*               root.tbl uexp.tbl upow.tbl                                */
  32 /* An ultimate power routine. Given two IEEE double machine numbers y,x    */
  33 /* it computes the correctly rounded (to nearest) value of x^y.            */
  34 /* Assumption: Machine arithmetic operations are performed in              */
  35 /* round to nearest mode of IEEE 754 standard.                             */
  36 /*                                                                         */
  37 /***************************************************************************/
  38 #include "endian.h"
  39 #include "upow.h"
  40 #include "dla.h"
  41 #include "mydefs.h"
  42 #include "MathLib.h"
  43 #include "upow.tbl"
  44 #include "math_private.h"
  45
  46
  47 double __exp1(double x, double xx, double error);
  48 static double log1(double x, double *delta, double *error);
  49 static double my_log2(double x, double *delta, double *error);
  50 double __slowpow(double x, double y,double z);
  51 static double power1(double x, double y);
  52 static int checkint(double x);
  53
  54 /***************************************************************************/
  55 /* An ultimate power routine. Given two IEEE double machine numbers y,x    */
  56 /* it computes the correctly rounded (to nearest) value of X^y.            */
  57 /***************************************************************************/
  58 double __ieee754_pow(double x, double y) {
  59   double z,a,aa,error, t,a1,a2,y1,y2;
  60 #if 0
  61   double gor=1.0;
  62 #endif
  63   mynumber u,v;
  64   int k;
  65   int4 qx,qy;
  66   v.x=y;
  67   u.x=x;
  68   if (v.i[LOW_HALF] == 0) { /* of y */
  69     qx = u.i[HIGH_HALF]&0x7fffffff;
  70     /* Checking  if x is not too small to compute */
  71     if (((qx==0x7ff00000)&&(u.i[LOW_HALF]!=0))||(qx>0x7ff00000)) return NaNQ.x;
  72     if (y == 1.0) return x;
  73     if (y == 2.0) return x*x;
  74     if (y == -1.0) return 1.0/x;
  75     if (y == 0) return 1.0;
  76   }
  77   /* else */
  78   if(((u.i[HIGH_HALF]>0 && u.i[HIGH_HALF]<0x7ff00000)||        /* x>0 and not x->0 */
  79        (u.i[HIGH_HALF]==0 && u.i[LOW_HALF]!=0))  &&
  80                                       /*   2^-1023< x<= 2^-1023 * 0x1.0000ffffffff */
  81       (v.i[HIGH_HALF]&0x7fffffff) < 0x4ff00000) {              /* if y<-1 or y>1   */
  82     z = log1(x,&aa,&error);                                 /* x^y  =e^(y log (X)) */
  83     t = y*134217729.0;
  84     y1 = t - (t-y);
  85     y2 = y - y1;
  86     t = z*134217729.0;
  87     a1 = t - (t-z);
  88     a2 = (z - a1)+aa;
  89     a = y1*a1;
  90     aa = y2*a1 + y*a2;
  91     a1 = a+aa;
  92     a2 = (a-a1)+aa;
  93     error = error*ABS(y);
  94     t = __exp1(a1,a2,1.9e16*error);     /* return -10 or 0 if wasn't computed exactly */
  95     return (t>0)?t:power1(x,y);
  96   }
  97
  98   if (x == 0) {
  99     if (((v.i[HIGH_HALF] & 0x7fffffff) == 0x7ff00000 && v.i[LOW_HALF] != 0)
 100         || (v.i[HIGH_HALF] & 0x7fffffff) > 0x7ff00000)
 101       return y;
 102     if (ABS(y) > 1.0e20) return (y>0)?0:INF.x;
 103     k = checkint(y);
 104     if (k == -1)
 105       return y < 0 ? 1.0/x : x;
 106     else
 107       return y < 0 ? 1.0/ABS(x) : 0.0;                               /* return 0 */
 108   }
 109   /* if x<0 */
 110   if (u.i[HIGH_HALF] < 0) {
 111     k = checkint(y);
 112     if (k==0) {
 113       if ((v.i[HIGH_HALF] & 0x7fffffff) == 0x7ff00000 && v.i[LOW_HALF] == 0) {
 114         if (x == -1.0) return 1.0;
 115         else if (x > -1.0) return v.i[HIGH_HALF] < 0 ? INF.x : 0.0;
 116         else return v.i[HIGH_HALF] < 0 ? 0.0 : INF.x;
 117       }
 118       else if (u.i[HIGH_HALF] == 0xfff00000 && u.i[LOW_HALF] == 0)
 119         return y < 0 ? 0.0 : INF.x;
 120       return NaNQ.x;                              /* y not integer and x<0 */
 121     }
 122     else if (u.i[HIGH_HALF] == 0xfff00000 && u.i[LOW_HALF] == 0)
 123       {
 124         if (k < 0)
 125           return y < 0 ? nZERO.x : nINF.x;
 126         else
 127           return y < 0 ? 0.0 : INF.x;
 128       }
 129     return (k==1)?__ieee754_pow(-x,y):-__ieee754_pow(-x,y); /* if y even or odd */
 130   }
 131   /* x>0 */
 132   qx = u.i[HIGH_HALF]&0x7fffffff;  /*   no sign   */
 133   qy = v.i[HIGH_HALF]&0x7fffffff;  /*   no sign   */
 134
 135   if (qx > 0x7ff00000 || (qx == 0x7ff00000 && u.i[LOW_HALF] != 0)) return NaNQ.x;
 136                                                                  /*  if 0<x<2^-0x7fe */
 137   if (qy > 0x7ff00000 || (qy == 0x7ff00000 && v.i[LOW_HALF] != 0))
 138     return x == 1.0 ? 1.0 : NaNQ.x;
 139                                                                  /*  if y<2^-0x7fe   */
 140
 141   if (qx == 0x7ff00000)                              /* x= 2^-0x3ff */
 142     {if (y == 0) return NaNQ.x;
 143     return (y>0)?x:0; }
 144
 145   if (qy > 0x45f00000 && qy < 0x7ff00000) {
 146     if (x == 1.0) return 1.0;
 147     if (y>0) return (x>1.0)?INF.x:0;
 148     if (y<0) return (x<1.0)?INF.x:0;
 149   }
 150
 151   if (x == 1.0) return 1.0;
 152   if (y>0) return (x>1.0)?INF.x:0;
 153   if (y<0) return (x<1.0)?INF.x:0;
 154   return 0;     /* unreachable, to make the compiler happy */
 155 }
 156
 157 /**************************************************************************/
 158 /* Computing x^y using more accurate but more slow log routine            */
 159 /**************************************************************************/
 160 static double power1(double x, double y) {
 161   double z,a,aa,error, t,a1,a2,y1,y2;
 162   z = my_log2(x,&aa,&error);
 163   t = y*134217729.0;
 164   y1 = t - (t-y);
 165   y2 = y - y1;
 166   t = z*134217729.0;
 167   a1 = t - (t-z);
 168   a2 = z - a1;
 169   a = y*z;
 170   aa = ((y1*a1-a)+y1*a2+y2*a1)+y2*a2+aa*y;
 171   a1 = a+aa;
 172   a2 = (a-a1)+aa;
 173   error = error*ABS(y);
 174   t = __exp1(a1,a2,1.9e16*error);
 175   return (t >= 0)?t:__slowpow(x,y,z);
 176 }
 177
 178 /****************************************************************************/
 179 /* Computing log(x) (x is left argument). The result is the returned double */
 180 /* + the parameter delta.                                                   */
 181 /* The result is bounded by error (rightmost argument)                      */
 182 /****************************************************************************/
 183 static double log1(double x, double *delta, double *error) {
 184   int i,j,m;
 185 #if 0
 186   int n;
 187 #endif
 188   double uu,vv,eps,nx,e,e1,e2,t,t1,t2,res,add=0;
 189 #if 0
 190   double cor;
 191 #endif
 192   mynumber u,v;
 193 #ifdef BIG_ENDI
 194   mynumber
 195 /**/ two52          = {{0x43300000, 0x00000000}}; /* 2**52         */
 196 #else
 197 #ifdef LITTLE_ENDI
 198   mynumber
 199 /**/ two52          = {{0x00000000, 0x43300000}}; /* 2**52         */
 200 #endif
 201 #endif
 202
 203   u.x = x;
 204   m = u.i[HIGH_HALF];
 205   *error = 0;
 206   *delta = 0;
 207   if (m < 0x00100000)             /*  1<x<2^-1007 */
 208     { x = x*t52.x; add = -52.0; u.x = x; m = u.i[HIGH_HALF];}
 209
 210   if ((m&0x000fffff) < 0x0006a09e)
 211     {u.i[HIGH_HALF] = (m&0x000fffff)|0x3ff00000; two52.i[LOW_HALF]=(m>>20); }
 212   else
 213     {u.i[HIGH_HALF] = (m&0x000fffff)|0x3fe00000; two52.i[LOW_HALF]=(m>>20)+1; }
 214
 215   v.x = u.x + bigu.x;
 216   uu = v.x - bigu.x;
 217   i = (v.i[LOW_HALF]&0x000003ff)<<2;
 218   if (two52.i[LOW_HALF] == 1023)         /* nx = 0              */
 219   {
 220       if (i > 1192 && i < 1208)          /* |x-1| < 1.5*2**-10  */
 221       {
 222           t = x - 1.0;
 223           t1 = (t+5.0e6)-5.0e6;
 224           t2 = t-t1;
 225           e1 = t - 0.5*t1*t1;
 226           e2 = t*t*t*(r3+t*(r4+t*(r5+t*(r6+t*(r7+t*r8)))))-0.5*t2*(t+t1);
 227           res = e1+e2;
 228           *error = 1.0e-21*ABS(t);
 229           *delta = (e1-res)+e2;
 230           return res;
 231       }                  /* |x-1| < 1.5*2**-10  */
 232       else
 233       {
 234           v.x = u.x*(ui.x[i]+ui.x[i+1])+bigv.x;
 235           vv = v.x-bigv.x;
 236           j = v.i[LOW_HALF]&0x0007ffff;
 237           j = j+j+j;
 238           eps = u.x - uu*vv;
 239           e1 = eps*ui.x[i];
 240           e2 = eps*(ui.x[i+1]+vj.x[j]*(ui.x[i]+ui.x[i+1]));
 241           e = e1+e2;
 242           e2 =  ((e1-e)+e2);
 243           t=ui.x[i+2]+vj.x[j+1];
 244           t1 = t+e;
 245           t2 = (((t-t1)+e)+(ui.x[i+3]+vj.x[j+2]))+e2+e*e*(p2+e*(p3+e*p4));
 246           res=t1+t2;
 247           *error = 1.0e-24;
 248           *delta = (t1-res)+t2;
 249           return res;
 250       }
 251   }   /* nx = 0 */
 252   else                            /* nx != 0   */
 253   {
 254       eps = u.x - uu;
 255       nx = (two52.x - two52e.x)+add;
 256       e1 = eps*ui.x[i];
 257       e2 = eps*ui.x[i+1];
 258       e=e1+e2;
 259       e2 = (e1-e)+e2;
 260       t=nx*ln2a.x+ui.x[i+2];
 261       t1=t+e;
 262       t2=(((t-t1)+e)+nx*ln2b.x+ui.x[i+3]+e2)+e*e*(q2+e*(q3+e*(q4+e*(q5+e*q6))));
 263       res = t1+t2;
 264       *error = 1.0e-21;
 265       *delta = (t1-res)+t2;
 266       return res;
 267   }                                /* nx != 0   */
 268 }
 269
 270 /****************************************************************************/
 271 /* More slow but more accurate routine of log                               */
 272 /* Computing log(x)(x is left argument).The result is return double + delta.*/
 273 /* The result is bounded by error (right argument)                           */
 274 /****************************************************************************/
 275 static double my_log2(double x, double *delta, double *error) {
 276   int i,j,m;
 277 #if 0
 278   int n;
 279 #endif
 280   double uu,vv,eps,nx,e,e1,e2,t,t1,t2,res,add=0;
 281 #if 0
 282   double cor;
 283 #endif
 284   double ou1,ou2,lu1,lu2,ov,lv1,lv2,a,a1,a2;
 285   double y,yy,z,zz,j1,j2,j3,j4,j5,j6,j7,j8;
 286   mynumber u,v;
 287 #ifdef BIG_ENDI
 288   mynumber
 289 /**/ two52          = {{0x43300000, 0x00000000}}; /* 2**52         */
 290 #else
 291 #ifdef LITTLE_ENDI
 292   mynumber
 293 /**/ two52          = {{0x00000000, 0x43300000}}; /* 2**52         */
 294 #endif
 295 #endif
 296
 297   u.x = x;
 298   m = u.i[HIGH_HALF];
 299   *error = 0;
 300   *delta = 0;
 301   add=0;
 302   if (m<0x00100000) {  /* x < 2^-1022 */
 303     x = x*t52.x;  add = -52.0; u.x = x; m = u.i[HIGH_HALF]; }
 304
 305   if ((m&0x000fffff) < 0x0006a09e)
 306     {u.i[HIGH_HALF] = (m&0x000fffff)|0x3ff00000; two52.i[LOW_HALF]=(m>>20); }
 307   else
 308     {u.i[HIGH_HALF] = (m&0x000fffff)|0x3fe00000; two52.i[LOW_HALF]=(m>>20)+1; }
 309
 310   v.x = u.x + bigu.x;
 311   uu = v.x - bigu.x;
 312   i = (v.i[LOW_HALF]&0x000003ff)<<2;
 313   /*------------------------------------- |x-1| < 2**-11-------------------------------  */
 314   if ((two52.i[LOW_HALF] == 1023)  && (i == 1200))
 315   {
 316       t = x - 1.0;
 317       EMULV(t,s3,y,yy,j1,j2,j3,j4,j5);
 318       ADD2(-0.5,0,y,yy,z,zz,j1,j2);
 319       MUL2(t,0,z,zz,y,yy,j1,j2,j3,j4,j5,j6,j7,j8);
 320       MUL2(t,0,y,yy,z,zz,j1,j2,j3,j4,j5,j6,j7,j8);
 321
 322       e1 = t+z;
 323       e2 = (((t-e1)+z)+zz)+t*t*t*(ss3+t*(s4+t*(s5+t*(s6+t*(s7+t*s8)))));
 324       res = e1+e2;
 325       *error = 1.0e-25*ABS(t);
 326       *delta = (e1-res)+e2;
 327       return res;
 328   }
 329   /*----------------------------- |x-1| > 2**-11  --------------------------  */
 330   else
 331   {          /*Computing log(x) according to log table                        */
 332       nx = (two52.x - two52e.x)+add;
 333       ou1 = ui.x[i];
 334       ou2 = ui.x[i+1];
 335       lu1 = ui.x[i+2];
 336       lu2 = ui.x[i+3];
 337       v.x = u.x*(ou1+ou2)+bigv.x;
 338       vv = v.x-bigv.x;
 339       j = v.i[LOW_HALF]&0x0007ffff;
 340       j = j+j+j;
 341       eps = u.x - uu*vv;
 342       ov  = vj.x[j];
 343       lv1 = vj.x[j+1];
 344       lv2 = vj.x[j+2];
 345       a = (ou1+ou2)*(1.0+ov);
 346       a1 = (a+1.0e10)-1.0e10;
 347       a2 = a*(1.0-a1*uu*vv);
 348       e1 = eps*a1;
 349       e2 = eps*a2;
 350       e = e1+e2;
 351       e2 = (e1-e)+e2;
 352       t=nx*ln2a.x+lu1+lv1;
 353       t1 = t+e;
 354       t2 = (((t-t1)+e)+(lu2+lv2+nx*ln2b.x+e2))+e*e*(p2+e*(p3+e*p4));
 355       res=t1+t2;
 356       *error = 1.0e-27;
 357       *delta = (t1-res)+t2;
 358       return res;
 359   }
 360 }
 361
 362 /**********************************************************************/
 363 /* Routine receives a double x and checks if it is an integer. If not */
 364 /* it returns 0, else it returns 1 if even or -1 if odd.              */
 365 /**********************************************************************/
 366 static int checkint(double x) {
 367   union {int4 i[2]; double x;} u;
 368   int k,m,n;
 369 #if 0
 370   int l;
 371 #endif
 372   u.x = x;
 373   m = u.i[HIGH_HALF]&0x7fffffff;    /* no sign */
 374   if (m >= 0x7ff00000) return 0;    /*  x is +/-inf or NaN  */
 375   if (m >= 0x43400000) return 1;    /*  |x| >= 2**53   */
 376   if (m < 0x40000000) return 0;     /* |x| < 2,  can not be 0 or 1  */
 377   n = u.i[LOW_HALF];
 378   k = (m>>20)-1023;                 /*  1 <= k <= 52   */
 379   if (k == 52) return (n&1)? -1:1;  /* odd or even*/
 380   if (k>20) {
 381     if (n<<(k-20)) return 0;        /* if not integer */
 382     return (n<<(k-21))?-1:1;
 383   }
 384   if (n) return 0;                  /*if  not integer*/
 385   if (k == 20) return (m&1)? -1:1;
 386   if (m<<(k+12)) return 0;
 387   return (m<<(k+11))?-1:1;
 388 }