diff options
Diffstat (limited to 'libm/double/expn.c')
-rw-r--r-- | libm/double/expn.c | 208 |
1 files changed, 0 insertions, 208 deletions
diff --git a/libm/double/expn.c b/libm/double/expn.c deleted file mode 100644 index 89b6b139e..000000000 --- a/libm/double/expn.c +++ /dev/null @@ -1,208 +0,0 @@ -/* 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 ); -} - |