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