fermi_dirac.c

开放gsl矩阵运算

  /* Now sum the terms where eta() is negative.
   * The argument of eta() must be odd as well,
   * so it is convenient to transform the series
   * as follows:
   *
   * Sum[ eta(j+1-n) x^n / n!, {n,j+4,Infinity}]
   * = x^j / j! Sum[ eta(1-2m) x^(2m) j! / (2m+j)! , {m,2,Infinity}]
   *
   * We do not need to do this sum if j is large enough.
   */
  if(j < 32) {
    int m;
    gsl_sf_result jfact;
    double sum2;
    double pre2;

    gsl_sf_fact_e((unsigned int)j, &jfact);
    pre2 = gsl_sf_pow_int(x, j) / jfact.val;

    gsl_sf_eta_int_e(-3, &eta_factor);
    pow_factor = x*x*x*x / ((j+4)*(j+3)*(j+2)*(j+1));
    sum2 = eta_factor.val * pow_factor;

    for(m=3; m<24; m++) {
      gsl_sf_eta_int_e(1-2*m, &eta_factor);
      pow_factor *= x*x / ((j+2*m)*(j+2*m-1));
      sum2 += eta_factor.val * pow_factor;
    }

    sum += pre2 * sum2;
  }

  result->val = sum;
  result->err = 2.0 * GSL_DBL_EPSILON * fabs(sum);

  return GSL_SUCCESS;
}


/* series of hypergeometric functions for integer j > 0, x > 0
 * [Goano (7)]
 */
static
int
fd_UMseries_int(const int j, const double x, gsl_sf_result * result)
{
  const int nmax = 2000;
  double pre;
  double lnpre_val;
  double lnpre_err;
  double sum_even_val = 1.0;
  double sum_even_err = 0.0;
  double sum_odd_val  = 0.0;
  double sum_odd_err  = 0.0;
  int stat_sum;
  int stat_e;
  int stat_h = GSL_SUCCESS;
  int n;

  if(x < 500.0 && j < 80) {
    double p = gsl_sf_pow_int(x, j+1);
    gsl_sf_result g;
    gsl_sf_fact_e(j+1, &g); /* Gamma(j+2) */
    lnpre_val = 0.0;
    lnpre_err = 0.0;
    pre   = p/g.val;
  }
  else {
    double lnx = log(x);
    gsl_sf_result lg;
    gsl_sf_lngamma_e(j + 2.0, &lg);
    lnpre_val = (j+1.0)*lnx - lg.val;
    lnpre_err = 2.0 * GSL_DBL_EPSILON * fabs((j+1.0)*lnx) + lg.err;
    pre = 1.0;
  }

  /* Add up the odd terms of the sum.
   */
  for(n=1; n<nmax; n+=2) {
    double del_val;
    double del_err;
    gsl_sf_result U;
    gsl_sf_result M;
    int stat_h_U = gsl_sf_hyperg_U_int_e(1, j+2, n*x, &U);
    int stat_h_F = gsl_sf_hyperg_1F1_int_e(1, j+2, -n*x, &M);
    stat_h = GSL_ERROR_SELECT_3(stat_h, stat_h_U, stat_h_F);
    del_val = ((j+1.0)*U.val - M.val);
    del_err = (fabs(j+1.0)*U.err + M.err);
    sum_odd_val += del_val;
    sum_odd_err += del_err;
    if(fabs(del_val/sum_odd_val) < GSL_DBL_EPSILON) break;
  }

  /* Add up the even terms of the sum.
   */
  for(n=2; n<nmax; n+=2) {
    double del_val;
    double del_err;
    gsl_sf_result U;
    gsl_sf_result M;
    int stat_h_U = gsl_sf_hyperg_U_int_e(1, j+2, n*x, &U);
    int stat_h_F = gsl_sf_hyperg_1F1_int_e(1, j+2, -n*x, &M);
    stat_h = GSL_ERROR_SELECT_3(stat_h, stat_h_U, stat_h_F);
    del_val = ((j+1.0)*U.val - M.val);
    del_err = (fabs(j+1.0)*U.err + M.err);
    sum_even_val -= del_val;
    sum_even_err += del_err;
    if(fabs(del_val/sum_even_val) < GSL_DBL_EPSILON) break;
  }

  stat_sum = ( n >= nmax ? GSL_EMAXITER : GSL_SUCCESS );
  stat_e   = gsl_sf_exp_mult_err_e(lnpre_val, lnpre_err,
                                      pre*(sum_even_val + sum_odd_val),
				      pre*(sum_even_err + sum_odd_err),
				      result);
  result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);

  return GSL_ERROR_SELECT_3(stat_e, stat_h, stat_sum);
}


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

