diff options
Diffstat (limited to 'libm/double/powi.c')
-rw-r--r-- | libm/double/powi.c | 186 |
1 files changed, 0 insertions, 186 deletions
diff --git a/libm/double/powi.c b/libm/double/powi.c deleted file mode 100644 index 46d9a1400..000000000 --- a/libm/double/powi.c +++ /dev/null @@ -1,186 +0,0 @@ -/* powi.c - * - * Real raised to integer power - * - * - * - * SYNOPSIS: - * - * double x, y, powi(); - * 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 -*/ - -#include <math.h> -#ifdef ANSIPROT -extern double log ( double ); -extern double frexp ( double, int * ); -extern int signbit ( double ); -#else -double log(), frexp(); -int signbit(); -#endif -extern double NEGZERO, INFINITY, MAXNUM, MAXLOG, MINLOG, LOGE2; - -double powi( x, nn ) -double x; -int nn; -{ -int n, e, sign, asign, lx; -double w, y, s; - -/* See pow.c for these tests. */ -if( x == 0.0 ) - { - if( nn == 0 ) - return( 1.0 ); - else if( nn < 0 ) - return( INFINITY ); - 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 = frexp( 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 * LOGE2; - } -else - { - s = LOGE2 * e; - } - -if( s > MAXLOG ) - { - mtherr( "powi", OVERFLOW ); - y = INFINITY; - goto done; - } - -#if DENORMAL -if( s < MINLOG ) - { - 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 < (-MAXLOG+2.0)) && (sign < 0) ) - { - x = 1.0/x; - sign = -sign; - } -#else -/* do not produce denormal answer */ -if( s < -MAXLOG ) - 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 = NEGZERO; - else - y = -y; - } -return(y); -} |