| /* pow.c |
| * |
| * Power function |
| * |
| * |
| * |
| * SYNOPSIS: |
| * |
| * double x, y, z, pow(); |
| * |
| * z = pow( x, y ); |
| * |
| * |
| * |
| * DESCRIPTION: |
| * |
| * Computes x raised to the yth power. Analytically, |
| * |
| * x**y = exp( y log(x) ). |
| * |
| * Following Cody and Waite, this program uses a lookup table |
| * of 2**-i/16 and pseudo extended precision arithmetic to |
| * obtain an extra three bits of accuracy in both the logarithm |
| * and the exponential. |
| * |
| * |
| * |
| * ACCURACY: |
| * |
| * Relative error: |
| * arithmetic domain # trials peak rms |
| * IEEE -26,26 30000 4.2e-16 7.7e-17 |
| * IEEE 0,8700 30000 1.5e-14 2.1e-15 |
| * 0.99 < x < 1.01, 0 < y < 8700, uniformly distributed. |
| * |
| * |
| * ERROR MESSAGES: |
| * |
| * message condition value returned |
| * pow overflow x**y > MAXNUM INFINITY |
| * pow underflow x**y < 1/MAXNUM 0.0 |
| * pow domain x<0 and y noninteger 0.0 |
| * |
| */ |
| |
| /* |
| Cephes Math Library Release 2.8: June, 2000 |
| Copyright 1984, 1995, 2000 by Stephen L. Moshier |
| */ |
| |
| /* |
| Modified for mingw |
| 2002-09-27 Danny Smith <dannysmith@users.sourceforge.net> |
| */ |
| |
| #ifdef __MINGW32__ |
| #include "cephes_mconf.h" |
| #else |
| #include "mconf.h" |
| static char fname[] = {"pow"}; |
| #endif |
| |
| #ifndef _SET_ERRNO |
| #define _SET_ERRNO(x) |
| #endif |
| |
| #define SQRTH 0.70710678118654752440 |
| |
| #ifdef UNK |
| static uD P[4] = { |
| { { 4.97778295871696322025E-1 } }, |
| { { 3.73336776063286838734E0 } }, |
| { { 7.69994162726912503298E0 } }, |
| { { 4.66651806774358464979E0 } } |
| }; |
| static uD Q[4] = { |
| { { 9.33340916416696166113E0 } }, |
| { { 2.79999886606328401649E1 } }, |
| { { 3.35994905342304405431E1 } }, |
| { { 1.39995542032307539578E1 } } |
| }; |
| /* 2^(-i/16), IEEE precision */ |
| static uD A[17] = { |
| { { 1.00000000000000000000E0 } }, |
| { { 9.57603280698573700036E-1 } }, |
| { { 9.17004043204671215328E-1 } }, |
| { { 8.78126080186649726755E-1 } }, |
| { { 8.40896415253714502036E-1 } }, |
| { { 8.05245165974627141736E-1 } }, |
| { { 7.71105412703970372057E-1 } }, |
| { { 7.38413072969749673113E-1 } }, |
| { { 7.07106781186547572737E-1 } }, |
| { { 6.77127773468446325644E-1 } }, |
| { { 6.48419777325504820276E-1 } }, |
| { { 6.20928906036742001007E-1 } }, |
| { { 5.94603557501360513449E-1 } }, |
| { { 5.69394317378345782288E-1 } }, |
| { { 5.45253866332628844837E-1 } }, |
| { { 5.22136891213706877402E-1 } }, |
| { { 5.00000000000000000000E-1 } } |
| }; |
| static uD B[9] = { |
| { { 0.00000000000000000000E0 } }, |
| { { 1.64155361212281360176E-17 } }, |
| { { 4.09950501029074826006E-17 } }, |
| { { 3.97491740484881042808E-17 } }, |
| { { -4.83364665672645672553E-17 } }, |
| { { 1.26912513974441574796E-17 } }, |
| { { 1.99100761573282305549E-17 } }, |
| { { -1.52339103990623557348E-17 } }, |
| { { 0.00000000000000000000E0 } } |
| }; |
| static uD R[7] = { |
| { { 1.49664108433729301083E-5 } }, |
| { { 1.54010762792771901396E-4 } }, |
| { { 1.33335476964097721140E-3 } }, |
| { { 9.61812908476554225149E-3 } }, |
| { { 5.55041086645832347466E-2 } }, |
| { { 2.40226506959099779976E-1 } }, |
| { { 6.93147180559945308821E-1 } } |
| }; |
| |
| #define douba(k) A[k].d |
| #define doubb(k) B[k].d |
| #define MEXP 16383.0 |
| #ifdef DENORMAL |
| #define MNEXP -17183.0 |
| #else |
| #define MNEXP -16383.0 |
| #endif |
| #endif |
| |
| #ifdef IBMPC |
| static const uD P[4] = { |
| { { 0x5cf0,0x7f5b,0xdb99,0x3fdf } }, |
| { { 0xdf15,0xea9e,0xddef,0x400d } }, |
| { { 0xeb6f,0x7f78,0xccbd,0x401e } }, |
| { { 0x9b74,0xb65c,0xaa83,0x4012 } } |
| }; |
| static const uD Q[4] = { |
| { { 0x914e,0x9b20,0xaab4,0x4022 } }, |
| { { 0xc9f5,0x41c1,0xffff,0x403b } }, |
| { { 0x6402,0x1b17,0xccbc,0x4040 } }, |
| { { 0xe92e,0x918a,0xffc5,0x402b } } |
| }; |
| static const uD A[17] = { |
| { { 0x0000,0x0000,0x0000,0x3ff0 } }, |
| { { 0x90da,0xa2a4,0xa4af,0x3fee } }, |
| { { 0xa487,0xdcfb,0x5818,0x3fed } }, |
| { { 0x529c,0xdd85,0x199b,0x3fec } }, |
| { { 0xd3ad,0x995a,0xe89f,0x3fea } }, |
| { { 0xf090,0x82a3,0xc491,0x3fe9 } }, |
| { { 0xa0db,0x422a,0xace5,0x3fe8 } }, |
| { { 0x0187,0x73eb,0xa114,0x3fe7 } }, |
| { { 0x3bcd,0x667f,0xa09e,0x3fe6 } }, |
| { { 0x5429,0xdd48,0xab07,0x3fe5 } }, |
| { { 0x2a27,0xd536,0xbfda,0x3fe4 } }, |
| { { 0x3422,0x4c12,0xdea6,0x3fe3 } }, |
| { { 0xb715,0x0a31,0x06fe,0x3fe3 } }, |
| { { 0x6238,0x6e75,0x387a,0x3fe2 } }, |
| { { 0x517b,0x3c7d,0x72b8,0x3fe1 } }, |
| { { 0x890f,0x6cf9,0xb558,0x3fe0 } }, |
| { { 0x0000,0x0000,0x0000,0x3fe0 } } |
| }; |
| static const uD B[9] = { |
| { { 0x0000,0x0000,0x0000,0x0000 } }, |
| { { 0x3707,0xd75b,0xed02,0x3c72 } }, |
| { { 0xcc81,0x345d,0xa1cd,0x3c87 } }, |
| { { 0x4b27,0x5686,0xe9f1,0x3c86 } }, |
| { { 0x6456,0x13b2,0xdd34,0xbc8b } }, |
| { { 0x42e2,0xafec,0x4397,0x3c6d } }, |
| { { 0x82e4,0xd231,0xf46a,0x3c76 } }, |
| { { 0x8a76,0xb9d7,0x9041,0xbc71 } }, |
| { { 0x0000,0x0000,0x0000,0x0000 } } |
| }; |
| static const uD R[7] = { |
| { { 0x937f,0xd7f2,0x6307,0x3eef } }, |
| { { 0x9259,0x60fc,0x2fbe,0x3f24 } }, |
| { { 0xef1d,0xc84a,0xd87e,0x3f55 } }, |
| { { 0x33b7,0x6ef1,0xb2ab,0x3f83 } }, |
| { { 0x1a92,0xd704,0x6b08,0x3fac } }, |
| { { 0xc56d,0xff82,0xbfbd,0x3fce } }, |
| { { 0x39ef,0xfefa,0x2e42,0x3fe6 } } |
| }; |
| |
| #define douba(k) (A[(k)].d) |
| #define doubb(k) (B[(k)].d) |
| #define MEXP 16383.0 |
| #ifdef DENORMAL |
| #define MNEXP -17183.0 |
| #else |
| #define MNEXP -16383.0 |
| #endif |
| #endif |
| |
| #ifdef MIEEE |
| static uD P[4] = { |
| { { 0x3fdf,0xdb99,0x7f5b,0x5cf0 } }, |
| { { 0x400d,0xddef,0xea9e,0xdf15 } }, |
| { { 0x401e,0xccbd,0x7f78,0xeb6f } }, |
| { { 0x4012,0xaa83,0xb65c,0x9b74 } } |
| }; |
| static uD Q[4] = { |
| { { 0x4022,0xaab4,0x9b20,0x914e } }, |
| { { 0x403b,0xffff,0x41c1,0xc9f5 } }, |
| { { 0x4040,0xccbc,0x1b17,0x6402 } }, |
| { { 0x402b,0xffc5,0x918a,0xe92e } } |
| }; |
| static unsigned short A[17] = { |
| { { 0x3ff0,0x0000,0x0000,0x0000 } }, |
| { { 0x3fee,0xa4af,0xa2a4,0x90da } }, |
| { { 0x3fed,0x5818,0xdcfb,0xa487 } }, |
| { { 0x3fec,0x199b,0xdd85,0x529c } }, |
| { { 0x3fea,0xe89f,0x995a,0xd3ad } }, |
| { { 0x3fe9,0xc491,0x82a3,0xf090 } }, |
| { { 0x3fe8,0xace5,0x422a,0xa0db } }, |
| { { 0x3fe7,0xa114,0x73eb,0x0187 } }, |
| { { 0x3fe6,0xa09e,0x667f,0x3bcd } }, |
| { { 0x3fe5,0xab07,0xdd48,0x5429 } }, |
| { { 0x3fe4,0xbfda,0xd536,0x2a27 } }, |
| { { 0x3fe3,0xdea6,0x4c12,0x3422 } }, |
| { { 0x3fe3,0x06fe,0x0a31,0xb715 } }, |
| { { 0x3fe2,0x387a,0x6e75,0x6238 } }, |
| { { 0x3fe1,0x72b8,0x3c7d,0x517b } }, |
| { { 0x3fe0,0xb558,0x6cf9,0x890f } }, |
| { { 0x3fe0,0x0000,0x0000,0x0000 } } |
| }; |
| static uD B[9] = { |
| { { 0x0000,0x0000,0x0000,0x0000 } }, |
| { { 0x3c72,0xed02,0xd75b,0x3707 } }, |
| { { 0x3c87,0xa1cd,0x345d,0xcc81 } }, |
| { { 0x3c86,0xe9f1,0x5686,0x4b27 } }, |
| { { 0xbc8b,0xdd34,0x13b2,0x6456 } }, |
| { { 0x3c6d,0x4397,0xafec,0x42e2 } }, |
| { { 0x3c76,0xf46a,0xd231,0x82e4 } }, |
| { { 0xbc71,0x9041,0xb9d7,0x8a76 } }, |
| { { 0x0000,0x0000,0x0000,0x0000 } } |
| }; |
| static uD R[7] = { |
| { { 0x3eef,0x6307,0xd7f2,0x937f } }, |
| { { 0x3f24,0x2fbe,0x60fc,0x9259 } }, |
| { { 0x3f55,0xd87e,0xc84a,0xef1d } }, |
| { { 0x3f83,0xb2ab,0x6ef1,0x33b7 } }, |
| { { 0x3fac,0x6b08,0xd704,0x1a92 } }, |
| { { 0x3fce,0xbfbd,0xff82,0xc56d } }, |
| { { 0x3fe6,0x2e42,0xfefa,0x39ef } } |
| }; |
| |
| #define douba(k) (A[(k)].d) |
| #define doubb(k) (B[(k)].d) |
| #define MEXP 16383.0 |
| #ifdef DENORMAL |
| #define MNEXP -17183.0 |
| #else |
| #define MNEXP -16383.0 |
| #endif |
| #endif |
| |
| /* log2(e) - 1 */ |
| #define LOG2EA 0.44269504088896340736 |
| |
| #define F W |
| #define Fa Wa |
| #define Fb Wb |
| #define G W |
| #define Ga Wa |
| #define Gb u |
| #define H W |
| #define Ha Wb |
| #define Hb Wb |
| |
| #ifdef __MINGW32__ |
| static __inline__ double reduc( double ); |
| extern double __powi ( double, int ); |
| extern double pow ( double x, double y); |
| |
| #else /* __MINGW32__ */ |
| |
| extern double floor ( double ); |
| extern double fabs ( double ); |
| extern double frexp ( double, int * ); |
| extern double ldexp ( double, int ); |
| extern double polevl ( double, void *, int ); |
| extern double p1evl ( double, uD *, int ); |
| extern double __powi ( double, int ); |
| extern int signbit ( double ); |
| extern int isnan ( double ); |
| extern int isfinite ( double ); |
| static double reduc ( double ); |
| extern double MAXNUM; |
| #ifdef INFINITIES |
| extern double INFINITY; |
| #endif |
| #ifdef NANS |
| extern double NAN; |
| #endif |
| #ifdef MINUSZERO |
| extern double NEGZERO; |
| #endif |
| |
| #endif /* __MINGW32__ */ |
| |
| double pow(double x, double y) |
| { |
| double w, z, W, Wa, Wb, ya, yb, u; |
| /* double F, Fa, Fb, G, Ga, Gb, H, Ha, Hb */ |
| double aw, ay, wy; |
| int e, i, nflg, iyflg, yoddint; |
| |
| if (y == 0.0) |
| return (1.0); |
| #ifdef NANS |
| if (isnan(x) || isnan(y)) |
| { |
| _SET_ERRNO (EDOM); |
| return (x + y); |
| } |
| #endif |
| if (y == 1.0) |
| return (x); |
| |
| #ifdef INFINITIES |
| if (!isfinite(y) && (x == 1.0 || x == -1.0)) |
| { |
| mtherr( "pow", DOMAIN ); |
| #ifdef NANS |
| return( NAN ); |
| #else |
| return( INFINITY ); |
| #endif |
| } |
| #endif |
| |
| if (x == 1.0) |
| return (1.0); |
| |
| if (y >= MAXNUM) |
| { |
| _SET_ERRNO (ERANGE); |
| #ifdef INFINITIES |
| if (x > 1.0) |
| return (INFINITY); |
| #else |
| if (x > 1.0) |
| return (MAXNUM); |
| #endif |
| if (x > 0.0 && x < 1.0) |
| return (0.0); |
| if (x < -1.0) |
| { |
| #ifdef INFINITIES |
| return (INFINITY); |
| #else |
| return (MAXNUM); |
| #endif |
| } |
| if (x > -1.0 && x < 0.0) |
| return (0.0); |
| } |
| if (y <= -MAXNUM) |
| { |
| _SET_ERRNO (ERANGE); |
| if (x > 1.0) |
| return (0.0); |
| #ifdef INFINITIES |
| if (x > 0.0 && x < 1.0) |
| return (INFINITY); |
| #else |
| if (x > 0.0 && x < 1.0) |
| return (MAXNUM); |
| #endif |
| if (x < -1.0) |
| return (0.0); |
| #ifdef INFINITIES |
| if (x > -1.0 && x < 0.0) |
| return (INFINITY); |
| #else |
| if (x > -1.0 && x < 0.0) |
| return (MAXNUM); |
| #endif |
| } |
| if (x > MAXNUM) |
| { |
| #if INFINITIES |
| if (y > 0.0) |
| return (INFINITY); |
| #else |
| if (y > 0.0) |
| return (MAXNUM); |
| #endif |
| return (0.0); |
| } |
| /* Set iyflg to 1 if y is an integer. */ |
| iyflg = 0; |
| w = floor(y); |
| if (w == y) |
| iyflg = 1; |
| |
| /* Test for odd integer y. */ |
| yoddint = 0; |
| if (iyflg) |
| { |
| ya = fabs(y); |
| ya = floor(0.5 * ya); |
| yb = 0.5 * fabs(w); |
| if (ya != yb) |
| yoddint = 1; |
| } |
| |
| if (x < -MAXNUM) |
| { |
| if (y > 0.0) |
| { |
| #ifdef INFINITIES |
| if (yoddint) |
| return (-INFINITY); |
| return (INFINITY); |
| #else |
| if (yoddint) |
| return (-MAXNUM); |
| return (MAXNUM); |
| #endif |
| } |
| if (y < 0.0) |
| { |
| #ifdef MINUSZERO |
| if (yoddint) |
| return (NEGZERO); |
| #endif |
| return (0.0); |
| } |
| } |
| |
| nflg = 0; /* flag = 1 if x<0 raised to integer power */ |
| if (x <= 0.0) |
| { |
| if (x == 0.0) |
| { |
| if (y < 0.0) |
| { |
| #ifdef MINUSZERO |
| if (signbit(x) && yoddint) |
| return (-INFINITY); |
| #endif |
| #ifdef INFINITIES |
| return (INFINITY); |
| #else |
| return (MAXNUM); |
| #endif |
| } |
| if (y > 0.0) |
| { |
| #ifdef MINUSZERO |
| if (signbit(x) && yoddint) |
| return (NEGZERO); |
| #endif |
| return (0.0); |
| } |
| return (1.0); |
| } |
| else |
| { |
| if (iyflg == 0) |
| { /* noninteger power of negative number */ |
| mtherr(fname, DOMAIN); |
| _SET_ERRNO (EDOM); |
| #ifdef NANS |
| return (NAN); |
| #else |
| return (0.0L); |
| #endif |
| } |
| nflg = 1; |
| } |
| } |
| |
| /* Integer power of an integer. */ |
| |
| if (iyflg) |
| { |
| i = w; |
| w = floor(x); |
| if( (w == x) && (fabs(y) < 32768.0) ) |
| { |
| w = __powi(x, (int) y); |
| return (w); |
| } |
| } |
| |
| if (nflg) |
| x = fabs(x); |
| |
| /* For results close to 1, use a series expansion. */ |
| w = x - 1.0; |
| aw = fabs(w); |
| ay = fabs(y); |
| wy = w * y; |
| ya = fabs(wy); |
| if ((aw <= 1.0e-3 && ay <= 1.0) |
| || (ya <= 1.0e-3 && ay >= 1.0)) |
| { |
| z = (((((w*(y-5.)/720. + 1./120.)*w*(y-4.) + 1./24.)*w*(y-3.) |
| + 1./6.)*w*(y-2.) + 0.5)*w*(y-1.) )*wy + wy + 1.; |
| goto done; |
| } |
| /* These are probably too much trouble. */ |
| #if 0 |
| w = y * log(x); |
| if (aw > 1.0e-3 && fabs(w) < 1.0e-3) |
| { |
| z = ((((((w/7. + 1.)*w/6. + 1.)*w/5. + 1.)*w/4. + 1.)*w/3. + 1.)*w/2. + 1.)*w + 1.; |
| goto done; |
| } |
| |
| if (ya <= 1.0e-3 && aw <= 1.0e-4) |
| { |
| z = (((((wy*1./720. |
| + (-w*1./48. + 1./120.) )*wy |
| + ((w*17./144. - 1./12.)*w + 1./24.) )*wy |
| + (((-w*5./16. + 7./24.)*w - 1./4.)*w + 1./6.) )*wy |
| + ((((w*137./360. - 5./12.)*w + 11./24.)*w - 1./2.)*w + 1./2.) )*wy |
| + (((((-w*1./6. + 1./5.)*w - 1./4)*w + 1./3.)*w -1./2.)*w ) )*wy |
| + wy + 1.0; |
| goto done; |
| } |
| #endif |
| |
| /* separate significand from exponent */ |
| x = frexp(x, &e); |
| |
| #if 0 |
| /* For debugging, check for gross overflow. */ |
| if ((e * y) > (MEXP + 1024)) |
| goto overflow; |
| #endif |
| |
| /* Find significand of x in antilog table A[]. */ |
| i = 1; |
| if (x <= douba(9)) |
| i = 9; |
| if (x <= douba(i + 4)) |
| i += 4; |
| if (x <= douba(i + 2)) |
| i += 2; |
| if (x >= douba(1)) |
| i = -1; |
| i += 1; |
| |
| /* Find (x - A[i])/A[i] |
| * in order to compute log(x/A[i]): |
| * |
| * log(x) = log( a x/a ) = log(a) + log(x/a) |
| * |
| * log(x/a) = log(1+v), v = x/a - 1 = (x-a)/a |
| */ |
| x -= douba(i); |
| x -= doubb(i/2); |
| x /= douba(i); |
| |
| /* rational approximation for log(1+v): |
| * |
| * log(1+v) = v - v**2/2 + v**3 P(v) / Q(v) |
| */ |
| z = x*x; |
| w = x * ( z * polevl( x, P, 3 ) / p1evl( x, Q, 4 ) ); |
| w = w - ldexp( z, -1 ); /* w - 0.5 * z */ |
| |
| /* Convert to base 2 logarithm: |
| * multiply by log2(e) |
| */ |
| w = w + LOG2EA * w; |
| /* Note x was not yet added in |
| * to above rational approximation, |
| * so do it now, while multiplying |
| * by log2(e). |
| */ |
| z = w + LOG2EA * x; |
| z = z + x; |
| |
| /* Compute exponent term of the base 2 logarithm. */ |
| w = -i; |
| w = ldexp( w, -4 ); /* divide by 16 */ |
| w += e; |
| /* Now base 2 log of x is w + z. */ |
| |
| /* Multiply base 2 log by y, in extended precision. */ |
| |
| /* separate y into large part ya |
| * and small part yb less than 1/16 |
| */ |
| ya = reduc(y); |
| yb = y - ya; |
| |
| F = z * y + w * yb; |
| Fa = reduc(F); |
| Fb = F - Fa; |
| |
| G = Fa + w * ya; |
| Ga = reduc(G); |
| Gb = G - Ga; |
| |
| H = Fb + Gb; |
| Ha = reduc(H); |
| w = ldexp(Ga + Ha, 4); |
| |
| /* Test the power of 2 for overflow */ |
| if (w > MEXP) |
| { |
| #ifndef INFINITIES |
| mtherr(fname, OVERFLOW); |
| #endif |
| #ifdef INFINITIES |
| if (nflg && yoddint) |
| return (-INFINITY); |
| return (INFINITY); |
| #else |
| if (nflg && yoddint) |
| return (-MAXNUM); |
| return (MAXNUM); |
| #endif |
| } |
| |
| if (w < (MNEXP - 1)) |
| { |
| #ifndef DENORMAL |
| mtherr(fname, UNDERFLOW); |
| #endif |
| #ifdef MINUSZERO |
| if (nflg && yoddint) |
| return (NEGZERO); |
| #endif |
| return (0.0); |
| } |
| |
| e = w; |
| Hb = H - Ha; |
| |
| if (Hb > 0.0) |
| { |
| e += 1; |
| Hb -= 0.0625; |
| } |
| |
| /* Now the product y * log2(x) = Hb + e/16.0. |
| * |
| * Compute base 2 exponential of Hb, |
| * where -0.0625 <= Hb <= 0. |
| */ |
| z = Hb * polevl(Hb, R, 6); /* z = 2**Hb - 1 */ |
| |
| /* Express e/16 as an integer plus a negative number of 16ths. |
| * Find lookup table entry for the fractional power of 2. |
| */ |
| if (e < 0) |
| i = 0; |
| else |
| i = 1; |
| i = e/16 + i; |
| e = 16*i - e; |
| w = douba(e); |
| z = w + w * z; /* 2**-e * ( 1 + (2**Hb-1) ) */ |
| z = ldexp(z, i); /* multiply by integer power of 2 */ |
| |
| done: |
| /* Negate if odd integer power of negative number */ |
| if (nflg && yoddint) |
| { |
| #ifdef MINUSZERO |
| if (z == 0.0) |
| z = NEGZERO; |
| else |
| #endif |
| z = -z; |
| } |
| return (z); |
| } |
| |
| |
| /* Find a multiple of 1/16 that is within 1/16 of x. */ |
| static double reduc(double x) |
| { |
| double t; |
| |
| t = ldexp(x, 4); |
| t = floor(t); |
| t = ldexp(t, -4); |
| return (t); |
| } |
| |