/* powi.c * * Real raised to integer power * * * * SYNOPSIS: * * float x, y, __powif(); * int n; * * y = powi( 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 * DEC .04,26 -26,26 100000 2.7e-16 4.3e-17 * IEEE .04,26 -26,26 50000 2.0e-15 3.8e-16 * IEEE 1,2 -1022,1023 50000 8.6e-14 1.6e-14 * * Returns MAXNUM on overflow, zero on underflow. * */ /* powi.c */ /* Cephes Math Library Release 2.8: June, 2000 Copyright 1984, 1995, 2000 by Stephen L. Moshier */ /* Modified for float from powi.c and adapted to mingw 2002-10-01 Danny Smith */ #ifdef __MINGW32__ #include "cephes_mconf.h" #else #include "mconf.h" #ifdef ANSIPROT extern float logf ( float ); extern float frexpf ( float, int * ); extern int signbitf ( float ); #else float logf(), frexpf(); int signbitf(); #endif extern float NEGZEROF, INFINITYF, MAXNUMF, MAXLOGF, MINLOGF, LOGE2F; #endif /* __MINGW32__ */ #ifndef _SET_ERRNO #define _SET_ERRNO(x) #endif float __powif( float x, int nn ) { int n, e, sign, asign, lx; float w, y, s; /* See pow.c for these tests. */ if( x == 0.0F ) { if( nn == 0 ) return( 1.0F ); else if( nn < 0 ) return( INFINITYF ); else { if( nn & 1 ) return( x ); else return( 0.0 ); } } if( nn == 0 ) return( 1.0 ); if( nn == -1 ) return( 1.0/x ); if( x < 0.0 ) { asign = -1; x = -x; } else asign = 0; if( nn < 0 ) { sign = -1; n = -nn; } else { sign = 1; n = nn; } /* Even power will be positive. */ if( (n & 1) == 0 ) asign = 0; /* Overflow detection */ /* Calculate approximate logarithm of answer */ s = frexpf( x, &lx ); e = (lx - 1)*n; if( (e == 0) || (e > 64) || (e < -64) ) { s = (s - 7.0710678118654752e-1) / (s + 7.0710678118654752e-1); s = (2.9142135623730950 * s - 0.5 + lx) * nn * LOGE2F; } else { s = LOGE2F * e; } if( s > MAXLOGF ) { mtherr( "__powif", OVERFLOW ); _SET_ERRNO(ERANGE); y = INFINITYF; goto done; } #if DENORMAL if( s < MINLOGF ) { y = 0.0; goto done; } /* 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 < (-MAXLOGF+2.0)) && (sign < 0) ) { x = 1.0/x; sign = -sign; } #else /* do not produce denormal answer */ if( s < -MAXLOGF ) return(0.0); #endif /* First bit of the power */ if( n & 1 ) y = x; else y = 1.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; } if( sign < 0 ) y = 1.0/y; done: if( asign ) { /* odd power of negative number */ if( y == 0.0 ) y = NEGZEROF; else y = -y; } return(y); }