/* [Goano (4)] */
int gsl_sf_fermi_dirac_m1_e(const double x, gsl_sf_result * result)
{
  if(x < GSL_LOG_DBL_MIN) {
    UNDERFLOW_ERROR(result);
  }
  else if(x < 0.0) {
    const double ex = exp(x);
    result->val = ex/(1.0+ex);
    result->err = 2.0 * (fabs(x) + 1.0) * GSL_DBL_EPSILON * fabs(result->val);
    return GSL_SUCCESS;
  }
  else {
    double ex = exp(-x);
    result->val = 1.0/(1.0 + ex);
    result->err = 2.0 * GSL_DBL_EPSILON * (x + 1.0) * ex;
    return GSL_SUCCESS;
  }
}


/* [Goano (3)] */
int gsl_sf_fermi_dirac_0_e(const double x, gsl_sf_result * result)
{
  if(x < GSL_LOG_DBL_MIN) {
    UNDERFLOW_ERROR(result);
  }
  else if(x < -5.0) {
    double ex  = exp(x);
    double ser = 1.0 - ex*(0.5 - ex*(1.0/3.0 - ex*(1.0/4.0 - ex*(1.0/5.0 - ex/6.0))));
    result->val = ex * ser;
    result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
    return GSL_SUCCESS;
  }
  else if(x < 10.0) {
    result->val = log(1.0 + exp(x));
    result->err = fabs(x * GSL_DBL_EPSILON);
    return GSL_SUCCESS;
  }
  else {
    double ex = exp(-x);
    result->val = x + ex * (1.0 - 0.5*ex + ex*ex/3.0 - ex*ex*ex/4.0);
    result->err = (x + ex) * GSL_DBL_EPSILON;
    return GSL_SUCCESS;
  }
}


int gsl_sf_fermi_dirac_1_e(const double x, gsl_sf_result * result)
{
  if(x < GSL_LOG_DBL_MIN) {
    UNDERFLOW_ERROR(result);
  }
  else if(x < -1.0) {
    /* series [Goano (6)]
     */
    double ex   = exp(x);
    double term = ex;
    double sum  = term;
    int n;
    for(n=2; n<100 ; n++) {
      double rat = (n-1.0)/n;
      term *= -ex * rat * rat;
      sum  += term;
      if(fabs(term/sum) < GSL_DBL_EPSILON) break;
    }
    result->val = sum;
    result->err = 2.0 * fabs(sum) * GSL_DBL_EPSILON;
    return GSL_SUCCESS;
  }
  else if(x < 1.0) {
    return cheb_eval_e(&fd_1_a_cs, x, result);
  }
  else if(x < 4.0) {
    double t = 2.0/3.0*(x-1.0) - 1.0;
    return cheb_eval_e(&fd_1_b_cs, t, result);
  }
  else if(x < 10.0) {
    double t = 1.0/3.0*(x-4.0) - 1.0;
    return cheb_eval_e(&fd_1_c_cs, t, result);
  }
  else if(x < 30.0) {
    double t = 0.1*x - 2.0;
    gsl_sf_result c;
    cheb_eval_e(&fd_1_d_cs, t, &c);
    result->val  = c.val * x*x;
    result->err  = c.err * x*x + GSL_DBL_EPSILON * fabs(result->val);
    return GSL_SUCCESS;
  }
  else if(x < 1.0/GSL_SQRT_DBL_EPSILON) {
    double t = 60.0/x - 1.0;
    gsl_sf_result c;
    cheb_eval_e(&fd_1_e_cs, t, &c);
    result->val  = c.val * x*x;
    result->err  = c.err * x*x + GSL_DBL_EPSILON * fabs(result->val);
    return GSL_SUCCESS;
  }
  else if(x < GSL_SQRT_DBL_MAX) {
    result->val = 0.5 * x*x;
    result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
    return GSL_SUCCESS;
  }
  else {
    OVERFLOW_ERROR(result);
  }
}


