summaryrefslogtreecommitdiff
path: root/libm/double/powi.c
diff options
context:
space:
mode:
authorEric Andersen <andersen@codepoet.org>2001-05-10 00:40:28 +0000
committerEric Andersen <andersen@codepoet.org>2001-05-10 00:40:28 +0000
commit1077fa4d772832f77a677ce7fb7c2d513b959e3f (patch)
tree579bee13fb0b58d2800206366ec2caecbb15f3fc /libm/double/powi.c
parent22358dd7ce7bb49792204b698f01a6f69b9c8e08 (diff)
uClibc now has a math library. muahahahaha!
-Erik
Diffstat (limited to 'libm/double/powi.c')
-rw-r--r--libm/double/powi.c186
1 files changed, 186 insertions, 0 deletions
diff --git a/libm/double/powi.c b/libm/double/powi.c
new file mode 100644
index 000000000..46d9a1400
--- /dev/null
+++ b/libm/double/powi.c
@@ -0,0 +1,186 @@
+/* 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);
+}