3 * Real raised to integer power
18 * Returns argument x raised to the nth power.
19 * The routine efficiently decomposes n as a sum of powers of
20 * two. The desired power is a product of two-to-the-kth
21 * powers of x. Thus to compute the 32767 power of x requires
22 * 28 multiplications instead of 32767 multiplications.
30 * arithmetic x domain n domain # trials peak rms
31 * IEEE .04,26 -26,26 100000 1.1e-6 2.0e-7
32 * IEEE 1,2 -128,128 100000 1.1e-5 1.0e-6
34 * Returns MAXNUMF on overflow, zero on underflow.
41 Cephes Math Library Release 2.2: June, 1992
42 Copyright 1984, 1987, 1989 by Stephen L. Moshier
43 Direct inquiries to 30 Frost Street, Cambridge, MA 02140
47 extern float MAXNUMF, MAXLOGF, MINLOGF, LOGE2F;
49 float frexpf( float, int * );
51 float powif( float x, int nn )
53 int n, e, sign, asign, lx;
93 /* Overflow detection */
95 /* Calculate approximate logarithm of answer */
98 if( (e == 0) || (e > 64) || (e < -64) )
100 s = (s - 7.0710678118654752e-1) / (s + 7.0710678118654752e-1);
101 s = (2.9142135623730950 * s - 0.5 + lx) * nn * LOGE2F;
110 mtherr( "powi", OVERFLOW );
118 /* Handle tiny denormal answer, but with less accuracy
119 * since roundoff error in 1.0/x will be amplified.
120 * The precise demarcation should be the gradual underflow threshold.
122 if( s < (-MAXLOGF+2.0) )
128 /* First bit of the power */
142 w = w * w; /* arg to the 2-to-the-kth power */
143 if( n & 1 ) /* if that bit is set, then include in product */
152 y = -y; /* odd power of negative number */