Back to index

glibc  2.9
e_pow.c
Go to the documentation of this file.
00001 /*
00002  * IBM Accurate Mathematical Library
00003  * written by International Business Machines Corp.
00004  * Copyright (C) 2001, 2002, 2004 Free Software Foundation
00005  *
00006  * This program is free software; you can redistribute it and/or modify
00007  * it under the terms of the GNU Lesser General Public License as published by
00008  * the Free Software Foundation; either version 2.1 of the License, or
00009  * (at your option) any later version.
00010  *
00011  * This program is distributed in the hope that it will be useful,
00012  * but WITHOUT ANY WARRANTY; without even the implied warranty of
00013  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00014  * GNU Lesser General Public License for more details.
00015  *
00016  * You should have received a copy of the GNU Lesser General Public License
00017  * along with this program; if not, write to the Free Software
00018  * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
00019  */
00020 /***************************************************************************/
00021 /*  MODULE_NAME: upow.c                                                    */
00022 /*                                                                         */
00023 /*  FUNCTIONS: upow                                                        */
00024 /*             power1                                                      */
00025 /*             my_log2                                                        */
00026 /*             log1                                                        */
00027 /*             checkint                                                    */
00028 /* FILES NEEDED: dla.h endian.h mpa.h mydefs.h                             */
00029 /*               halfulp.c mpexp.c mplog.c slowexp.c slowpow.c mpa.c       */
00030 /*                          uexp.c  upow.c                                */
00031 /*               root.tbl uexp.tbl upow.tbl                                */
00032 /* An ultimate power routine. Given two IEEE double machine numbers y,x    */
00033 /* it computes the correctly rounded (to nearest) value of x^y.            */
00034 /* Assumption: Machine arithmetic operations are performed in              */
00035 /* round to nearest mode of IEEE 754 standard.                             */
00036 /*                                                                         */
00037 /***************************************************************************/
00038 #include "endian.h"
00039 #include "upow.h"
00040 #include "dla.h"
00041 #include "mydefs.h"
00042 #include "MathLib.h"
00043 #include "upow.tbl"
00044 #include "math_private.h"
00045 
00046 
00047 double __exp1(double x, double xx, double error);
00048 static double log1(double x, double *delta, double *error);
00049 static double my_log2(double x, double *delta, double *error);
00050 double __slowpow(double x, double y,double z);
00051 static double power1(double x, double y);
00052 static int checkint(double x);
00053 
00054 /***************************************************************************/
00055 /* An ultimate power routine. Given two IEEE double machine numbers y,x    */
00056 /* it computes the correctly rounded (to nearest) value of X^y.            */
00057 /***************************************************************************/
00058 double __ieee754_pow(double x, double y) {
00059   double z,a,aa,error, t,a1,a2,y1,y2;
00060 #if 0
00061   double gor=1.0;
00062 #endif
00063   mynumber u,v;
00064   int k;
00065   int4 qx,qy;
00066   v.x=y;
00067   u.x=x;
00068   if (v.i[LOW_HALF] == 0) { /* of y */
00069     qx = u.i[HIGH_HALF]&0x7fffffff;
00070     /* Checking  if x is not too small to compute */
00071     if (((qx==0x7ff00000)&&(u.i[LOW_HALF]!=0))||(qx>0x7ff00000)) return NaNQ.x;
00072     if (y == 1.0) return x;
00073     if (y == 2.0) return x*x;
00074     if (y == -1.0) return 1.0/x;
00075     if (y == 0) return 1.0;
00076   }
00077   /* else */
00078   if(((u.i[HIGH_HALF]>0 && u.i[HIGH_HALF]<0x7ff00000)||        /* x>0 and not x->0 */
00079        (u.i[HIGH_HALF]==0 && u.i[LOW_HALF]!=0))  &&
00080                                       /*   2^-1023< x<= 2^-1023 * 0x1.0000ffffffff */
00081       (v.i[HIGH_HALF]&0x7fffffff) < 0x4ff00000) {              /* if y<-1 or y>1   */
00082     z = log1(x,&aa,&error);                                 /* x^y  =e^(y log (X)) */
00083     t = y*134217729.0;
00084     y1 = t - (t-y);
00085     y2 = y - y1;
00086     t = z*134217729.0;
00087     a1 = t - (t-z);
00088     a2 = (z - a1)+aa;
00089     a = y1*a1;
00090     aa = y2*a1 + y*a2;
00091     a1 = a+aa;
00092     a2 = (a-a1)+aa;
00093     error = error*ABS(y);
00094     t = __exp1(a1,a2,1.9e16*error);     /* return -10 or 0 if wasn't computed exactly */
00095     return (t>0)?t:power1(x,y);
00096   }
00097 
00098   if (x == 0) {
00099     if (((v.i[HIGH_HALF] & 0x7fffffff) == 0x7ff00000 && v.i[LOW_HALF] != 0)
00100        || (v.i[HIGH_HALF] & 0x7fffffff) > 0x7ff00000)
00101       return y;
00102     if (ABS(y) > 1.0e20) return (y>0)?0:INF.x;
00103     k = checkint(y);
00104     if (k == -1)
00105       return y < 0 ? 1.0/x : x;
00106     else
00107       return y < 0 ? 1.0/ABS(x) : 0.0;                               /* return 0 */
00108   }
00109 
00110   qx = u.i[HIGH_HALF]&0x7fffffff;  /*   no sign   */
00111   qy = v.i[HIGH_HALF]&0x7fffffff;  /*   no sign   */
00112 
00113   if (qx >= 0x7ff00000 && (qx > 0x7ff00000 || u.i[LOW_HALF] != 0)) return NaNQ.x;
00114   if (qy >= 0x7ff00000 && (qy > 0x7ff00000 || v.i[LOW_HALF] != 0))
00115     return x == 1.0 ? 1.0 : NaNQ.x;
00116 
00117   /* if x<0 */
00118   if (u.i[HIGH_HALF] < 0) {
00119     k = checkint(y);
00120     if (k==0) {
00121       if (qy == 0x7ff00000) {
00122        if (x == -1.0) return 1.0;
00123        else if (x > -1.0) return v.i[HIGH_HALF] < 0 ? INF.x : 0.0;
00124        else return v.i[HIGH_HALF] < 0 ? 0.0 : INF.x;
00125       }
00126       else if (qx == 0x7ff00000)
00127        return y < 0 ? 0.0 : INF.x;
00128       return NaNQ.x;                              /* y not integer and x<0 */
00129     }
00130     else if (qx == 0x7ff00000)
00131       {
00132        if (k < 0)
00133          return y < 0 ? nZERO.x : nINF.x;
00134        else
00135          return y < 0 ? 0.0 : INF.x;
00136       }
00137     return (k==1)?__ieee754_pow(-x,y):-__ieee754_pow(-x,y); /* if y even or odd */
00138   }
00139   /* x>0 */
00140 
00141   if (qx == 0x7ff00000)                              /* x= 2^-0x3ff */
00142     {if (y == 0) return NaNQ.x;
00143     return (y>0)?x:0; }
00144 
00145   if (qy > 0x45f00000 && qy < 0x7ff00000) {
00146     if (x == 1.0) return 1.0;
00147     if (y>0) return (x>1.0)?INF.x:0;
00148     if (y<0) return (x<1.0)?INF.x:0;
00149   }
00150 
00151   if (x == 1.0) return 1.0;
00152   if (y>0) return (x>1.0)?INF.x:0;
00153   if (y<0) return (x<1.0)?INF.x:0;
00154   return 0;     /* unreachable, to make the compiler happy */
00155 }
00156 
00157 /**************************************************************************/
00158 /* Computing x^y using more accurate but more slow log routine            */
00159 /**************************************************************************/
00160 static double power1(double x, double y) {
00161   double z,a,aa,error, t,a1,a2,y1,y2;
00162   z = my_log2(x,&aa,&error);
00163   t = y*134217729.0;
00164   y1 = t - (t-y);
00165   y2 = y - y1;
00166   t = z*134217729.0;
00167   a1 = t - (t-z);
00168   a2 = z - a1;
00169   a = y*z;
00170   aa = ((y1*a1-a)+y1*a2+y2*a1)+y2*a2+aa*y;
00171   a1 = a+aa;
00172   a2 = (a-a1)+aa;
00173   error = error*ABS(y);
00174   t = __exp1(a1,a2,1.9e16*error);
00175   return (t >= 0)?t:__slowpow(x,y,z);
00176 }
00177 
00178 /****************************************************************************/
00179 /* Computing log(x) (x is left argument). The result is the returned double */
00180 /* + the parameter delta.                                                   */
00181 /* The result is bounded by error (rightmost argument)                      */
00182 /****************************************************************************/
00183 static double log1(double x, double *delta, double *error) {
00184   int i,j,m;
00185 #if 0
00186   int n;
00187 #endif
00188   double uu,vv,eps,nx,e,e1,e2,t,t1,t2,res,add=0;
00189 #if 0
00190   double cor;
00191 #endif
00192   mynumber u,v;
00193 #ifdef BIG_ENDI
00194   mynumber
00195  two52          = {{0x43300000, 0x00000000}}; /* 2**52         */
00196 #else
00197 #ifdef LITTLE_ENDI
00198   mynumber
00199  two52          = {{0x00000000, 0x43300000}}; /* 2**52         */
00200 #endif
00201 #endif
00202 
00203   u.x = x;
00204   m = u.i[HIGH_HALF];
00205   *error = 0;
00206   *delta = 0;
00207   if (m < 0x00100000)             /*  1<x<2^-1007 */
00208     { x = x*t52.x; add = -52.0; u.x = x; m = u.i[HIGH_HALF];}
00209 
00210   if ((m&0x000fffff) < 0x0006a09e)
00211     {u.i[HIGH_HALF] = (m&0x000fffff)|0x3ff00000; two52.i[LOW_HALF]=(m>>20); }
00212   else
00213     {u.i[HIGH_HALF] = (m&0x000fffff)|0x3fe00000; two52.i[LOW_HALF]=(m>>20)+1; }
00214 
00215   v.x = u.x + bigu.x;
00216   uu = v.x - bigu.x;
00217   i = (v.i[LOW_HALF]&0x000003ff)<<2;
00218   if (two52.i[LOW_HALF] == 1023)         /* nx = 0              */
00219   {
00220       if (i > 1192 && i < 1208)          /* |x-1| < 1.5*2**-10  */
00221       {
00222          t = x - 1.0;
00223          t1 = (t+5.0e6)-5.0e6;
00224          t2 = t-t1;
00225          e1 = t - 0.5*t1*t1;
00226          e2 = t*t*t*(r3+t*(r4+t*(r5+t*(r6+t*(r7+t*r8)))))-0.5*t2*(t+t1);
00227          res = e1+e2;
00228          *error = 1.0e-21*ABS(t);
00229          *delta = (e1-res)+e2;
00230          return res;
00231       }                  /* |x-1| < 1.5*2**-10  */
00232       else
00233       {
00234          v.x = u.x*(ui.x[i]+ui.x[i+1])+bigv.x;
00235          vv = v.x-bigv.x;
00236          j = v.i[LOW_HALF]&0x0007ffff;
00237          j = j+j+j;
00238          eps = u.x - uu*vv;
00239          e1 = eps*ui.x[i];
00240          e2 = eps*(ui.x[i+1]+vj.x[j]*(ui.x[i]+ui.x[i+1]));
00241          e = e1+e2;
00242          e2 =  ((e1-e)+e2);
00243          t=ui.x[i+2]+vj.x[j+1];
00244          t1 = t+e;
00245          t2 = (((t-t1)+e)+(ui.x[i+3]+vj.x[j+2]))+e2+e*e*(p2+e*(p3+e*p4));
00246          res=t1+t2;
00247          *error = 1.0e-24;
00248          *delta = (t1-res)+t2;
00249          return res;
00250       }
00251   }   /* nx = 0 */
00252   else                            /* nx != 0   */
00253   {
00254       eps = u.x - uu;
00255       nx = (two52.x - two52e.x)+add;
00256       e1 = eps*ui.x[i];
00257       e2 = eps*ui.x[i+1];
00258       e=e1+e2;
00259       e2 = (e1-e)+e2;
00260       t=nx*ln2a.x+ui.x[i+2];
00261       t1=t+e;
00262       t2=(((t-t1)+e)+nx*ln2b.x+ui.x[i+3]+e2)+e*e*(q2+e*(q3+e*(q4+e*(q5+e*q6))));
00263       res = t1+t2;
00264       *error = 1.0e-21;
00265       *delta = (t1-res)+t2;
00266       return res;
00267   }                                /* nx != 0   */
00268 }
00269 
00270 /****************************************************************************/
00271 /* More slow but more accurate routine of log                               */
00272 /* Computing log(x)(x is left argument).The result is return double + delta.*/
00273 /* The result is bounded by error (right argument)                           */
00274 /****************************************************************************/
00275 static double my_log2(double x, double *delta, double *error) {
00276   int i,j,m;
00277 #if 0
00278   int n;
00279 #endif
00280   double uu,vv,eps,nx,e,e1,e2,t,t1,t2,res,add=0;
00281 #if 0
00282   double cor;
00283 #endif
00284   double ou1,ou2,lu1,lu2,ov,lv1,lv2,a,a1,a2;
00285   double y,yy,z,zz,j1,j2,j3,j4,j5,j6,j7,j8;
00286   mynumber u,v;
00287 #ifdef BIG_ENDI
00288   mynumber
00289  two52          = {{0x43300000, 0x00000000}}; /* 2**52         */
00290 #else
00291 #ifdef LITTLE_ENDI
00292   mynumber
00293  two52          = {{0x00000000, 0x43300000}}; /* 2**52         */
00294 #endif
00295 #endif
00296 
00297   u.x = x;
00298   m = u.i[HIGH_HALF];
00299   *error = 0;
00300   *delta = 0;
00301   add=0;
00302   if (m<0x00100000) {  /* x < 2^-1022 */
00303     x = x*t52.x;  add = -52.0; u.x = x; m = u.i[HIGH_HALF]; }
00304 
00305   if ((m&0x000fffff) < 0x0006a09e)
00306     {u.i[HIGH_HALF] = (m&0x000fffff)|0x3ff00000; two52.i[LOW_HALF]=(m>>20); }
00307   else
00308     {u.i[HIGH_HALF] = (m&0x000fffff)|0x3fe00000; two52.i[LOW_HALF]=(m>>20)+1; }
00309 
00310   v.x = u.x + bigu.x;
00311   uu = v.x - bigu.x;
00312   i = (v.i[LOW_HALF]&0x000003ff)<<2;
00313   /*------------------------------------- |x-1| < 2**-11-------------------------------  */
00314   if ((two52.i[LOW_HALF] == 1023)  && (i == 1200))
00315   {
00316       t = x - 1.0;
00317       EMULV(t,s3,y,yy,j1,j2,j3,j4,j5);
00318       ADD2(-0.5,0,y,yy,z,zz,j1,j2);
00319       MUL2(t,0,z,zz,y,yy,j1,j2,j3,j4,j5,j6,j7,j8);
00320       MUL2(t,0,y,yy,z,zz,j1,j2,j3,j4,j5,j6,j7,j8);
00321 
00322       e1 = t+z;
00323       e2 = (((t-e1)+z)+zz)+t*t*t*(ss3+t*(s4+t*(s5+t*(s6+t*(s7+t*s8)))));
00324       res = e1+e2;
00325       *error = 1.0e-25*ABS(t);
00326       *delta = (e1-res)+e2;
00327       return res;
00328   }
00329   /*----------------------------- |x-1| > 2**-11  --------------------------  */
00330   else
00331   {          /*Computing log(x) according to log table                        */
00332       nx = (two52.x - two52e.x)+add;
00333       ou1 = ui.x[i];
00334       ou2 = ui.x[i+1];
00335       lu1 = ui.x[i+2];
00336       lu2 = ui.x[i+3];
00337       v.x = u.x*(ou1+ou2)+bigv.x;
00338       vv = v.x-bigv.x;
00339       j = v.i[LOW_HALF]&0x0007ffff;
00340       j = j+j+j;
00341       eps = u.x - uu*vv;
00342       ov  = vj.x[j];
00343       lv1 = vj.x[j+1];
00344       lv2 = vj.x[j+2];
00345       a = (ou1+ou2)*(1.0+ov);
00346       a1 = (a+1.0e10)-1.0e10;
00347       a2 = a*(1.0-a1*uu*vv);
00348       e1 = eps*a1;
00349       e2 = eps*a2;
00350       e = e1+e2;
00351       e2 = (e1-e)+e2;
00352       t=nx*ln2a.x+lu1+lv1;
00353       t1 = t+e;
00354       t2 = (((t-t1)+e)+(lu2+lv2+nx*ln2b.x+e2))+e*e*(p2+e*(p3+e*p4));
00355       res=t1+t2;
00356       *error = 1.0e-27;
00357       *delta = (t1-res)+t2;
00358       return res;
00359   }
00360 }
00361 
00362 /**********************************************************************/
00363 /* Routine receives a double x and checks if it is an integer. If not */
00364 /* it returns 0, else it returns 1 if even or -1 if odd.              */
00365 /**********************************************************************/
00366 static int checkint(double x) {
00367   union {int4 i[2]; double x;} u;
00368   int k,m,n;
00369 #if 0
00370   int l;
00371 #endif
00372   u.x = x;
00373   m = u.i[HIGH_HALF]&0x7fffffff;    /* no sign */
00374   if (m >= 0x7ff00000) return 0;    /*  x is +/-inf or NaN  */
00375   if (m >= 0x43400000) return 1;    /*  |x| >= 2**53   */
00376   if (m < 0x40000000) return 0;     /* |x| < 2,  can not be 0 or 1  */
00377   n = u.i[LOW_HALF];
00378   k = (m>>20)-1023;                 /*  1 <= k <= 52   */
00379   if (k == 52) return (n&1)? -1:1;  /* odd or even*/
00380   if (k>20) {
00381     if (n<<(k-20)) return 0;        /* if not integer */
00382     return (n<<(k-21))?-1:1;
00383   }
00384   if (n) return 0;                  /*if  not integer*/
00385   if (k == 20) return (m&1)? -1:1;
00386   if (m<<(k+12)) return 0;
00387   return (m<<(k+11))?-1:1;
00388 }