powl.c

linux下用PCMCIA无线网卡虚拟无线AP的程序源码

extern long double polevll ( long double, void *, int );
extern long double p1evll ( long double, void *, int );
extern long double powil ( long double, int );
extern int isnanl ( long double );
extern int isfinitel ( long double );
static long double reducl( long double );
extern int signbitl ( long double );
#else
long double floorl(), fabsl(), frexpl(), ldexpl();
long double polevll(), p1evll(), powil();
static long double reducl();
int isnanl(), isfinitel(), signbitl();
#endif

#ifdef INFINITIES
extern long double INFINITYL;
#else
#define INFINITYL MAXNUML
#endif

#ifdef NANS
extern long double NANL;
#endif
#ifdef MINUSZERO
extern long double NEGZEROL;
#endif

long double powl( x, y )
long double x, y;
{
/* double F, Fa, Fb, G, Ga, Gb, H, Ha, Hb */
int i, nflg, iyflg, yoddint;
long e;

if( y == 0.0L )
	return( 1.0L );

#ifdef NANS
if( isnanl(x) )
	return( x );
if( isnanl(y) )
	return( y );
#endif

if( y == 1.0L )
	return( x );

#ifdef INFINITIES
if( !isfinitel(y) && (x == -1.0L || x == 1.0L) )
	{
	mtherr( "powl", DOMAIN );
#ifdef NANS
	return( NANL );
#else
	return( INFINITYL );
#endif
	}
#endif

if( x == 1.0L )
	return( 1.0L );

if( y >= MAXNUML )
	{
#ifdef INFINITIES
	if( x > 1.0L )
		return( INFINITYL );
#else
	if( x > 1.0L )
		return( MAXNUML );
#endif
	if( x > 0.0L && x < 1.0L )
		return( 0.0L );
#ifdef INFINITIES
	if( x < -1.0L )
		return( INFINITYL );
#else
	if( x < -1.0L )
		return( MAXNUML );
#endif
	if( x > -1.0L && x < 0.0L )
		return( 0.0L );
	}
if( y <= -MAXNUML )
	{
	if( x > 1.0L )
		return( 0.0L );
#ifdef INFINITIES
	if( x > 0.0L && x < 1.0L )
		return( INFINITYL );
#else
	if( x > 0.0L && x < 1.0L )
		return( MAXNUML );
#endif
	if( x < -1.0L )
		return( 0.0L );
#ifdef INFINITIES
	if( x > -1.0L && x < 0.0L )
		return( INFINITYL );
#else
	if( x > -1.0L && x < 0.0L )
		return( MAXNUML );
#endif
	}
if( x >= MAXNUML )
	{
#if INFINITIES
	if( y > 0.0L )
		return( INFINITYL );
#else
	if( y > 0.0L )
		return( MAXNUML );
#endif
	return( 0.0L );
	}

w = floorl(y);
/* Set iyflg to 1 if y is an integer.  */
iyflg = 0;
if( w == y )
	iyflg = 1;

/* Test for odd integer y.  */
yoddint = 0;
if( iyflg )
	{
	ya = fabsl(y);
	ya = floorl(0.5L * ya);
	yb = 0.5L * fabsl(w);
	if( ya != yb )
		yoddint = 1;
	}

if( x <= -MAXNUML )
	{
	if( y > 0.0L )
		{
#ifdef INFINITIES
		if( yoddint )
			return( -INFINITYL );
		return( INFINITYL );
#else
		if( yoddint )
			return( -MAXNUML );
		return( MAXNUML );
#endif
		}
	if( y < 0.0L )
		{
#ifdef MINUSZERO
		if( yoddint )
			return( NEGZEROL );
#endif
		return( 0.0 );
		}
 	}


nflg = 0;	/* flag = 1 if x<0 raised to integer power */
if( x <= 0.0L )
	{
	if( x == 0.0L )
		{
		if( y < 0.0 )
			{
#ifdef MINUSZERO
			if( signbitl(x) && yoddint )
				return( -INFINITYL );
#endif
#ifdef INFINITIES
			return( INFINITYL );
#else
			return( MAXNUML );
#endif
			}
		if( y > 0.0 )
			{
#ifdef MINUSZERO
			if( signbitl(x) && yoddint )
				return( NEGZEROL );
#endif
			return( 0.0 );
			}
		if( y == 0.0L )
			return( 1.0L );  /*   0**0   */
		else  
			return( 0.0L );  /*   0**y   */
		}
	else
		{
		if( iyflg == 0 )
			{ /* noninteger power of negative number */
			mtherr( fname, DOMAIN );
#ifdef NANS
			return(NANL);
#else
			return(0.0L);
#endif
			}
		nflg = 1;
		}
	}

/* Integer power of an integer.  */

if( iyflg )
	{
	i = w;
	w = floorl(x);
	if( (w == x) && (fabsl(y) < 32768.0) )
		{
		w = powil( x, (int) y );
		return( w );
		}
	}


if( nflg )
	x = fabsl(x);

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

/* find significand in antilog table A[] */
i = 1;
if( x <= douba(17) )
	i = 17;
if( x <= douba(i+8) )
	i += 8;
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 * polevll( x, P, 3 ) / p1evll( x, Q, 3 ) );
w = w - ldexpl( z, -1 );   /*  w - 0.5 * z  */

/* Convert to base 2 logarithm:
 * multiply by log2(e) = 1 + LOG2EA
 */
z = LOG2EA * w;
z += w;
z += LOG2EA * x;
z += x;

/* Compute exponent term of the base 2 logarithm. */
w = -i;
w = ldexpl( w, -LNXT );	/* divide by NXT */
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/NXT
 */
ya = reducl(y);
yb = y - ya;

/* (w+z)(ya+yb)
 * = w*ya + w*yb + z*y
 */
F = z * y  +  w * yb;
Fa = reducl(F);
Fb = F - Fa;

G = Fa + w * ya;
Ga = reducl(G);
Gb = G - Ga;

H = Fb + Gb;
Ha = reducl(H);
w = ldexpl( Ga+Ha, LNXT );

/* Test the power of 2 for overflow */
if( w > MEXP )
	{
/*	printf( "w = %.4Le ", w ); */
	mtherr( fname, OVERFLOW );
	return( MAXNUML );
	}

if( w < MNEXP )
	{
/*	printf( "w = %.4Le ", w ); */
	mtherr( fname, UNDERFLOW );
	return( 0.0L );
	}

e = w;
Hb = H - Ha;

if( Hb > 0.0L )
	{
	e += 1;
	Hb -= (1.0L/NXT);  /*0.0625L;*/
	}

/* Now the product y * log2(x)  =  Hb + e/NXT.
 *
 * Compute base 2 exponential of Hb,
 * where -0.0625 <= Hb <= 0.
 */
z = Hb * polevll( Hb, R, 6 );  /*    z  =  2**Hb - 1    */

/* Express e/NXT as an integer plus a negative number of (1/NXT)ths.
 * Find lookup table entry for the fractional power of 2.
 */
if( e < 0 )
	i = 0;
else
	i = 1;
i = e/NXT + i;
e = NXT*i - e;
w = douba( e );
z = w * z;      /*    2**-e * ( 1 + (2**Hb-1) )    */
z = z + w;
z = ldexpl( z, i );  /* multiply by integer power of 2 */

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

return( z );
}


/* Find a multiple of 1/NXT that is within 1/NXT of x. */
static long double reducl(x)
long double x;
{
long double t;

t = ldexpl( x, LNXT );
t = floorl( t );
t = ldexpl( t, -LNXT );
return(t);
}