int gsl_sf_fermi_dirac_2_e(const double x, gsl_sf_result * result)
{
  if(x < GSL_LOG_DBL_MIN) {
    UNDERFLOW_ERROR(result);
  }
  else if(x < -1.0) {
    /* series [Goano (6)]
     */
    double ex   = exp(x);
    double term = ex;
    double sum  = term;
    int n;
    for(n=2; n<100 ; n++) {
      double rat = (n-1.0)/n;
      term *= -ex * rat * rat * rat;
      sum  += term;
      if(fabs(term/sum) < GSL_DBL_EPSILON) break;
    }
    result->val = sum;
    result->err = 2.0 * GSL_DBL_EPSILON * fabs(sum);
    return GSL_SUCCESS;
  }
  else if(x < 1.0) {
    return cheb_eval_e(&fd_2_a_cs, x, result);
  }
  else if(x < 4.0) {
    double t = 2.0/3.0*(x-1.0) - 1.0;
    return cheb_eval_e(&fd_2_b_cs, t, result);
  }
  else if(x < 10.0) {
    double t = 1.0/3.0*(x-4.0) - 1.0;
    return cheb_eval_e(&fd_2_c_cs, t, result);
  }
  else if(x < 30.0) {
    double t = 0.1*x - 2.0;
    gsl_sf_result c;
    cheb_eval_e(&fd_2_d_cs, t, &c);
    result->val  = c.val * x*x*x;
    result->err  = c.err * x*x*x + 3.0 * GSL_DBL_EPSILON * fabs(result->val);
    return GSL_SUCCESS;
  }
  else if(x < 1.0/GSL_ROOT3_DBL_EPSILON) {
    double t = 60.0/x - 1.0;
    gsl_sf_result c;
    cheb_eval_e(&fd_2_e_cs, t, &c);
    result->val = c.val * x*x*x;
    result->err = c.err * x*x*x + 3.0 * GSL_DBL_EPSILON * fabs(result->val);
    return GSL_SUCCESS;
  }
  else if(x < GSL_ROOT3_DBL_MAX) {
    result->val = 1.0/6.0 * x*x*x;
    result->err = 3.0 * GSL_DBL_EPSILON * fabs(result->val);
    return GSL_SUCCESS;
  }
  else {
    OVERFLOW_ERROR(result);
  }
}


int gsl_sf_fermi_dirac_int_e(const int j, const double x, gsl_sf_result * result)
{
  if(j < -1) {
    return fd_nint(j, x, result);
  }
  else if (j == -1) {
    return gsl_sf_fermi_dirac_m1_e(x, result);
  }
  else if(j == 0) {
    return gsl_sf_fermi_dirac_0_e(x, result);
  }
  else if(j == 1) {
    return gsl_sf_fermi_dirac_1_e(x, result);
  }
  else if(j == 2) {
    return gsl_sf_fermi_dirac_2_e(x, result);
  }
  else if(x < 0.0) {
    return fd_neg(j, x, result);
  }
  else if(x == 0.0) {
    return gsl_sf_eta_int_e(j+1, result);
  }
  else if(x < 1.5) {
    return fd_series_int(j, x, result);
  }
  else {
    gsl_sf_result fasymp;
    int stat_asymp = fd_asymp(j, x, &fasymp);

    if(stat_asymp == GSL_SUCCESS) {
      result->val  = fasymp.val;
      result->err  = fasymp.err;
      result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
      return stat_asymp;
    }
    else {
      return fd_UMseries_int(j, x, result);
    }
  }
}


int gsl_sf_fermi_dirac_mhalf_e(const double x, gsl_sf_result * result)
{
  if(x < GSL_LOG_DBL_MIN) {
    UNDERFLOW_ERROR(result);
  }
  else if(x < -1.0) {
    /* series [Goano (6)]
     */
    double ex   = exp(x);
    double term = ex;
    double sum  = term;
    int n;
    for(n=2; n<200 ; n++) {
      double rat = (n-1.0)/n;
      term *= -ex * sqrt(rat);
      sum  += term;
      if(fabs(term/sum) < GSL_DBL_EPSILON) break;
    }
    result->val = sum;
    result->err = 2.0 * fabs(sum) * GSL_DBL_EPSILON;
    return GSL_SUCCESS;
  }
  else if(x < 1.0) {
    return cheb_eval_e(&fd_mhalf_a_cs, x, result);
  }
  else if(x < 4.0) {
    double t = 2.0/3.0*(x-1.0) - 1.0;
    return cheb_eval_e(&fd_mhalf_b_cs, t, result);
  }
  else if(x < 10.0) {
    double t = 1.0/3.0*(x-4.0) - 1.0;
    return cheb_eval_e(&fd_mhalf_c_cs, t, result);
  }
  else if(x < 30.0) {
    double rtx = sqrt(x);
    double t = 0.1*x - 2.0;
    gsl_sf_result c;
    cheb_eval_e(&fd_mhalf_d_cs, t, &c);
    result->val  = c.val * rtx;
    result->err  = c.err * rtx + 0.5 * GSL_DBL_EPSILON * fabs(result->val);
    return GSL_SUCCESS;
  }
  else {
    return fd_asymp(-0.5, x, result);
  }
}


