Back to index

lightning-sunbird  0.9+nobinonly
e_pow.c
Go to the documentation of this file.
00001 /* -*- Mode: C; tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 4 -*-
00002  *
00003  * ***** BEGIN LICENSE BLOCK *****
00004  * Version: MPL 1.1/GPL 2.0/LGPL 2.1
00005  *
00006  * The contents of this file are subject to the Mozilla Public License Version
00007  * 1.1 (the "License"); you may not use this file except in compliance with
00008  * the License. You may obtain a copy of the License at
00009  * http://www.mozilla.org/MPL/
00010  *
00011  * Software distributed under the License is distributed on an "AS IS" basis,
00012  * WITHOUT WARRANTY OF ANY KIND, either express or implied. See the License
00013  * for the specific language governing rights and limitations under the
00014  * License.
00015  *
00016  * The Original Code is Mozilla Communicator client code, released
00017  * March 31, 1998.
00018  *
00019  * The Initial Developer of the Original Code is
00020  * Sun Microsystems, Inc.
00021  * Portions created by the Initial Developer are Copyright (C) 1998
00022  * the Initial Developer. All Rights Reserved.
00023  *
00024  * Contributor(s):
00025  *
00026  * Alternatively, the contents of this file may be used under the terms of
00027  * either of the GNU General Public License Version 2 or later (the "GPL"),
00028  * or the GNU Lesser General Public License Version 2.1 or later (the "LGPL"),
00029  * in which case the provisions of the GPL or the LGPL are applicable instead
00030  * of those above. If you wish to allow use of your version of this file only
00031  * under the terms of either the GPL or the LGPL, and not to allow others to
00032  * use your version of this file under the terms of the MPL, indicate your
00033  * decision by deleting the provisions above and replace them with the notice
00034  * and other provisions required by the GPL or the LGPL. If you do not delete
00035  * the provisions above, a recipient may use your version of this file under
00036  * the terms of any one of the MPL, the GPL or the LGPL.
00037  *
00038  * ***** END LICENSE BLOCK ***** */
00039 
00040 /* @(#)e_pow.c 1.3 95/01/18 */
00041 /*
00042  * ====================================================
00043  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
00044  *
00045  * Developed at SunSoft, a Sun Microsystems, Inc. business.
00046  * Permission to use, copy, modify, and distribute this
00047  * software is freely granted, provided that this notice 
00048  * is preserved.
00049  * ====================================================
00050  */
00051 
00052 /* __ieee754_pow(x,y) return x**y
00053  *
00054  *                  n
00055  * Method:  Let x =  2   * (1+f)
00056  *     1. Compute and return log2(x) in two pieces:
00057  *            log2(x) = w1 + w2,
00058  *        where w1 has 53-24 = 29 bit trailing zeros.
00059  *     2. Perform y*log2(x) = n+y' by simulating muti-precision 
00060  *        arithmetic, where |y'|<=0.5.
00061  *     3. Return x**y = 2**n*exp(y'*log2)
00062  *
00063  * Special cases:
00064  *     1.  (anything) ** 0  is 1
00065  *     2.  (anything) ** 1  is itself
00066  *     3.  (anything) ** NAN is NAN
00067  *     4.  NAN ** (anything except 0) is NAN
00068  *     5.  +-(|x| > 1) **  +INF is +INF
00069  *     6.  +-(|x| > 1) **  -INF is +0
00070  *     7.  +-(|x| < 1) **  +INF is +0
00071  *     8.  +-(|x| < 1) **  -INF is +INF
00072  *     9.  +-1         ** +-INF is NAN
00073  *     10. +0 ** (+anything except 0, NAN)               is +0
00074  *     11. -0 ** (+anything except 0, NAN, odd integer)  is +0
00075  *     12. +0 ** (-anything except 0, NAN)               is +INF
00076  *     13. -0 ** (-anything except 0, NAN, odd integer)  is +INF
00077  *     14. -0 ** (odd integer) = -( +0 ** (odd integer) )
00078  *     15. +INF ** (+anything except 0,NAN) is +INF
00079  *     16. +INF ** (-anything except 0,NAN) is +0
00080  *     17. -INF ** (anything)  = -0 ** (-anything)
00081  *     18. (-anything) ** (integer) is (-1)**(integer)*(+anything**integer)
00082  *     19. (-anything except 0 and inf) ** (non-integer) is NAN
00083  *
00084  * Accuracy:
00085  *     pow(x,y) returns x**y nearly rounded. In particular
00086  *                   pow(integer,integer)
00087  *     always returns the correct integer provided it is 
00088  *     representable.
00089  *
00090  * Constants :
00091  * The hexadecimal values are the intended ones for the following 
00092  * constants. The decimal values may be used, provided that the 
00093  * compiler will convert from decimal to binary accurately enough 
00094  * to produce the hexadecimal values shown.
00095  */
00096 
00097 #include "fdlibm.h"
00098 
00099 #if defined(_MSC_VER)
00100 /* Microsoft Compiler */
00101 #pragma warning( disable : 4723 ) /* disables potential divide by 0 warning */
00102 #endif
00103 
00104 #ifdef __STDC__
00105 static const double 
00106 #else
00107 static double 
00108 #endif
00109 bp[] = {1.0, 1.5,},
00110 dp_h[] = { 0.0, 5.84962487220764160156e-01,}, /* 0x3FE2B803, 0x40000000 */
00111 dp_l[] = { 0.0, 1.35003920212974897128e-08,}, /* 0x3E4CFDEB, 0x43CFD006 */
00112 zero    =  0.0,
00113 one    =  1.0,
00114 two    =  2.0,
00115 two53  =  9007199254740992.0,      /* 0x43400000, 0x00000000 */
00116 really_big    =  1.0e300,
00117 tiny    =  1.0e-300,
00118        /* poly coefs for (3/2)*(log(x)-2s-2/3*s**3 */
00119 L1  =  5.99999999999994648725e-01, /* 0x3FE33333, 0x33333303 */
00120 L2  =  4.28571428578550184252e-01, /* 0x3FDB6DB6, 0xDB6FABFF */
00121 L3  =  3.33333329818377432918e-01, /* 0x3FD55555, 0x518F264D */
00122 L4  =  2.72728123808534006489e-01, /* 0x3FD17460, 0xA91D4101 */
00123 L5  =  2.30660745775561754067e-01, /* 0x3FCD864A, 0x93C9DB65 */
00124 L6  =  2.06975017800338417784e-01, /* 0x3FCA7E28, 0x4A454EEF */
00125 P1   =  1.66666666666666019037e-01, /* 0x3FC55555, 0x5555553E */
00126 P2   = -2.77777777770155933842e-03, /* 0xBF66C16C, 0x16BEBD93 */
00127 P3   =  6.61375632143793436117e-05, /* 0x3F11566A, 0xAF25DE2C */
00128 P4   = -1.65339022054652515390e-06, /* 0xBEBBBD41, 0xC5D26BF1 */
00129 P5   =  4.13813679705723846039e-08, /* 0x3E663769, 0x72BEA4D0 */
00130 lg2  =  6.93147180559945286227e-01, /* 0x3FE62E42, 0xFEFA39EF */
00131 lg2_h  =  6.93147182464599609375e-01, /* 0x3FE62E43, 0x00000000 */
00132 lg2_l  = -1.90465429995776804525e-09, /* 0xBE205C61, 0x0CA86C39 */
00133 ovt =  8.0085662595372944372e-0017, /* -(1024-log2(ovfl+.5ulp)) */
00134 cp    =  9.61796693925975554329e-01, /* 0x3FEEC709, 0xDC3A03FD =2/(3ln2) */
00135 cp_h  =  9.61796700954437255859e-01, /* 0x3FEEC709, 0xE0000000 =(float)cp */
00136 cp_l  = -7.02846165095275826516e-09, /* 0xBE3E2FE0, 0x145B01F5 =tail of cp_h*/
00137 ivln2    =  1.44269504088896338700e+00, /* 0x3FF71547, 0x652B82FE =1/ln2 */
00138 ivln2_h  =  1.44269502162933349609e+00, /* 0x3FF71547, 0x60000000 =24b 1/ln2*/
00139 ivln2_l  =  1.92596299112661746887e-08; /* 0x3E54AE0B, 0xF85DDF44 =1/ln2 tail*/
00140 
00141 #ifdef __STDC__
00142        double __ieee754_pow(double x, double y)
00143 #else
00144        double __ieee754_pow(x,y)
00145        double x, y;
00146 #endif
00147 {
00148         fd_twoints ux, uy, uz;
00149         double y1,t1,p_h,t,z,ax;
00150        double z_h,z_l,p_l;
00151        double t2,r,s,u,v,w;
00152        int i,j,k,yisint,n;
00153        int hx,hy,ix,iy;
00154        unsigned lx,ly;
00155 
00156         ux.d = x; uy.d = y;
00157        hx = __HI(ux); lx = __LO(ux);
00158        hy = __HI(uy); ly = __LO(uy);
00159        ix = hx&0x7fffffff;  iy = hy&0x7fffffff;
00160 
00161     /* y==zero: x**0 = 1 */
00162        if((iy|ly)==0) return one;  
00163 
00164     /* +-NaN return x+y */
00165        if(ix > 0x7ff00000 || ((ix==0x7ff00000)&&(lx!=0)) ||
00166           iy > 0x7ff00000 || ((iy==0x7ff00000)&&(ly!=0))) 
00167               return x+y;   
00168 
00169     /* determine if y is an odd int when x < 0
00170      * yisint = 0    ... y is not an integer
00171      * yisint = 1    ... y is an odd int
00172      * yisint = 2    ... y is an even int
00173      */
00174        yisint  = 0;
00175        if(hx<0) {    
00176            if(iy>=0x43400000) yisint = 2; /* even integer y */
00177            else if(iy>=0x3ff00000) {
00178               k = (iy>>20)-0x3ff;     /* exponent */
00179               if(k>20) {
00180                   j = ly>>(52-k);
00181                   if((j<<(52-k))==(int)ly) yisint = 2-(j&1);
00182               } else if(ly==0) {
00183                   j = iy>>(20-k);
00184                   if((j<<(20-k))==iy) yisint = 2-(j&1);
00185               }
00186            }         
00187        } 
00188 
00189     /* special value of y */
00190        if(ly==0) {   
00191            if (iy==0x7ff00000) {   /* y is +-inf */
00192                if(((ix-0x3ff00000)|lx)==0)
00193 #ifdef _WIN32
00194 /* VC++ optimizer reduces y - y to 0 */ 
00195                     return y / y;
00196 #else
00197                   return  y - y;   /* inf**+-1 is NaN */
00198 #endif
00199                else if (ix >= 0x3ff00000)/* (|x|>1)**+-inf = inf,0 */
00200                   return (hy>=0)? y: zero;
00201                else                /* (|x|<1)**-,+inf = inf,0 */
00202                   return (hy<0)?-y: zero;
00203            } 
00204            if(iy==0x3ff00000) {    /* y is  +-1 */
00205               if(hy<0) return one/x; else return x;
00206            }
00207            if(hy==0x40000000) return x*x; /* y is  2 */
00208            if(hy==0x3fe00000) {    /* y is  0.5 */
00209               if(hx>=0)     /* x >= +0 */
00210               return fd_sqrt(x);   
00211            }
00212        }
00213 
00214        ax   = fd_fabs(x);
00215     /* special value of x */
00216        if(lx==0) {
00217            if(ix==0x7ff00000||ix==0||ix==0x3ff00000){
00218               z = ax;                     /*x is +-0,+-inf,+-1*/
00219               if(hy<0) z = one/z;  /* z = (1/|x|) */
00220               if(hx<0) {
00221                   if(((ix-0x3ff00000)|yisint)==0) {
00222                      z = (z-z)/(z-z); /* (-1)**non-int is NaN */
00223                   } else if(yisint==1) {
00224 #ifdef HPUX
00225                         uz.d = z;
00226                      __HI(uz) ^= 1<<31; /* some HPUXes cannot negate 0.. */
00227                         z = uz.d;
00228 #else
00229                      z = -z;              /* (x<0)**odd = -(|x|**odd) */
00230 #endif
00231                      }
00232               }
00233               return z;
00234            }
00235        }
00236     
00237     /* (x<0)**(non-int) is NaN */
00238        if((((hx>>31)+1)|yisint)==0) return (x-x)/(x-x);
00239 
00240     /* |y| is really_big */
00241        if(iy>0x41e00000) { /* if |y| > 2**31 */
00242            if(iy>0x43f00000){      /* if |y| > 2**64, must o/uflow */
00243               if(ix<=0x3fefffff) return (hy<0)? really_big*really_big:tiny*tiny;
00244               if(ix>=0x3ff00000) return (hy>0)? really_big*really_big:tiny*tiny;
00245            }
00246        /* over/underflow if x is not close to one */
00247            if(ix<0x3fefffff) return (hy<0)? really_big*really_big:tiny*tiny;
00248            if(ix>0x3ff00000) return (hy>0)? really_big*really_big:tiny*tiny;
00249        /* now |1-x| is tiny <= 2**-20, suffice to compute 
00250           log(x) by x-x^2/2+x^3/3-x^4/4 */
00251            t = x-1;         /* t has 20 trailing zeros */
00252            w = (t*t)*(0.5-t*(0.3333333333333333333333-t*0.25));
00253            u = ivln2_h*t;   /* ivln2_h has 21 sig. bits */
00254            v = t*ivln2_l-w*ivln2;
00255            t1 = u+v;
00256             uz.d = t1;
00257            __LO(uz) = 0;
00258             t1 = uz.d;
00259            t2 = v-(t1-u);
00260        } else {
00261            double s_h,t_h;
00262            double s2,s_l,t_l;
00263            n = 0;
00264        /* take care subnormal number */
00265            if(ix<0x00100000) 
00266               {ax *= two53; n -= 53; uz.d = ax; ix = __HI(uz); }
00267            n  += ((ix)>>20)-0x3ff;
00268            j  = ix&0x000fffff;
00269        /* determine interval */
00270            ix = j|0x3ff00000;             /* normalize ix */
00271            if(j<=0x3988E) k=0;            /* |x|<sqrt(3/2) */
00272            else if(j<0xBB67A) k=1; /* |x|<sqrt(3)   */
00273            else {k=0;n+=1;ix -= 0x00100000;}
00274             uz.d = ax;
00275            __HI(uz) = ix;
00276             ax = uz.d;
00277 
00278        /* compute s = s_h+s_l = (x-1)/(x+1) or (x-1.5)/(x+1.5) */
00279            u = ax-bp[k];           /* bp[0]=1.0, bp[1]=1.5 */
00280            v = one/(ax+bp[k]);
00281            s = u*v;
00282            s_h = s;
00283             uz.d = s_h;
00284            __LO(uz) = 0;
00285             s_h = uz.d;
00286        /* t_h=ax+bp[k] High */
00287            t_h = zero;
00288             uz.d = t_h;
00289            __HI(uz)=((ix>>1)|0x20000000)+0x00080000+(k<<18); 
00290             t_h = uz.d;
00291            t_l = ax - (t_h-bp[k]);
00292            s_l = v*((u-s_h*t_h)-s_h*t_l);
00293        /* compute log(ax) */
00294            s2 = s*s;
00295            r = s2*s2*(L1+s2*(L2+s2*(L3+s2*(L4+s2*(L5+s2*L6)))));
00296            r += s_l*(s_h+s);
00297            s2  = s_h*s_h;
00298            t_h = 3.0+s2+r;
00299             uz.d = t_h;
00300            __LO(uz) = 0;
00301             t_h = uz.d;
00302            t_l = r-((t_h-3.0)-s2);
00303        /* u+v = s*(1+...) */
00304            u = s_h*t_h;
00305            v = s_l*t_h+t_l*s;
00306        /* 2/(3log2)*(s+...) */
00307            p_h = u+v;
00308             uz.d = p_h;
00309            __LO(uz) = 0;
00310             p_h = uz.d;
00311            p_l = v-(p_h-u);
00312            z_h = cp_h*p_h;         /* cp_h+cp_l = 2/(3*log2) */
00313            z_l = cp_l*p_h+p_l*cp+dp_l[k];
00314        /* log2(ax) = (s+..)*2/(3*log2) = n + dp_h + z_h + z_l */
00315            t = (double)n;
00316            t1 = (((z_h+z_l)+dp_h[k])+t);
00317             uz.d = t1;
00318            __LO(uz) = 0;
00319             t1 = uz.d;
00320            t2 = z_l-(((t1-t)-dp_h[k])-z_h);
00321        }
00322 
00323        s = one; /* s (sign of result -ve**odd) = -1 else = 1 */
00324        if((((hx>>31)+1)|(yisint-1))==0) s = -one;/* (-ve)**(odd int) */
00325 
00326     /* split up y into y1+y2 and compute (y1+y2)*(t1+t2) */
00327        y1  = y;
00328         uy.d = y1;
00329        __LO(uy) = 0;
00330         y1 = uy.d;
00331        p_l = (y-y1)*t1+y*t2;
00332        p_h = y1*t1;
00333        z = p_l+p_h;
00334         uz.d = z;
00335        j = __HI(uz);
00336        i = __LO(uz);
00337 
00338        if (j>=0x40900000) {                      /* z >= 1024 */
00339            if(((j-0x40900000)|i)!=0)                    /* if z > 1024 */
00340               return s*really_big*really_big;                  /* overflow */
00341            else {
00342               if(p_l+ovt>z-p_h) return s*really_big*really_big;       /* overflow */
00343            }
00344        } else if((j&0x7fffffff)>=0x4090cc00 ) {  /* z <= -1075 */
00345            if(((j-0xc090cc00)|i)!=0)             /* z < -1075 */
00346               return s*tiny*tiny;         /* underflow */
00347            else {
00348               if(p_l<=z-p_h) return s*tiny*tiny; /* underflow */
00349            }
00350        }
00351     /*
00352      * compute 2**(p_h+p_l)
00353      */
00354        i = j&0x7fffffff;
00355        k = (i>>20)-0x3ff;
00356        n = 0;
00357        if(i>0x3fe00000) {          /* if |z| > 0.5, set n = [z+0.5] */
00358            n = j+(0x00100000>>(k+1));
00359            k = ((n&0x7fffffff)>>20)-0x3ff;       /* new k for n */
00360            t = zero;
00361             uz.d = t;
00362            __HI(uz) = (n&~(0x000fffff>>k));
00363             t = uz.d;
00364            n = ((n&0x000fffff)|0x00100000)>>(20-k);
00365            if(j<0) n = -n;
00366            p_h -= t;
00367        } 
00368        t = p_l+p_h;
00369         uz.d = t;
00370        __LO(uz) = 0;
00371         t = uz.d;
00372        u = t*lg2_h;
00373        v = (p_l-(t-p_h))*lg2+t*lg2_l;
00374        z = u+v;
00375        w = v-(z-u);
00376        t  = z*z;
00377        t1  = z - t*(P1+t*(P2+t*(P3+t*(P4+t*P5))));
00378        r  = (z*t1)/(t1-two)-(w+z*w);
00379        z  = one-(r-z);
00380         uz.d = z;
00381        j  = __HI(uz);
00382        j += (n<<20);
00383        if((j>>20)<=0) z = fd_scalbn(z,n); /* subnormal output */
00384        else { uz.d = z; __HI(uz) += (n<<20); z = uz.d; }
00385        return s*z;
00386 }