123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207 |
- /* 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 );
- }
|