expint.c

math library from gnu

  const double xmaxt = -GSL_LOG_DBL_MIN;      /* XMAXT = -LOG (R1MACH(1)) */
  const double xmax  = xmaxt - log(xmaxt);    /* XMAX = XMAXT - LOG(XMAXT) */

  /* CHECK_POINTER(result) */

  if(x < -xmax && !scale) {
      OVERFLOW_ERROR(result);
  }
  else if(x <= -10.0) {
    const double s = 1.0/x * ( scale ? 1.0 : exp(-x) );
    gsl_sf_result result_c;
    cheb_eval_e(&AE11_cs, 20.0/x+1.0, &result_c);
    result->val  = s * (1.0 + result_c.val);
    result->err  = s * result_c.err;
    result->err += 2.0 * GSL_DBL_EPSILON * (fabs(x) + 1.0) * fabs(result->val);
    return GSL_SUCCESS;
  }
  else if(x <= -4.0) {
    const double s = 1.0/x * ( scale ? 1.0 : exp(-x) );
    gsl_sf_result result_c;
    cheb_eval_e(&AE12_cs, (40.0/x+7.0)/3.0, &result_c);
    result->val  = s * (1.0 + result_c.val);
    result->err  = s * result_c.err;
    result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
    return GSL_SUCCESS;
  }
  else if(x <= -1.0) {
    const double ln_term = -log(fabs(x));
    const double scale_factor = ( scale ? exp(x) : 1.0 );
    gsl_sf_result result_c;
    cheb_eval_e(&E11_cs, (2.0*x+5.0)/3.0, &result_c);
    result->val  = scale_factor * (ln_term + result_c.val);
    result->err  = scale_factor * (result_c.err + GSL_DBL_EPSILON * fabs(ln_term));
    result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
    return GSL_SUCCESS;
  }
  else if(x == 0.0) {
    DOMAIN_ERROR(result);
  }
  else if(x <= 1.0) {
    const double ln_term = -log(fabs(x));
    const double scale_factor = ( scale ? exp(x) : 1.0 );
    gsl_sf_result result_c;
    cheb_eval_e(&E12_cs, x, &result_c);
    result->val  = scale_factor * (ln_term - 0.6875 + x + result_c.val);
    result->err  = scale_factor * (result_c.err + GSL_DBL_EPSILON * fabs(ln_term));
    result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
    return GSL_SUCCESS;
  }
  else if(x <= 4.0) {
    const double s = 1.0/x * ( scale ? 1.0 : exp(-x) );
    gsl_sf_result result_c;
    cheb_eval_e(&AE13_cs, (8.0/x-5.0)/3.0, &result_c);
    result->val  = s * (1.0 + result_c.val);
    result->err  = s * result_c.err;
    result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
    return GSL_SUCCESS;
  }
  else if(x <= xmax || scale) {
    const double s = 1.0/x * ( scale ? 1.0 : exp(-x) );
    gsl_sf_result result_c;
    cheb_eval_e(&AE14_cs, 8.0/x-1.0, &result_c);
    result->val  = s * (1.0 +  result_c.val);
    result->err  = s * (GSL_DBL_EPSILON + result_c.err);
    result->err += 2.0 * (x + 1.0) * GSL_DBL_EPSILON * fabs(result->val);
    if(result->val == 0.0)
      UNDERFLOW_ERROR(result);
    else
      return GSL_SUCCESS;
  }
  else {
    UNDERFLOW_ERROR(result);
  }
}


static
int expint_E2_impl(const double x, gsl_sf_result * result, const int scale)
{
  const double xmaxt = -GSL_LOG_DBL_MIN;
  const double xmax  = xmaxt - log(xmaxt);

  /* CHECK_POINTER(result) */

  if(x < -xmax && !scale) {
    OVERFLOW_ERROR(result);
  }
  else if (x == 0.0) {
    result->val = (scale ? 1.0 : 1.0);
    result->err = 0.0;
    return GSL_SUCCESS;
  } else if(x < 100.0) {
    const double ex = ( scale ? 1.0 : exp(-x) );
    gsl_sf_result result_E1;
    int stat_E1 = expint_E1_impl(x, &result_E1, scale);
    result->val  = ex - x*result_E1.val;
    result->err  = GSL_DBL_EPSILON*ex + fabs(x) * result_E1.err;
    result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
    return stat_E1;
  }
  else if(x < xmax || scale) {
    const double s = ( scale ? 1.0 : exp(-x) );
    const double c1  = -2.0;
    const double c2  =  6.0;
    const double c3  = -24.0;
    const double c4  =  120.0;
    const double c5  = -720.0;
    const double c6  =  5040.0;
    const double c7  = -40320.0;
    const double c8  =  362880.0;
    const double c9  = -3628800.0;
    const double c10 =  39916800.0;
    const double c11 = -479001600.0;
    const double c12 =  6227020800.0;
    const double c13 = -87178291200.0;
    const double y = 1.0/x;
    const double sum6 = c6+y*(c7+y*(c8+y*(c9+y*(c10+y*(c11+y*(c12+y*c13))))));
    const double sum  = y*(c1+y*(c2+y*(c3+y*(c4+y*(c5+y*sum6)))));
    result->val = s * (1.0 + sum)/x;
    result->err = 2.0 * (x + 1.0) * GSL_DBL_EPSILON * result->val;
    if(result->val == 0.0)
      UNDERFLOW_ERROR(result);
    else
      return GSL_SUCCESS;
  }
  else {
    UNDERFLOW_ERROR(result);
  }
}

