summaryrefslogtreecommitdiff
path: root/libm/double/expn.c
diff options
context:
space:
mode:
Diffstat (limited to 'libm/double/expn.c')
-rw-r--r--libm/double/expn.c208
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 );
+}
+