summaryrefslogtreecommitdiff
path: root/libm/float/powf.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/float/powf.c
parent22358dd7ce7bb49792204b698f01a6f69b9c8e08 (diff)
uClibc now has a math library. muahahahaha!
-Erik
Diffstat (limited to 'libm/float/powf.c')
-rw-r--r--libm/float/powf.c338
1 files changed, 338 insertions, 0 deletions
diff --git a/libm/float/powf.c b/libm/float/powf.c
new file mode 100644
index 000000000..367a39ad4
--- /dev/null
+++ b/libm/float/powf.c
@@ -0,0 +1,338 @@
+/* powf.c
+ *
+ * Power function
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * float x, y, z, powf();
+ *
+ * z = powf( x, y );
+ *
+ *
+ *
+ * DESCRIPTION:
+ *
+ * Computes x raised to the yth power. Analytically,
+ *
+ * x**y = exp( y log(x) ).
+ *
+ * Following Cody and Waite, this program uses a lookup table
+ * of 2**-i/16 and pseudo extended precision arithmetic to
+ * obtain an extra three bits of accuracy in both the logarithm
+ * and the exponential.
+ *
+ *
+ *
+ * ACCURACY:
+ *
+ * Relative error:
+ * arithmetic domain # trials peak rms
+ * IEEE -10,10 100,000 1.4e-7 3.6e-8
+ * 1/10 < x < 10, x uniformly distributed.
+ * -10 < y < 10, y uniformly distributed.
+ *
+ *
+ * ERROR MESSAGES:
+ *
+ * message condition value returned
+ * powf overflow x**y > MAXNUMF MAXNUMF
+ * powf underflow x**y < 1/MAXNUMF 0.0
+ * powf domain x<0 and y noninteger 0.0
+ *
+ */
+
+/*
+Cephes Math Library Release 2.2: June, 1992
+Copyright 1984, 1987, 1988 by Stephen L. Moshier
+Direct inquiries to 30 Frost Street, Cambridge, MA 02140
+*/
+
+
+#include <math.h>
+static char fname[] = {"powf"};
+
+
+/* 2^(-i/16)
+ * The decimal values are rounded to 24-bit precision
+ */
+static float A[] = {
+ 1.00000000000000000000E0,
+ 9.57603275775909423828125E-1,
+ 9.17004048824310302734375E-1,
+ 8.78126084804534912109375E-1,
+ 8.40896427631378173828125E-1,
+ 8.05245161056518554687500E-1,
+ 7.71105408668518066406250E-1,
+ 7.38413095474243164062500E-1,
+ 7.07106769084930419921875E-1,
+ 6.77127778530120849609375E-1,
+ 6.48419797420501708984375E-1,
+ 6.20928883552551269531250E-1,
+ 5.94603538513183593750000E-1,
+ 5.69394290447235107421875E-1,
+ 5.45253872871398925781250E-1,
+ 5.22136867046356201171875E-1,
+ 5.00000000000000000000E-1
+};
+/* continuation, for even i only
+ * 2^(i/16) = A[i] + B[i/2]
+ */
+static float B[] = {
+ 0.00000000000000000000E0,
+-5.61963907099083340520586E-9,
+-1.23776636307969995237668E-8,
+ 4.03545234539989593104537E-9,
+ 1.21016171044789693621048E-8,
+-2.00949968760174979411038E-8,
+ 1.89881769396087499852802E-8,
+-6.53877009617774467211965E-9,
+ 0.00000000000000000000E0
+};
+
+/* 1 / A[i]
+ * The decimal values are full precision
+ */
+static float Ainv[] = {
+ 1.00000000000000000000000E0,
+ 1.04427378242741384032197E0,
+ 1.09050773266525765920701E0,
+ 1.13878863475669165370383E0,
+ 1.18920711500272106671750E0,
+ 1.24185781207348404859368E0,
+ 1.29683955465100966593375E0,
+ 1.35425554693689272829801E0,
+ 1.41421356237309504880169E0,
+ 1.47682614593949931138691E0,
+ 1.54221082540794082361229E0,
+ 1.61049033194925430817952E0,
+ 1.68179283050742908606225E0,
+ 1.75625216037329948311216E0,
+ 1.83400808640934246348708E0,
+ 1.91520656139714729387261E0,
+ 2.00000000000000000000000E0
+};
+
+#ifdef DEC
+#define MEXP 2032.0
+#define MNEXP -2032.0
+#else
+#define MEXP 2048.0
+#define MNEXP -2400.0
+#endif
+
+/* log2(e) - 1 */
+#define LOG2EA 0.44269504088896340736F
+extern float MAXNUMF;
+
+#define F W
+#define Fa Wa
+#define Fb Wb
+#define G W
+#define Ga Wa
+#define Gb u
+#define H W
+#define Ha Wb
+#define Hb Wb
+
+
+#ifdef ANSIC
+float floorf( float );
+float frexpf( float, int *);
+float ldexpf( float, int );
+float powif( float, int );
+#else
+float floorf(), frexpf(), ldexpf(), powif();
+#endif
+
+/* Find a multiple of 1/16 that is within 1/16 of x. */
+#define reduc(x) 0.0625 * floorf( 16 * (x) )
+
+#ifdef ANSIC
+float powf( float x, float y )
+#else
+float powf( x, y )
+float x, y;
+#endif
+{
+float u, w, z, W, Wa, Wb, ya, yb;
+/* float F, Fa, Fb, G, Ga, Gb, H, Ha, Hb */
+int e, i, nflg;
+
+
+nflg = 0; /* flag = 1 if x<0 raised to integer power */
+w = floorf(y);
+if( w < 0 )
+ z = -w;
+else
+ z = w;
+if( (w == y) && (z < 32768.0) )
+ {
+ i = w;
+ w = powif( x, i );
+ return( w );
+ }
+
+
+if( x <= 0.0F )
+ {
+ if( x == 0.0 )
+ {
+ if( y == 0.0 )
+ return( 1.0 ); /* 0**0 */
+ else
+ return( 0.0 ); /* 0**y */
+ }
+ else
+ {
+ if( w != y )
+ { /* noninteger power of negative number */
+ mtherr( fname, DOMAIN );
+ return(0.0);
+ }
+ nflg = 1;
+ if( x < 0 )
+ x = -x;
+ }
+ }
+
+/* separate significand from exponent */
+x = frexpf( x, &e );
+
+/* find significand in antilog table A[] */
+i = 1;
+if( x <= A[9] )
+ i = 9;
+if( x <= A[i+4] )
+ i += 4;
+if( x <= A[i+2] )
+ i += 2;
+if( x >= A[1] )
+ i = -1;
+i += 1;
+
+
+/* Find (x - A[i])/A[i]
+ * in order to compute log(x/A[i]):
+ *
+ * log(x) = log( a x/a ) = log(a) + log(x/a)
+ *
+ * log(x/a) = log(1+v), v = x/a - 1 = (x-a)/a
+ */
+x -= A[i];
+x -= B[ i >> 1 ];
+x *= Ainv[i];
+
+
+/* rational approximation for log(1+v):
+ *
+ * log(1+v) = v - 0.5 v^2 + v^3 P(v)
+ * Theoretical relative error of the approximation is 3.5e-11
+ * on the interval 2^(1/16) - 1 > v > 2^(-1/16) - 1
+ */
+z = x*x;
+w = (((-0.1663883081054895 * x
+ + 0.2003770364206271) * x
+ - 0.2500006373383951) * x
+ + 0.3333331095506474) * x * z;
+w -= 0.5 * z;
+
+/* Convert to base 2 logarithm:
+ * multiply by log2(e)
+ */
+w = w + LOG2EA * w;
+/* Note x was not yet added in
+ * to above rational approximation,
+ * so do it now, while multiplying
+ * by log2(e).
+ */
+z = w + LOG2EA * x;
+z = z + x;
+
+/* Compute exponent term of the base 2 logarithm. */
+w = -i;
+w *= 0.0625; /* divide by 16 */
+w += e;
+/* Now base 2 log of x is w + z. */
+
+/* Multiply base 2 log by y, in extended precision. */
+
+/* separate y into large part ya
+ * and small part yb less than 1/16
+ */
+ya = reduc(y);
+yb = y - ya;
+
+
+F = z * y + w * yb;
+Fa = reduc(F);
+Fb = F - Fa;
+
+G = Fa + w * ya;
+Ga = reduc(G);
+Gb = G - Ga;
+
+H = Fb + Gb;
+Ha = reduc(H);
+w = 16 * (Ga + Ha);
+
+/* Test the power of 2 for overflow */
+if( w > MEXP )
+ {
+ mtherr( fname, OVERFLOW );
+ return( MAXNUMF );
+ }
+
+if( w < MNEXP )
+ {
+ mtherr( fname, UNDERFLOW );
+ return( 0.0 );
+ }
+
+e = w;
+Hb = H - Ha;
+
+if( Hb > 0.0 )
+ {
+ e += 1;
+ Hb -= 0.0625;
+ }
+
+/* Now the product y * log2(x) = Hb + e/16.0.
+ *
+ * Compute base 2 exponential of Hb,
+ * where -0.0625 <= Hb <= 0.
+ * Theoretical relative error of the approximation is 2.8e-12.
+ */
+/* z = 2**Hb - 1 */
+z = ((( 9.416993633606397E-003 * Hb
+ + 5.549356188719141E-002) * Hb
+ + 2.402262883964191E-001) * Hb
+ + 6.931471791490764E-001) * Hb;
+
+/* Express e/16 as an integer plus a negative number of 16ths.
+ * Find lookup table entry for the fractional power of 2.
+ */
+if( e < 0 )
+ i = -( -e >> 4 );
+else
+ i = (e >> 4) + 1;
+e = (i << 4) - e;
+w = A[e];
+z = w + w * z; /* 2**-e * ( 1 + (2**Hb-1) ) */
+z = ldexpf( z, i ); /* multiply by integer power of 2 */
+
+if( nflg )
+ {
+/* For negative x,
+ * find out if the integer exponent
+ * is odd or even.
+ */
+ w = 2 * floorf( (float) 0.5 * w );
+ if( w != y )
+ z = -z; /* odd exponent */
+ }
+
+return( z );
+}