static
int expint_En_impl(const int n, const double x, gsl_sf_result * result, const int scale)
{
  if (n < 0) {
    DOMAIN_ERROR(result);
  } else if (n == 0) {
    if (x == 0) {
      DOMAIN_ERROR(result);
    } else {
      result->val = (scale ? 1.0 : exp(-x)) / x;
      result->err = 2 * GSL_DBL_EPSILON * fabs(result->val);
      CHECK_UNDERFLOW(result);
      return GSL_SUCCESS;
    }
  } else if (n == 1) {
    return expint_E1_impl(x, result, scale);
  } else if (n == 2) {
    return expint_E2_impl(x, result, scale);
  } else { 
    if(x < 0) {
      DOMAIN_ERROR(result);
    }
    if (x == 0) {
      result->val = (scale ? exp(x) : 1 ) * (1/(n-1.0));
      result->err = 2 * GSL_DBL_EPSILON * fabs(result->val);
      CHECK_UNDERFLOW(result);
      return GSL_SUCCESS;
    } else {
      gsl_sf_result result_g;
      double prefactor = pow(x, n-1);
      int status = gsl_sf_gamma_inc_e (1-n, x, &result_g);
      double scale_factor = ( scale ? exp(x) : 1.0 );
      result->val = scale_factor * prefactor * result_g.val;
      result->err = 2 * GSL_DBL_EPSILON * fabs(result->val);
      result->err += 2 * fabs(scale_factor * prefactor * result_g.err);
      if (status == GSL_SUCCESS) CHECK_UNDERFLOW(result);
      return status;
    }
  }
}

/*-*-*-*-*-*-*-*-*-*-*-* Functions with Error Codes *-*-*-*-*-*-*-*-*-*-*-*/


int gsl_sf_expint_E1_e(const double x, gsl_sf_result * result)
{
  return expint_E1_impl(x, result, 0);
}


int gsl_sf_expint_E1_scaled_e(const double x, gsl_sf_result * result)
{
  return expint_E1_impl(x, result, 1);
}


int gsl_sf_expint_E2_e(const double x, gsl_sf_result * result)
{
  return expint_E2_impl(x, result, 0);
}


int gsl_sf_expint_E2_scaled_e(const double x, gsl_sf_result * result)
{
  return expint_E2_impl(x, result, 1);
}

int gsl_sf_expint_En_e(const int n, const double x, gsl_sf_result * result)
{
  return expint_En_impl(n, x, result, 0);
}


int gsl_sf_expint_En_scaled_e(const int n, const double x, gsl_sf_result * result)
{
  return expint_En_impl(n, x, result, 1);
}


int gsl_sf_expint_Ei_e(const double x, gsl_sf_result * result)
{
  /* CHECK_POINTER(result) */

  {
    int status = gsl_sf_expint_E1_e(-x, result);
    result->val = -result->val;
    return status;
  }
}


int gsl_sf_expint_Ei_scaled_e(const double x, gsl_sf_result * result)
{
  /* CHECK_POINTER(result) */

  {
    int status = gsl_sf_expint_E1_scaled_e(-x, result);
    result->val = -result->val;
    return status;
  }
}


#if 0
static double recurse_En(int n, double x, double E1)
{
  int i;
  double En = E1;
  double ex = exp(-x);
  for(i=2; i<=n; i++) {
    En = 1./(i-1) * (ex - x * En);
  }
  return En;
}
#endif


/*-*-*-*-*-*-*-*-*-* Functions w/ Natural Prototypes *-*-*-*-*-*-*-*-*-*-*/

#include "eval.h"

double gsl_sf_expint_E1(const double x)
{
  EVAL_RESULT(gsl_sf_expint_E1_e(x, &result));
}

double gsl_sf_expint_E1_scaled(const double x)
{
  EVAL_RESULT(gsl_sf_expint_E1_scaled_e(x, &result));
}

double gsl_sf_expint_E2(const double x)
{
  EVAL_RESULT(gsl_sf_expint_E2_e(x, &result));
}

double gsl_sf_expint_E2_scaled(const double x)
{
  EVAL_RESULT(gsl_sf_expint_E2_scaled_e(x, &result));
}

double gsl_sf_expint_En(const int n, const double x)
{
  EVAL_RESULT(gsl_sf_expint_En_e(n, x, &result));
}

double gsl_sf_expint_En_scaled(const int n, const double x)
{
  EVAL_RESULT(gsl_sf_expint_En_scaled_e(n, x, &result));
}

double gsl_sf_expint_Ei(const double x)
{
  EVAL_RESULT(gsl_sf_expint_Ei_e(x, &result));
}

double gsl_sf_expint_Ei_scaled(const double x)
{
  EVAL_RESULT(gsl_sf_expint_Ei_scaled_e(x, &result));
}