diff options
author | Eric Andersen <andersen@codepoet.org> | 2001-11-22 14:04:29 +0000 |
---|---|---|
committer | Eric Andersen <andersen@codepoet.org> | 2001-11-22 14:04:29 +0000 |
commit | 7ce331c01ce6eb7b3f5c715a38a24359da9c6ee2 (patch) | |
tree | 3a7e8476e868ae15f4da1b7ce26b2db6f434468c /libm/float/powif.c | |
parent | c117dd5fb183afb1a4790a6f6110d88704be6bf8 (diff) |
Totally rework the math library, this time based on the MacOs X
math library (which is itself based on the math lib from FreeBSD).
-Erik
Diffstat (limited to 'libm/float/powif.c')
-rw-r--r-- | libm/float/powif.c | 156 |
1 files changed, 0 insertions, 156 deletions
diff --git a/libm/float/powif.c b/libm/float/powif.c deleted file mode 100644 index d226896ba..000000000 --- a/libm/float/powif.c +++ /dev/null @@ -1,156 +0,0 @@ -/* powif.c - * - * Real raised to integer power - * - * - * - * SYNOPSIS: - * - * float x, y, powif(); - * int n; - * - * y = powif( 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 .04,26 -26,26 100000 1.1e-6 2.0e-7 - * IEEE 1,2 -128,128 100000 1.1e-5 1.0e-6 - * - * Returns MAXNUMF on overflow, zero on underflow. - * - */ - -/* powi.c */ - -/* -Cephes Math Library Release 2.2: June, 1992 -Copyright 1984, 1987, 1989 by Stephen L. Moshier -Direct inquiries to 30 Frost Street, Cambridge, MA 02140 -*/ - -#include <math.h> -extern float MAXNUMF, MAXLOGF, MINLOGF, LOGE2F; - -float frexpf( float, int * ); - -float powif( float x, int nn ) -{ -int n, e, sign, asign, lx; -float w, y, s; - -if( x == 0.0 ) - { - if( nn == 0 ) - return( 1.0 ); - else if( nn < 0 ) - return( MAXNUMF ); - else - return( 0.0 ); - } - -if( nn == 0 ) - return( 1.0 ); - - -if( x < 0.0 ) - { - 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 = 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( "powi", OVERFLOW ); - y = MAXNUMF; - goto done; - } - -if( s < MINLOGF ) - return(0.0); - -/* 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) ) - { - x = 1.0/x; - sign = 0; - } - -/* First bit of the power */ -if( n & 1 ) - y = x; - -else - { - y = 1.0; - 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.0/y; -return(y); -} |