summaryrefslogtreecommitdiff
path: root/libm/double/powi.c
diff options
context:
space:
mode:
authorEric Andersen <andersen@codepoet.org>2001-11-22 14:04:29 +0000
committerEric Andersen <andersen@codepoet.org>2001-11-22 14:04:29 +0000
commit7ce331c01ce6eb7b3f5c715a38a24359da9c6ee2 (patch)
tree3a7e8476e868ae15f4da1b7ce26b2db6f434468c /libm/double/powi.c
parentc117dd5fb183afb1a4790a6f6110d88704be6bf8 (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/double/powi.c')
-rw-r--r--libm/double/powi.c186
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);
-}