/*							qexpn.c
 *
 *		Exponential integral En
 *
 *
 *
 * SYNOPSIS:
 *
 * int qexpn( n, x, y );
 * int n;
 * QELT *x, *y;
 *
 * qexpn( n, x, y );
 *
 *
 *
 * DESCRIPTION:
 *
 * Evaluates the exponential integral
 *
 *                 inf.
 *                   -
 *                  | |   -xt
 *                  |    e
 *      E (x)  =    |    ----  dt.
 *       n          |      n
 *                | |     t
 *                 -
 *                  1
 *
 *
 * Both n and x must be nonnegative.
 *
 *
 * ACCURACY:
 * 
 * Series expansions are truncated at less than full working precision.
 *
 */

/*							expn.c	*/

/* Cephes Math Library Release 1.1:  March, 1985
 * Copyright 1985 by Stephen L. Moshier */

/* Exponential integral	*/
#include "qhead.h"
extern QELT qone[];
extern QELT qeul[];

extern double MAXNUM;
extern int qgamma();

int qexpn( n, x, yy )
int n;
QELT x[], yy[];
{
static QELT ans[NQ], r[NQ], t[NQ], yk[NQ], xk[NQ], qn[NQ];
static QELT pk[NQ], pkm1[NQ], pkm2[NQ], qk[NQ], qkm1[NQ], qkm2[NQ];
static QELT psi[NQ], z[NQ];
int i, k;
long ln;
double temp;

if( n < 0 )
	{
	mtherr("qexpn", DOMAIN );
	goto overf;
	}

if( x[0] != 0 )
	{
	mtherr("qexpn", DOMAIN );
	goto overf;
	}

if( x[1] < 3 )
	{
	if( n < 2 )
		goto overf;
	else
		{
		ln = n - 1;
		ltoq( &ln, ans );
		qdiv( ans, qone, yy );
		return 0;
		}
	}

if( n == 0 )
	{
	qexp( x, ans );		/* exp(-x)/x */
	qmul( ans, x, ans );
	qdiv( ans, qone, yy );
	return 0;
	}

/*							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] > (QELT) (EXPONE+1) )
	goto cfrac;

/*							expn.c	*/

/*		Power series expansion		*/


qlog( x, psi );		/* psi = -EUL - log(x) */
qneg( psi );
qsub( qeul, psi, psi );

for( i=1; i<n; i++ )
	{
	ln = i;
	ltoq( &ln, qn );
	qdiv( qn, qone, qn );	/* psi = psi + 1.0/i */
	qadd( qn, psi, psi );
	}

qmov( x, z );		/* z = -x */
qneg( z );
qclear( xk );		/* xk = 0.0 */
qmov( qone, yk );	/* yk = 1.0 */
ln = n;
ltoq( &ln, qn );
qsub( qn, qone, pk );	/* pk = 1.0 - n */
if( n == 1 )
	qclear( ans );	/* ans = 0.0 */
else
	qdiv( pk, qone, ans );	/* ans = 1.0/pk */

do
	{
	qadd( qone, xk, xk );	/* xk += 1.0 */
	qdiv( xk, z, qn );	/* yk *= z/xk */
	qmul( qn, yk, yk );
	qadd( qone, pk, pk );	/* pk += 1.0 */
	if( pk[1] > 10 )
		{
		qdiv( pk, yk, t );	/* ans += yk/pk */
		qadd( t, ans, ans );
		}
	else
		qmov( qone, t );	/* t = 1.0 */
	}
while( ((int) t[1] - (int) ans[1]) > -70 );

qtoe( xk, (unsigned short *) &temp );
k = temp;	/* k = xk */
ln = n;
ltoq( &ln, t );	/* t = n */
/* ans = (powi(z, n-1) * psi / gamma(t)) - ans */
qgamma( t, t );
qlog( x, qn );
ln = n - 1;
ltoq( &ln, yk );
qmul( yk, qn, qn );
qexp( qn, qn );
if( ((n-1) & 1) != 0 )
	qneg( qn );
qmul( psi, qn, qn );
qdiv( t, qn, qn );
qsub( ans, qn, ans );
goto done;


/*							expn.c	*/
/*		continued fraction		*/
cfrac:
k = 1;
qmov( qone, pkm2 );	/* pkm2 = 1.0 */
qmov( x, qkm2 );	/* qkm2 = x   */
qmov( qone, pkm1 );	/* pkm1 = 1.0 */
ln = n;
ltoq( &ln, qn );
qadd( qn, x, qkm1 );	/* qkm1 = x + n */
qdiv( qkm1, pkm1, ans );	/* ans = pkm1/qkm1 */

do
	{
	k += 1;
	if( k & 1 )
		{
		qmov( qone, yk );	/* yk = 1.0 */

		ln = n + (k-1)/2;	/* xk = n + (k-1)/2 */
		ltoq( &ln, xk );
		}
	else
		{
		qmov( x, yk );		/* yk = x */
		ln = k/2;		/* xk = k/2 */
		ltoq( &ln, xk );
		}
	qmul( yk, pkm1, qn );	/* pk = pkm1 * yk  +  pkm2 * xk */
	qmul( xk, pkm2, pk );
	qadd( qn, pk, pk );	
	qmul( yk, qkm1, qn );	/* qk = qkm1 * yk  +  qkm2 * xk */
	qmul( xk, qkm2, qk );
	qadd( qn, qk, qk );
	if( qk[1] > 2 )
		{
		qdiv( qk, pk, r );	/* r = pk/qk */
		qsub( r, ans, t );	/* t = abs( (ans - r)/r ) */
		qmov( r, ans );		/* ans = r */
		}
	else
		qmov( qone, t );	/* t = 1.0 */
	qmov( pkm1, pkm2 );		/* pkm2 = pkm1 */
	qmov( pk, pkm1 );		/* pkm1 = pk   */
	qmov( qkm1, qkm2 );		/* qkm2 = qkm1 */
	qmov( qk, qkm1 );		/* qkm1 = qk   */
if( pk[1] > (QELT) (EXPONE + 64) )
	{
	pkm2[1] -= 64;
	pkm1[1] -= 64;
	qkm2[1] -= 64;
	qkm1[1] -= 64;
	}
	}
while( ((int) t[1] - (int) r[1]) > -70 );

qexp( x, qn );
qdiv( qn, ans, ans );	/* ans *= exp( -x ) */

done:
qmov( ans, yy );
return 0;

overf:
mtherr( "qexpn", OVERFLOW );
return 0;
}
