src/libm/e_pow.c
changeset 2756 a98604b691c8
child 3337 9ac6f0782dd6
     1.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     1.2 +++ b/src/libm/e_pow.c	Mon Sep 15 06:33:23 2008 +0000
     1.3 @@ -0,0 +1,342 @@
     1.4 +/* @(#)e_pow.c 5.1 93/09/24 */
     1.5 +/*
     1.6 + * ====================================================
     1.7 + * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
     1.8 + *
     1.9 + * Developed at SunPro, a Sun Microsystems, Inc. business.
    1.10 + * Permission to use, copy, modify, and distribute this
    1.11 + * software is freely granted, provided that this notice
    1.12 + * is preserved.
    1.13 + * ====================================================
    1.14 + */
    1.15 +
    1.16 +#if defined(LIBM_SCCS) && !defined(lint)
    1.17 +static char rcsid[] = "$NetBSD: e_pow.c,v 1.9 1995/05/12 04:57:32 jtc Exp $";
    1.18 +#endif
    1.19 +
    1.20 +/* __ieee754_pow(x,y) return x**y
    1.21 + *
    1.22 + *		      n
    1.23 + * Method:  Let x =  2   * (1+f)
    1.24 + *	1. Compute and return log2(x) in two pieces:
    1.25 + *		log2(x) = w1 + w2,
    1.26 + *	   where w1 has 53-24 = 29 bit trailing zeros.
    1.27 + *	2. Perform y*log2(x) = n+y' by simulating muti-precision
    1.28 + *	   arithmetic, where |y'|<=0.5.
    1.29 + *	3. Return x**y = 2**n*exp(y'*log2)
    1.30 + *
    1.31 + * Special cases:
    1.32 + *	1.  (anything) ** 0  is 1
    1.33 + *	2.  (anything) ** 1  is itself
    1.34 + *	3.  (anything) ** NAN is NAN
    1.35 + *	4.  NAN ** (anything except 0) is NAN
    1.36 + *	5.  +-(|x| > 1) **  +INF is +INF
    1.37 + *	6.  +-(|x| > 1) **  -INF is +0
    1.38 + *	7.  +-(|x| < 1) **  +INF is +0
    1.39 + *	8.  +-(|x| < 1) **  -INF is +INF
    1.40 + *	9.  +-1         ** +-INF is NAN
    1.41 + *	10. +0 ** (+anything except 0, NAN)               is +0
    1.42 + *	11. -0 ** (+anything except 0, NAN, odd integer)  is +0
    1.43 + *	12. +0 ** (-anything except 0, NAN)               is +INF
    1.44 + *	13. -0 ** (-anything except 0, NAN, odd integer)  is +INF
    1.45 + *	14. -0 ** (odd integer) = -( +0 ** (odd integer) )
    1.46 + *	15. +INF ** (+anything except 0,NAN) is +INF
    1.47 + *	16. +INF ** (-anything except 0,NAN) is +0
    1.48 + *	17. -INF ** (anything)  = -0 ** (-anything)
    1.49 + *	18. (-anything) ** (integer) is (-1)**(integer)*(+anything**integer)
    1.50 + *	19. (-anything except 0 and inf) ** (non-integer) is NAN
    1.51 + *
    1.52 + * Accuracy:
    1.53 + *	pow(x,y) returns x**y nearly rounded. In particular
    1.54 + *			pow(integer,integer)
    1.55 + *	always returns the correct integer provided it is
    1.56 + *	representable.
    1.57 + *
    1.58 + * Constants :
    1.59 + * The hexadecimal values are the intended ones for the following
    1.60 + * constants. The decimal values may be used, provided that the
    1.61 + * compiler will convert from decimal to binary accurately enough
    1.62 + * to produce the hexadecimal values shown.
    1.63 + */
    1.64 +
    1.65 +#include "math.h"
    1.66 +#include "math_private.h"
    1.67 +
    1.68 +libm_hidden_proto(scalbn)
    1.69 +    libm_hidden_proto(fabs)
    1.70 +#ifdef __STDC__
    1.71 +     static const double
    1.72 +#else
    1.73 +     static double
    1.74 +#endif
    1.75 +       bp[] = { 1.0, 1.5, }, dp_h[] = {
    1.76 +     0.0, 5.84962487220764160156e-01,}, /* 0x3FE2B803, 0x40000000 */
    1.77 +
    1.78 +         dp_l[] = {
    1.79 +     0.0, 1.35003920212974897128e-08,}, /* 0x3E4CFDEB, 0x43CFD006 */
    1.80 +
    1.81 +         zero = 0.0, one = 1.0, two = 2.0, two53 = 9007199254740992.0,  /* 0x43400000, 0x00000000 */
    1.82 +         huge = 1.0e300, tiny = 1.0e-300,
    1.83 +         /* poly coefs for (3/2)*(log(x)-2s-2/3*s**3 */
    1.84 +         L1 = 5.99999999999994648725e-01,       /* 0x3FE33333, 0x33333303 */
    1.85 +         L2 = 4.28571428578550184252e-01,       /* 0x3FDB6DB6, 0xDB6FABFF */
    1.86 +         L3 = 3.33333329818377432918e-01,       /* 0x3FD55555, 0x518F264D */
    1.87 +         L4 = 2.72728123808534006489e-01,       /* 0x3FD17460, 0xA91D4101 */
    1.88 +         L5 = 2.30660745775561754067e-01,       /* 0x3FCD864A, 0x93C9DB65 */
    1.89 +         L6 = 2.06975017800338417784e-01,       /* 0x3FCA7E28, 0x4A454EEF */
    1.90 +         P1 = 1.66666666666666019037e-01,       /* 0x3FC55555, 0x5555553E */
    1.91 +         P2 = -2.77777777770155933842e-03,      /* 0xBF66C16C, 0x16BEBD93 */
    1.92 +         P3 = 6.61375632143793436117e-05,       /* 0x3F11566A, 0xAF25DE2C */
    1.93 +         P4 = -1.65339022054652515390e-06,      /* 0xBEBBBD41, 0xC5D26BF1 */
    1.94 +         P5 = 4.13813679705723846039e-08,       /* 0x3E663769, 0x72BEA4D0 */
    1.95 +         lg2 = 6.93147180559945286227e-01,      /* 0x3FE62E42, 0xFEFA39EF */
    1.96 +         lg2_h = 6.93147182464599609375e-01,    /* 0x3FE62E43, 0x00000000 */
    1.97 +         lg2_l = -1.90465429995776804525e-09,   /* 0xBE205C61, 0x0CA86C39 */
    1.98 +         ovt = 8.0085662595372944372e-0017,     /* -(1024-log2(ovfl+.5ulp)) */
    1.99 +         cp = 9.61796693925975554329e-01,       /* 0x3FEEC709, 0xDC3A03FD =2/(3ln2) */
   1.100 +         cp_h = 9.61796700954437255859e-01,     /* 0x3FEEC709, 0xE0000000 =(float)cp */
   1.101 +         cp_l = -7.02846165095275826516e-09,    /* 0xBE3E2FE0, 0x145B01F5 =tail of cp_h */
   1.102 +         ivln2 = 1.44269504088896338700e+00,    /* 0x3FF71547, 0x652B82FE =1/ln2 */
   1.103 +         ivln2_h = 1.44269502162933349609e+00,  /* 0x3FF71547, 0x60000000 =24b 1/ln2 */
   1.104 +         ivln2_l = 1.92596299112661746887e-08;  /* 0x3E54AE0B, 0xF85DDF44 =1/ln2 tail */
   1.105 +
   1.106 +#ifdef __STDC__
   1.107 +     double attribute_hidden __ieee754_pow(double x, double y)
   1.108 +#else
   1.109 +     double attribute_hidden __ieee754_pow(x, y)
   1.110 +     double x, y;
   1.111 +#endif
   1.112 +     {
   1.113 +         double z, ax, z_h, z_l, p_h, p_l;
   1.114 +         double y1, t1, t2, r, s, t, u, v, w;
   1.115 +         int32_t i, j, k, yisint, n;
   1.116 +         int32_t hx, hy, ix, iy;
   1.117 +         u_int32_t lx, ly;
   1.118 +
   1.119 +         EXTRACT_WORDS(hx, lx, x);
   1.120 +         EXTRACT_WORDS(hy, ly, y);
   1.121 +         ix = hx & 0x7fffffff;
   1.122 +         iy = hy & 0x7fffffff;
   1.123 +
   1.124 +         /* y==zero: x**0 = 1 */
   1.125 +         if ((iy | ly) == 0)
   1.126 +             return one;
   1.127 +
   1.128 +         /* +-NaN return x+y */
   1.129 +         if (ix > 0x7ff00000 || ((ix == 0x7ff00000) && (lx != 0)) ||
   1.130 +             iy > 0x7ff00000 || ((iy == 0x7ff00000) && (ly != 0)))
   1.131 +             return x + y;
   1.132 +
   1.133 +         /* determine if y is an odd int when x < 0
   1.134 +          * yisint = 0       ... y is not an integer
   1.135 +          * yisint = 1       ... y is an odd int
   1.136 +          * yisint = 2       ... y is an even int
   1.137 +          */
   1.138 +         yisint = 0;
   1.139 +         if (hx < 0) {
   1.140 +             if (iy >= 0x43400000)
   1.141 +                 yisint = 2;    /* even integer y */
   1.142 +             else if (iy >= 0x3ff00000) {
   1.143 +                 k = (iy >> 20) - 0x3ff;        /* exponent */
   1.144 +                 if (k > 20) {
   1.145 +                     j = ly >> (52 - k);
   1.146 +                     if ((j << (52 - k)) == ly)
   1.147 +                         yisint = 2 - (j & 1);
   1.148 +                 } else if (ly == 0) {
   1.149 +                     j = iy >> (20 - k);
   1.150 +                     if ((j << (20 - k)) == iy)
   1.151 +                         yisint = 2 - (j & 1);
   1.152 +                 }
   1.153 +             }
   1.154 +         }
   1.155 +
   1.156 +         /* special value of y */
   1.157 +         if (ly == 0) {
   1.158 +             if (iy == 0x7ff00000) {    /* y is +-inf */
   1.159 +                 if (((ix - 0x3ff00000) | lx) == 0)
   1.160 +                     return y - y;      /* inf**+-1 is NaN */
   1.161 +                 else if (ix >= 0x3ff00000)     /* (|x|>1)**+-inf = inf,0 */
   1.162 +                     return (hy >= 0) ? y : zero;
   1.163 +                 else           /* (|x|<1)**-,+inf = inf,0 */
   1.164 +                     return (hy < 0) ? -y : zero;
   1.165 +             }
   1.166 +             if (iy == 0x3ff00000) {    /* y is  +-1 */
   1.167 +                 if (hy < 0)
   1.168 +                     return one / x;
   1.169 +                 else
   1.170 +                     return x;
   1.171 +             }
   1.172 +             if (hy == 0x40000000)
   1.173 +                 return x * x;  /* y is  2 */
   1.174 +             if (hy == 0x3fe00000) {    /* y is  0.5 */
   1.175 +                 if (hx >= 0)   /* x >= +0 */
   1.176 +                     return __ieee754_sqrt(x);
   1.177 +             }
   1.178 +         }
   1.179 +
   1.180 +         ax = fabs(x);
   1.181 +         /* special value of x */
   1.182 +         if (lx == 0) {
   1.183 +             if (ix == 0x7ff00000 || ix == 0 || ix == 0x3ff00000) {
   1.184 +                 z = ax;        /*x is +-0,+-inf,+-1 */
   1.185 +                 if (hy < 0)
   1.186 +                     z = one / z;       /* z = (1/|x|) */
   1.187 +                 if (hx < 0) {
   1.188 +                     if (((ix - 0x3ff00000) | yisint) == 0) {
   1.189 +                         z = (z - z) / (z - z); /* (-1)**non-int is NaN */
   1.190 +                     } else if (yisint == 1)
   1.191 +                         z = -z;        /* (x<0)**odd = -(|x|**odd) */
   1.192 +                 }
   1.193 +                 return z;
   1.194 +             }
   1.195 +         }
   1.196 +
   1.197 +         /* (x<0)**(non-int) is NaN */
   1.198 +         if (((((u_int32_t) hx >> 31) - 1) | yisint) == 0)
   1.199 +             return (x - x) / (x - x);
   1.200 +
   1.201 +         /* |y| is huge */
   1.202 +         if (iy > 0x41e00000) { /* if |y| > 2**31 */
   1.203 +             if (iy > 0x43f00000) {     /* if |y| > 2**64, must o/uflow */
   1.204 +                 if (ix <= 0x3fefffff)
   1.205 +                     return (hy < 0) ? huge * huge : tiny * tiny;
   1.206 +                 if (ix >= 0x3ff00000)
   1.207 +                     return (hy > 0) ? huge * huge : tiny * tiny;
   1.208 +             }
   1.209 +             /* over/underflow if x is not close to one */
   1.210 +             if (ix < 0x3fefffff)
   1.211 +                 return (hy < 0) ? huge * huge : tiny * tiny;
   1.212 +             if (ix > 0x3ff00000)
   1.213 +                 return (hy > 0) ? huge * huge : tiny * tiny;
   1.214 +             /* now |1-x| is tiny <= 2**-20, suffice to compute
   1.215 +                log(x) by x-x^2/2+x^3/3-x^4/4 */
   1.216 +             t = x - 1;         /* t has 20 trailing zeros */
   1.217 +             w = (t * t) * (0.5 - t * (0.3333333333333333333333 - t * 0.25));
   1.218 +             u = ivln2_h * t;   /* ivln2_h has 21 sig. bits */
   1.219 +             v = t * ivln2_l - w * ivln2;
   1.220 +             t1 = u + v;
   1.221 +             SET_LOW_WORD(t1, 0);
   1.222 +             t2 = v - (t1 - u);
   1.223 +         } else {
   1.224 +             double s2, s_h, s_l, t_h, t_l;
   1.225 +             n = 0;
   1.226 +             /* take care subnormal number */
   1.227 +             if (ix < 0x00100000) {
   1.228 +                 ax *= two53;
   1.229 +                 n -= 53;
   1.230 +                 GET_HIGH_WORD(ix, ax);
   1.231 +             }
   1.232 +             n += ((ix) >> 20) - 0x3ff;
   1.233 +             j = ix & 0x000fffff;
   1.234 +             /* determine interval */
   1.235 +             ix = j | 0x3ff00000;       /* normalize ix */
   1.236 +             if (j <= 0x3988E)
   1.237 +                 k = 0;         /* |x|<sqrt(3/2) */
   1.238 +             else if (j < 0xBB67A)
   1.239 +                 k = 1;         /* |x|<sqrt(3)   */
   1.240 +             else {
   1.241 +                 k = 0;
   1.242 +                 n += 1;
   1.243 +                 ix -= 0x00100000;
   1.244 +             }
   1.245 +             SET_HIGH_WORD(ax, ix);
   1.246 +
   1.247 +             /* compute s = s_h+s_l = (x-1)/(x+1) or (x-1.5)/(x+1.5) */
   1.248 +             u = ax - bp[k];    /* bp[0]=1.0, bp[1]=1.5 */
   1.249 +             v = one / (ax + bp[k]);
   1.250 +             s = u * v;
   1.251 +             s_h = s;
   1.252 +             SET_LOW_WORD(s_h, 0);
   1.253 +             /* t_h=ax+bp[k] High */
   1.254 +             t_h = zero;
   1.255 +             SET_HIGH_WORD(t_h,
   1.256 +                           ((ix >> 1) | 0x20000000) + 0x00080000 + (k << 18));
   1.257 +             t_l = ax - (t_h - bp[k]);
   1.258 +             s_l = v * ((u - s_h * t_h) - s_h * t_l);
   1.259 +             /* compute log(ax) */
   1.260 +             s2 = s * s;
   1.261 +             r = s2 * s2 * (L1 +
   1.262 +                            s2 * (L2 +
   1.263 +                                  s2 * (L3 +
   1.264 +                                        s2 * (L4 + s2 * (L5 + s2 * L6)))));
   1.265 +             r += s_l * (s_h + s);
   1.266 +             s2 = s_h * s_h;
   1.267 +             t_h = 3.0 + s2 + r;
   1.268 +             SET_LOW_WORD(t_h, 0);
   1.269 +             t_l = r - ((t_h - 3.0) - s2);
   1.270 +             /* u+v = s*(1+...) */
   1.271 +             u = s_h * t_h;
   1.272 +             v = s_l * t_h + t_l * s;
   1.273 +             /* 2/(3log2)*(s+...) */
   1.274 +             p_h = u + v;
   1.275 +             SET_LOW_WORD(p_h, 0);
   1.276 +             p_l = v - (p_h - u);
   1.277 +             z_h = cp_h * p_h;  /* cp_h+cp_l = 2/(3*log2) */
   1.278 +             z_l = cp_l * p_h + p_l * cp + dp_l[k];
   1.279 +             /* log2(ax) = (s+..)*2/(3*log2) = n + dp_h + z_h + z_l */
   1.280 +             t = (double) n;
   1.281 +             t1 = (((z_h + z_l) + dp_h[k]) + t);
   1.282 +             SET_LOW_WORD(t1, 0);
   1.283 +             t2 = z_l - (((t1 - t) - dp_h[k]) - z_h);
   1.284 +         }
   1.285 +
   1.286 +         s = one;               /* s (sign of result -ve**odd) = -1 else = 1 */
   1.287 +         if (((((u_int32_t) hx >> 31) - 1) | (yisint - 1)) == 0)
   1.288 +             s = -one;          /* (-ve)**(odd int) */
   1.289 +
   1.290 +         /* split up y into y1+y2 and compute (y1+y2)*(t1+t2) */
   1.291 +         y1 = y;
   1.292 +         SET_LOW_WORD(y1, 0);
   1.293 +         p_l = (y - y1) * t1 + y * t2;
   1.294 +         p_h = y1 * t1;
   1.295 +         z = p_l + p_h;
   1.296 +         EXTRACT_WORDS(j, i, z);
   1.297 +         if (j >= 0x40900000) { /* z >= 1024 */
   1.298 +             if (((j - 0x40900000) | i) != 0)   /* if z > 1024 */
   1.299 +                 return s * huge * huge;        /* overflow */
   1.300 +             else {
   1.301 +                 if (p_l + ovt > z - p_h)
   1.302 +                     return s * huge * huge;    /* overflow */
   1.303 +             }
   1.304 +         } else if ((j & 0x7fffffff) >= 0x4090cc00) {   /* z <= -1075 */
   1.305 +             if (((j - 0xc090cc00) | i) != 0)   /* z < -1075 */
   1.306 +                 return s * tiny * tiny;        /* underflow */
   1.307 +             else {
   1.308 +                 if (p_l <= z - p_h)
   1.309 +                     return s * tiny * tiny;    /* underflow */
   1.310 +             }
   1.311 +         }
   1.312 +         /*
   1.313 +          * compute 2**(p_h+p_l)
   1.314 +          */
   1.315 +         i = j & 0x7fffffff;
   1.316 +         k = (i >> 20) - 0x3ff;
   1.317 +         n = 0;
   1.318 +         if (i > 0x3fe00000) {  /* if |z| > 0.5, set n = [z+0.5] */
   1.319 +             n = j + (0x00100000 >> (k + 1));
   1.320 +             k = ((n & 0x7fffffff) >> 20) - 0x3ff;      /* new k for n */
   1.321 +             t = zero;
   1.322 +             SET_HIGH_WORD(t, n & ~(0x000fffff >> k));
   1.323 +             n = ((n & 0x000fffff) | 0x00100000) >> (20 - k);
   1.324 +             if (j < 0)
   1.325 +                 n = -n;
   1.326 +             p_h -= t;
   1.327 +         }
   1.328 +         t = p_l + p_h;
   1.329 +         SET_LOW_WORD(t, 0);
   1.330 +         u = t * lg2_h;
   1.331 +         v = (p_l - (t - p_h)) * lg2 + t * lg2_l;
   1.332 +         z = u + v;
   1.333 +         w = v - (z - u);
   1.334 +         t = z * z;
   1.335 +         t1 = z - t * (P1 + t * (P2 + t * (P3 + t * (P4 + t * P5))));
   1.336 +         r = (z * t1) / (t1 - two) - (w + z * w);
   1.337 +         z = one - (r - z);
   1.338 +         GET_HIGH_WORD(j, z);
   1.339 +         j += (n << 20);
   1.340 +         if ((j >> 20) <= 0)
   1.341 +             z = scalbn(z, n);  /* subnormal output */
   1.342 +         else
   1.343 +             SET_HIGH_WORD(z, j);
   1.344 +         return s * z;
   1.345 +     }