Category : C Source Code
Archive   : CEPHES22.ZIP
Filename : POW.C

 
Output of file : POW.C contained in archive : CEPHES22.ZIP
/* 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 18000 4.4e-17 8.9e-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
* DEC 0,8700 1000 1.3e-15 2.1e-16
* 0.99 < x < 1.01, 0 < y < 8700, uniformly distributed.
*
*
* ERROR MESSAGES:
*
* message condition value returned
* pow overflow x**y > MAXNUM MAXNUM
* pow underflow x**y < 1/MAXNUM 0.0
* pow domain x<0 and y noninteger 0.0
*
*/

/*
Cephes Math Library Release 2.1: December, 1988
Copyright 1984, 1987, 1988 by Stephen L. Moshier
Direct inquiries to 30 Frost Street, Cambridge, MA 02140
*/


#include "mconf.h"
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 2031.0
#define MNEXP -2031.0
#endif

#ifdef DEC
static short P[] = {
0037776,0156313,0175332,0163602,
0040556,0167577,0052366,0174245,
0040766,0062753,0175707,0055564,
0040625,0052035,0131344,0155636,
};
static short Q[] = {
/*0040200,0000000,0000000,0000000,*/
0041025,0052644,0154404,0105155,
0041337,0177772,0007016,0047646,
0041406,0062740,0154273,0020020,
0041137,0177054,0106127,0044555,
};
static 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 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 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 short P[] = {
0x5cf0,0x7f5b,0xdb99,0x3fdf,
0xdf15,0xea9e,0xddef,0x400d,
0xeb6f,0x7f78,0xccbd,0x401e,
0x9b74,0xb65c,0xaa83,0x4012,
};
static short Q[] = {
/*0x0000,0x0000,0x0000,0x3ff0,*/
0x914e,0x9b20,0xaab4,0x4022,
0xc9f5,0x41c1,0xffff,0x403b,
0x6402,0x1b17,0xccbc,0x4040,
0xe92e,0x918a,0xffc5,0x402b,
};
static 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 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 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
/* The following if denormal numbers are supported, else -MEXP: */
#define MNEXP -17183.0
#endif

#ifdef MIEEE
static short P[] = {
0x3fdf,0xdb99,0x7f5b,0x5cf0,
0x400d,0xddef,0xea9e,0xdf15,
0x401e,0xccbd,0x7f78,0xeb6f,
0x4012,0xaa83,0xb65c,0x9b74
};
static short Q[] = {
0x4022,0xaab4,0x9b20,0x914e,
0x403b,0xffff,0x41c1,0xc9f5,
0x4040,0xccbc,0x1b17,0x6402,
0x402b,0xffc5,0x918a,0xe92e
};
static 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 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 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
#define MNEXP -17183.0
#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

extern double MAXNUM;

double pow( x, y )
double x, y;
{
double w, z, W, Wa, Wb, ya, yb;
double u, v;
/* double F, Fa, Fb, G, Ga, Gb, H, Ha, Hb */
int e, i, nflg;
double floor(), fabs(), frexp(), ldexp();
double reduc(), polevl(), p1evl(), powi();


nflg = 0; /* flag = 1 if x<0 raised to integer power */
w = floor(y);
if( (w == y) && (fabs(w) < 32768.0) )
{
i = w;
w = powi( x, i );
return( w );
}


if( x <= 0.0 )
{
if( x == 0.0 )
{
if( y == 0.0 )
return( 1.0 ); /* 0**0 */
else
return( 0.0 ); /* 0**y */
}
else
{
if( w != y )
{ /* noninteger power of negative number */
mtherr( fname, DOMAIN );
return(0.0);
}
nflg = 1;
x = fabs(x);
}
}

/* separate significand from exponent */
x = frexp( x, &e );

/* find significand 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 )
{
mtherr( fname, OVERFLOW );
return( MAXNUM );
}

if( w < MNEXP )
{
mtherr( fname, UNDERFLOW );
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 */

if( nflg )
{
/* For negative x,
* find out if the integer exponent
* is odd or even.
*/
w = ldexp( y, -1 );
w = floor(w);
w = ldexp( w, 1 );
if( w != y )
z = -z; /* odd exponent */
}

return( z );
}


/* Find a multiple of 1/16 that is within 1/16 of x. */
static double reduc(x)
double x;
{
double t;
double ldexp(), floor();

t = ldexp( x, 4 );
t = floor( t );
t = ldexp( t, -4 );
return(t);
}


  3 Responses to “Category : C Source Code
Archive   : CEPHES22.ZIP
Filename : POW.C

  1. Very nice! Thank you for this wonderful archive. I wonder why I found it only now. Long live the BBS file archives!

  2. This is so awesome! 😀 I’d be cool if you could download an entire archive of this at once, though.

  3. But one thing that puzzles me is the “mtswslnkmcjklsdlsbdmMICROSOFT” string. There is an article about it here. It is definitely worth a read: http://www.os2museum.com/wp/mtswslnk/