summaryrefslogtreecommitdiff
path: root/libm/float/expnf.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/expnf.c
parent22358dd7ce7bb49792204b698f01a6f69b9c8e08 (diff)
uClibc now has a math library. muahahahaha!
-Erik
Diffstat (limited to 'libm/float/expnf.c')
-rw-r--r--libm/float/expnf.c207
1 files changed, 207 insertions, 0 deletions
diff --git a/libm/float/expnf.c b/libm/float/expnf.c
new file mode 100644
index 000000000..ebf0ccb3e
--- /dev/null
+++ b/libm/float/expnf.c
@@ -0,0 +1,207 @@
+/* expnf.c
+ *
+ * Exponential integral En
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * int n;
+ * float x, y, expnf();
+ *
+ * y = expnf( n, x );
+ *
+ *
+ *
+ * DESCRIPTION:
+ *
+ * Evaluates the exponential integral
+ *
+ * inf.
+ * -
+ * | | -xt
+ * | e
+ * E (x) = | ---- dt.
+ * n | n
+ * | | t
+ * -
+ * 1
+ *
+ *
+ * Both n and x must be nonnegative.
+ *
+ * The routine employs either a power series, a continued
+ * fraction, or an asymptotic formula depending on the
+ * relative values of n and x.
+ *
+ * ACCURACY:
+ *
+ * Relative error:
+ * arithmetic domain # trials peak rms
+ * IEEE 0, 30 10000 5.6e-7 1.2e-7
+ *
+ */
+
+/* expn.c */
+
+/* Cephes Math Library Release 2.2: July, 1992
+ * Copyright 1985, 1992 by Stephen L. Moshier
+ * Direct inquiries to 30 Frost Street, Cambridge, MA 02140 */
+
+#include <math.h>
+
+#define EUL 0.57721566490153286060
+#define BIG 16777216.
+extern float MAXNUMF, MACHEPF, MAXLOGF;
+#ifdef ANSIC
+float powf(float, float), gammaf(float), logf(float), expf(float);
+#else
+float powf(), gammaf(), logf(), expf();
+#endif
+#define fabsf(x) ( (x) < 0 ? -(x) : (x) )
+
+
+float expnf( int n, float xx )
+{
+float x, ans, r, t, yk, xk;
+float pk, pkm1, pkm2, qk, qkm1, qkm2;
+float psi, z;
+int i, k;
+static float big = BIG;
+
+
+x = xx;
+if( n < 0 )
+ goto domerr;
+
+if( x < 0 )
+ {
+domerr: mtherr( "expnf", DOMAIN );
+ return( MAXNUMF );
+ }
+
+if( x > MAXLOGF )
+ return( 0.0 );
+
+if( x == 0.0 )
+ {
+ if( n < 2 )
+ {
+ mtherr( "expnf", SING );
+ return( MAXNUMF );
+ }
+ else
+ return( 1.0/(n-1.0) );
+ }
+
+if( n == 0 )
+ return( expf(-x)/x );
+
+/* expn.c */
+/* Expansion for large n */
+
+if( n > 5000 )
+ {
+ xk = x + n;
+ yk = 1.0 / (xk * xk);
+ t = n;
+ ans = yk * t * (6.0 * x * x - 8.0 * t * x + t * t);
+ ans = yk * (ans + t * (t - 2.0 * x));
+ ans = yk * (ans + t);
+ ans = (ans + 1.0) * expf( -x ) / xk;
+ goto done;
+ }
+
+if( x > 1.0 )
+ goto cfrac;
+
+/* expn.c */
+
+/* Power series expansion */
+
+psi = -EUL - logf(x);
+for( i=1; i<n; i++ )
+ psi = psi + 1.0/i;
+
+z = -x;
+xk = 0.0;
+yk = 1.0;
+pk = 1.0 - n;
+if( n == 1 )
+ ans = 0.0;
+else
+ ans = 1.0/pk;
+do
+ {
+ xk += 1.0;
+ yk *= z/xk;
+ pk += 1.0;
+ if( pk != 0.0 )
+ {
+ ans += yk/pk;
+ }
+ if( ans != 0.0 )
+ t = fabsf(yk/ans);
+ else
+ t = 1.0;
+ }
+while( t > MACHEPF );
+k = xk;
+t = n;
+r = n - 1;
+ans = (powf(z, r) * psi / gammaf(t)) - ans;
+goto done;
+
+/* expn.c */
+/* continued fraction */
+cfrac:
+k = 1;
+pkm2 = 1.0;
+qkm2 = x;
+pkm1 = 1.0;
+qkm1 = x + n;
+ans = pkm1/qkm1;
+
+do
+ {
+ k += 1;
+ if( k & 1 )
+ {
+ yk = 1.0;
+ xk = n + (k-1)/2;
+ }
+ else
+ {
+ yk = x;
+ xk = k/2;
+ }
+ pk = pkm1 * yk + pkm2 * xk;
+ qk = qkm1 * yk + qkm2 * xk;
+ if( qk != 0 )
+ {
+ r = pk/qk;
+ t = fabsf( (ans - r)/r );
+ ans = r;
+ }
+ else
+ t = 1.0;
+ pkm2 = pkm1;
+ pkm1 = pk;
+ qkm2 = qkm1;
+ qkm1 = qk;
+if( fabsf(pk) > big )
+ {
+ pkm2 *= MACHEPF;
+ pkm1 *= MACHEPF;
+ qkm2 *= MACHEPF;
+ qkm1 *= MACHEPF;
+ }
+ }
+while( t > MACHEPF );
+
+ans *= expf( -x );
+
+done:
+return( ans );
+}
+