1 files changed, 208 insertions, 0 deletions
diff --git a/libm/double/expn.c b/libm/double/expn.c
new file mode 100644
index 000000000..89b6b139e
--- /dev/null
+++ b/libm/double/expn.c
@@ -0,0 +1,208 @@
+/*							expn.c
+ *
+ *		Exponential integral En
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * int n;
+ * double x, y, expn();
+ *
+ * y = expn( 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
+ *    DEC       0, 30        5000       2.0e-16     4.6e-17
+ *    IEEE      0, 30       10000       1.7e-15     3.6e-16
+ *
+ */
+
+/*							expn.c	*/
+
+/* Cephes Math Library Release 2.8:  June, 2000
+   Copyright 1985, 2000 by Stephen L. Moshier */
+
+#include <math.h>
+#ifdef ANSIPROT
+extern double pow ( double, double );
+extern double gamma ( double );
+extern double log ( double );
+extern double exp ( double );
+extern double fabs ( double );
+#else
+double pow(), gamma(), log(), exp(), fabs();
+#endif
+#define EUL 0.57721566490153286060
+#define BIG  1.44115188075855872E+17
+extern double MAXNUM, MACHEP, MAXLOG;
+
+double expn( n, x )
+int n;
+double x;
+{
+double ans, r, t, yk, xk;
+double pk, pkm1, pkm2, qk, qkm1, qkm2;
+double psi, z;
+int i, k;
+static double big = BIG;
+
+if( n < 0 )
+	goto domerr;
+
+if( x < 0 )
+	{
+domerr:	mtherr( "expn", DOMAIN );
+	return( MAXNUM );
+	}
+
+if( x > MAXLOG )
+	return( 0.0 );
+
+if( x == 0.0 )
+	{
+	if( n < 2 )
+		{
+		mtherr( "expn", SING );
+		return( MAXNUM );
+		}
+	else
+		return( 1.0/(n-1.0) );
+	}
+
+if( n == 0 )
+	return( exp(-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) * exp( -x ) / xk;
+	goto done;
+	}
+
+if( x > 1.0 )
+	goto cfrac;
+
+/*							expn.c	*/
+
+/*		Power series expansion		*/
+
+psi = -EUL - log(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 = fabs(yk/ans);
+	else
+		t = 1.0;
+	}
+while( t > MACHEP );
+k = xk;
+t = n;
+r = n - 1;
+ans = (pow(z, r) * psi / gamma(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 = fabs( (ans - r)/r );
+		ans = r;
+		}
+	else
+		t = 1.0;
+	pkm2 = pkm1;
+	pkm1 = pk;
+	qkm2 = qkm1;
+	qkm1 = qk;
+if( fabs(pk) > big )
+		{
+		pkm2 /= big;
+		pkm1 /= big;
+		qkm2 /= big;
+		qkm1 /= big;
+		}
+	}
+while( t > MACHEP );
+
+ans *= exp( -x );
+
+done:
+return( ans );
+}
+