/*							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 "mconf.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);
}