view release on metacpan or  search on metacpan

libmd/expn.c  view on Meta::CPAN

/*							md_expn.c
 *
 *		Exponential integral En
 *
 *
 *
 * SYNOPSIS:
 *
 * int n;
 * double x, y, md_expn();
 *
 * y = md_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
 *
 */

/*							md_expn.c	*/

/* Cephes Math Library Release 2.8:  June, 2000
   Copyright 1985, 2000 by Stephen L. Moshier */

#include "mconf.h"
#ifdef ANSIPROT
extern double md_pow ( double, double );
extern double md_gamma ( double );
extern double md_log ( double );
extern double md_exp ( double );
extern double md_fabs ( double );
#else
double md_pow(), md_gamma(), md_log(), md_exp(), md_fabs();
#endif
#define EUL 0.57721566490153286060
#define BIG  1.44115188075855872E+17
extern double MAXNUM, MACHEP, MAXLOG;

double md_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( "md_expn", DOMAIN );
	return( MAXNUM );
	}

if( x > MAXLOG )
	return( 0.0 );

if( x == 0.0 )
	{
	if( n < 2 )
		{
		mtherr( "md_expn", SING );
		return( MAXNUM );
		}
	else