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 );
+}
+