From 1077fa4d772832f77a677ce7fb7c2d513b959e3f Mon Sep 17 00:00:00 2001 From: Eric Andersen Date: Thu, 10 May 2001 00:40:28 +0000 Subject: uClibc now has a math library. muahahahaha! -Erik --- libm/double/pow.c | 756 ++++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 756 insertions(+) create mode 100644 libm/double/pow.c (limited to 'libm/double/pow.c') diff --git a/libm/double/pow.c b/libm/double/pow.c new file mode 100644 index 000000000..768ad1062 --- /dev/null +++ b/libm/double/pow.c @@ -0,0 +1,756 @@ +/* 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 + * DEC -26,26 60000 4.8e-17 9.1e-18 + * 1/26 < x < 26, with log(x) uniformly distributed. + * -26 < y < 26, y uniformly distributed. + * 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 +*/ + + +#include +static char fname[] = {"pow"}; + +#define SQRTH 0.70710678118654752440 + +#ifdef UNK +static double P[] = { + 4.97778295871696322025E-1, + 3.73336776063286838734E0, + 7.69994162726912503298E0, + 4.66651806774358464979E0 +}; +static double Q[] = { +/* 1.00000000000000000000E0, */ + 9.33340916416696166113E0, + 2.79999886606328401649E1, + 3.35994905342304405431E1, + 1.39995542032307539578E1 +}; +/* 2^(-i/16), IEEE precision */ +static double A[] = { + 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 double B[] = { + 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 double R[] = { + 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] +#define doubb(k) B[k] +#define MEXP 16383.0 +#ifdef DENORMAL +#define MNEXP -17183.0 +#else +#define MNEXP -16383.0 +#endif +#endif + +#ifdef DEC +static unsigned short P[] = { +0037776,0156313,0175332,0163602, +0040556,0167577,0052366,0174245, +0040766,0062753,0175707,0055564, +0040625,0052035,0131344,0155636, +}; +static unsigned short Q[] = { +/*0040200,0000000,0000000,0000000,*/ +0041025,0052644,0154404,0105155, +0041337,0177772,0007016,0047646, +0041406,0062740,0154273,0020020, +0041137,0177054,0106127,0044555, +}; +static unsigned short A[] = { +0040200,0000000,0000000,0000000, +0040165,0022575,0012444,0103314, +0040152,0140306,0163735,0022071, +0040140,0146336,0166052,0112341, +0040127,0042374,0145326,0116553, +0040116,0022214,0012437,0102201, +0040105,0063452,0010525,0003333, +0040075,0004243,0117530,0006067, +0040065,0002363,0031771,0157145, +0040055,0054076,0165102,0120513, +0040045,0177326,0124661,0050471, +0040036,0172462,0060221,0120422, +0040030,0033760,0050615,0134251, +0040021,0141723,0071653,0010703, +0040013,0112701,0161752,0105727, +0040005,0125303,0063714,0044173, +0040000,0000000,0000000,0000000 +}; +static unsigned short B[] = { +0000000,0000000,0000000,0000000, +0021473,0040265,0153315,0140671, +0121074,0062627,0042146,0176454, +0121413,0003524,0136332,0066212, +0121767,0046404,0166231,0012553, +0121257,0015024,0002357,0043574, +0021736,0106532,0043060,0056206, +0121310,0020334,0165705,0035326, +0000000,0000000,0000000,0000000 +}; + +static unsigned short R[] = { +0034173,0014076,0137624,0115771, +0035041,0076763,0003744,0111311, +0035656,0141766,0041127,0074351, +0036435,0112533,0073611,0116664, +0037143,0054106,0134040,0152223, +0037565,0176757,0176026,0025551, +0040061,0071027,0173721,0147572 +}; + +/* +static double R[] = { +0.14928852680595608186e-4, +0.15400290440989764601e-3, +0.13333541313585784703e-2, +0.96181290595172416964e-2, +0.55504108664085595326e-1, +0.24022650695909537056e0, +0.69314718055994529629e0 +}; +*/ +#define douba(k) (*(double *)&A[(k)<<2]) +#define doubb(k) (*(double *)&B[(k)<<2]) +#define MEXP 2031.0 +#define MNEXP -2031.0 +#endif + +#ifdef IBMPC +static unsigned short P[] = { +0x5cf0,0x7f5b,0xdb99,0x3fdf, +0xdf15,0xea9e,0xddef,0x400d, +0xeb6f,0x7f78,0xccbd,0x401e, +0x9b74,0xb65c,0xaa83,0x4012, +}; +static unsigned short Q[] = { +/*0x0000,0x0000,0x0000,0x3ff0,*/ +0x914e,0x9b20,0xaab4,0x4022, +0xc9f5,0x41c1,0xffff,0x403b, +0x6402,0x1b17,0xccbc,0x4040, +0xe92e,0x918a,0xffc5,0x402b, +}; +static unsigned short A[] = { +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 unsigned short B[] = { +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 unsigned short R[] = { +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) (*(double *)&A[(k)<<2]) +#define doubb(k) (*(double *)&B[(k)<<2]) +#define MEXP 16383.0 +#ifdef DENORMAL +#define MNEXP -17183.0 +#else +#define MNEXP -16383.0 +#endif +#endif + +#ifdef MIEEE +static unsigned short P[] = { +0x3fdf,0xdb99,0x7f5b,0x5cf0, +0x400d,0xddef,0xea9e,0xdf15, +0x401e,0xccbd,0x7f78,0xeb6f, +0x4012,0xaa83,0xb65c,0x9b74 +}; +static unsigned short Q[] = { +0x4022,0xaab4,0x9b20,0x914e, +0x403b,0xffff,0x41c1,0xc9f5, +0x4040,0xccbc,0x1b17,0x6402, +0x402b,0xffc5,0x918a,0xe92e +}; +static unsigned short A[] = { +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 unsigned short B[] = { +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 unsigned short R[] = { +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) (*(double *)&A[(k)<<2]) +#define doubb(k) (*(double *)&B[(k)<<2]) +#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 ANSIPROT +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, void *, int ); +extern double powi ( double, int ); +extern int signbit ( double ); +extern int isnan ( double ); +extern int isfinite ( double ); +static double reduc ( double ); +#else +double floor(), fabs(), frexp(), ldexp(); +double polevl(), p1evl(), powi(); +int signbit(), isnan(), isfinite(); +static double reduc(); +#endif +extern double MAXNUM; +#ifdef INFINITIES +extern double INFINITY; +#endif +#ifdef NANS +extern double NAN; +#endif +#ifdef MINUSZERO +extern double NEGZERO; +#endif + +double pow( x, y ) +double x, 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) ) + return( x ); +if( isnan(y) ) + return( 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 ) + { +#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 ) + { + 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 ); +#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(x) +double x; +{ +double t; + +t = ldexp( x, 4 ); +t = floor( t ); +t = ldexp( t, -4 ); +return(t); +} -- cgit v1.2.3