int gsl_sf_fermi_dirac_half_e(const double x, gsl_sf_result * result)
{
  if(x < GSL_LOG_DBL_MIN) {
    UNDERFLOW_ERROR(result);
  }
  else if(x < -1.0) {
    /* series [Goano (6)]
     */
    double ex   = exp(x);
    double term = ex;
    double sum  = term;
    int n;
    for(n=2; n<100 ; n++) {
      double rat = (n-1.0)/n;
      term *= -ex * rat * sqrt(rat);
      sum  += term;
      if(fabs(term/sum) < GSL_DBL_EPSILON) break;
    }
    result->val = sum;
    result->err = 2.0 * fabs(sum) * GSL_DBL_EPSILON;
    return GSL_SUCCESS;
  }
  else if(x < 1.0) {
    return cheb_eval_e(&fd_half_a_cs, x, result);
  }
  else if(x < 4.0) {
    double t = 2.0/3.0*(x-1.0) - 1.0;
    return cheb_eval_e(&fd_half_b_cs, t, result);
  }
  else if(x < 10.0) {
    double t = 1.0/3.0*(x-4.0) - 1.0;
    return cheb_eval_e(&fd_half_c_cs, t, result);
  }
  else if(x < 30.0) {
    double x32 = x*sqrt(x);
    double t = 0.1*x - 2.0;
    gsl_sf_result c;
    cheb_eval_e(&fd_half_d_cs, t, &c);
    result->val = c.val * x32;
    result->err = c.err * x32 + 1.5 * GSL_DBL_EPSILON * fabs(result->val);
    return GSL_SUCCESS;
  }
  else {
    return fd_asymp(0.5, x, result);
  }
}


int gsl_sf_fermi_dirac_3half_e(const double x, gsl_sf_result * result)
{
  if(x < GSL_LOG_DBL_MIN) {
    UNDERFLOW_ERROR(result);
  }
  else if(x < -1.0) {
    /* series [Goano (6)]
     */
    double ex   = exp(x);
    double term = ex;
    double sum  = term;
    int n;
    for(n=2; n<100 ; n++) {
      double rat = (n-1.0)/n;
      term *= -ex * rat * rat * sqrt(rat);
      sum  += term;
      if(fabs(term/sum) < GSL_DBL_EPSILON) break;
    }
    result->val = sum;
    result->err = 2.0 * fabs(sum) * GSL_DBL_EPSILON;
    return GSL_SUCCESS;
  }
  else if(x < 1.0) {
    return cheb_eval_e(&fd_3half_a_cs, x, result);
  }
  else if(x < 4.0) {
    double t = 2.0/3.0*(x-1.0) - 1.0;
    return cheb_eval_e(&fd_3half_b_cs, t, result);
  }
  else if(x < 10.0) {
    double t = 1.0/3.0*(x-4.0) - 1.0;
    return cheb_eval_e(&fd_3half_c_cs, t, result);
  }
  else if(x < 30.0) {
    double x52 = x*x*sqrt(x);
    double t = 0.1*x - 2.0;
    gsl_sf_result c;
    cheb_eval_e(&fd_3half_d_cs, t, &c);
    result->val = c.val * x52;
    result->err = c.err * x52 + 2.5 * GSL_DBL_EPSILON * fabs(result->val);
    return GSL_SUCCESS;
  }
  else {
    return fd_asymp(1.5, x, result);
  }
}

/* [Goano p. 222] */
int gsl_sf_fermi_dirac_inc_0_e(const double x, const double b, gsl_sf_result * result)
{
  if(b < 0.0) {
    DOMAIN_ERROR(result);
  }
  else {
    double arg = b - x;
    gsl_sf_result f0;
    int status = gsl_sf_fermi_dirac_0_e(arg, &f0);
    result->val = f0.val - arg;
    result->err = f0.err + GSL_DBL_EPSILON * (fabs(x) + fabs(b));
    return status;
  }
}



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

#include "eval.h"

double gsl_sf_fermi_dirac_m1(const double x)
{
  EVAL_RESULT(gsl_sf_fermi_dirac_m1_e(x, &result));
}

double gsl_sf_fermi_dirac_0(const double x)
{
  EVAL_RESULT(gsl_sf_fermi_dirac_0_e(x, &result));
}

double gsl_sf_fermi_dirac_1(const double x)
{
  EVAL_RESULT(gsl_sf_fermi_dirac_1_e(x, &result));
}

double gsl_sf_fermi_dirac_2(const double x)
{
  EVAL_RESULT(gsl_sf_fermi_dirac_2_e(x, &result));
}

double gsl_sf_fermi_dirac_int(const int j, const double x)
{
  EVAL_RESULT(gsl_sf_fermi_dirac_int_e(j, x, &result));
}

double gsl_sf_fermi_dirac_mhalf(const double x)
{
  EVAL_RESULT(gsl_sf_fermi_dirac_mhalf_e(x, &result));
}

double gsl_sf_fermi_dirac_half(const double x)
{
  EVAL_RESULT(gsl_sf_fermi_dirac_half_e(x, &result));
}

double gsl_sf_fermi_dirac_3half(const double x)
{
  EVAL_RESULT(gsl_sf_fermi_dirac_3half_e(x, &result));
}

double gsl_sf_fermi_dirac_inc_0(const double x, const double b)
{
  EVAL_RESULT(gsl_sf_fermi_dirac_inc_0_e(x, b, &result));
}