/*							powil.c * *	Real raised to integer power, long double precision * * * * SYNOPSIS: * * long double x, y, powil(); * int n; * * y = powil( x, n ); * * * * DESCRIPTION: * * Returns argument x raised to the nth power. * The routine efficiently decomposes n as a sum of powers of * two. The desired power is a product of two-to-the-kth * powers of x.  Thus to compute the 32767 power of x requires * 28 multiplications instead of 32767 multiplications. * * * * ACCURACY: * * *                      Relative error: * arithmetic   x domain   n domain  # trials      peak         rms *    IEEE     .001,1000  -1022,1023  100,000     7.5e-32     1.4e-32 *    IEEE     .99,1.01     0,8700    100,000     4.6e-31     9.1e-32 * * Returns MAXNUM on overflow, zero on underflow. * *//*							powil.c	*//*Cephes Math Library Release 2.2:  December, 1990Copyright 1984, 1990 by Stephen L. MoshierDirect inquiries to 30 Frost Street, Cambridge, MA 02140*/#include "mconf.h"extern long double MAXNUML, MAXLOGL, MINLOGL;extern long double LOGE2L;long double powil( x, nn )long double x;int nn;{long double w, y;long double s;int n, e, sign, asign, lx;long double frexpl();if( x == 0.0L )	{	if( nn == 0 )		return( 1.0L );	else if( nn < 0 )		return( MAXNUML );	else		return( 0.0L );	}if( nn == 0 )	return( 1.0L );if( x < 0.0L )	{	asign = -1;	x = -x;	}else	asign = 0;if( nn < 0 )	{	sign = -1;	n = -nn;/*	x = 1.0/x;*/	}else	{	sign = 0;	n = nn;	}/* Overflow detection *//* Calculate approximate logarithm of answer */s = x;s = frexpl( s, &lx );e = (lx - 1)*n;if( (e == 0) || (e > 64) || (e < -64) )	{	s = (s - 7.0710678118654752e-1L) / (s +  7.0710678118654752e-1L);	s = (2.9142135623730950L * s - 0.5L + lx) * nn * LOGE2L;	}else	{	s = LOGE2L * e;	}if( s > MAXLOGL )	{	mtherr( "powil", OVERFLOW );	y = MAXNUML;	goto done;	}if( s < MINLOGL )	return(0.0L);/* Handle tiny denormal answer, but with less accuracy * since roundoff error in 1.0/x will be amplified. * The precise demarcation should be the gradual underflow threshold. */if( s < (-MAXLOGL+2.0L) )	{	x = 1.0L/x;	sign = 0;	}/* First bit of the power */if( n & 1 )	y = x;		else	{	y = 1.0L;	asign = 0;	}w = x;n >>= 1;while( n )	{	w = w * w;	/* arg to the 2-to-the-kth power */	if( n & 1 )	/* if that bit is set, then include in product */		y *= w;	n >>= 1;	}done:if( asign )	y = -y; /* odd power of negative number */if( sign )	y = 1.0L/y;return(y);}