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