gamma.c

开放gsl矩阵运算

int
gsl_sf_gamma_e(const double x, gsl_sf_result * result)
{
  if(x < 0.5) {
    int rint_x = (int)floor(x+0.5);
    double f_x = x - rint_x;
    double sgn = ( GSL_IS_EVEN(rint_x) ? 1.0 : -1.0 );
    double sin_term = sgn * sin(M_PI * f_x) / M_PI;

    if(sin_term == 0.0) {
      DOMAIN_ERROR(result);
    }
    else if(x > -169.0) {
      gsl_sf_result g;
      gamma_xgthalf(1.0-x, &g);
      if(fabs(sin_term) * g.val * GSL_DBL_MIN < 1.0) {
        result->val  = 1.0/(sin_term * g.val);
	result->err  = fabs(g.err/g.val) * fabs(result->val);
	result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
        return GSL_SUCCESS;
      }
      else {
        UNDERFLOW_ERROR(result);
      }
    }
    else {
      /* It is hard to control it here.
       * We can only exponentiate the
       * logarithm and eat the loss of
       * precision.
       */
      gsl_sf_result lng;
      double sgn;
      int stat_lng = gsl_sf_lngamma_sgn_e(x, &lng, &sgn);
      int stat_e   = gsl_sf_exp_mult_err_e(lng.val, lng.err, sgn, 0.0, result);
      return GSL_ERROR_SELECT_2(stat_e, stat_lng);
    }
  }
  else {
    return gamma_xgthalf(x, result);
  }
}


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

  if(x <= 0.0) {
    DOMAIN_ERROR(result);
  }
  else if(x < 0.5) {
    gsl_sf_result lg;
    const int stat_lg = gsl_sf_lngamma_e(x, &lg);
    const double lx = log(x);
    const double c  = 0.5*(M_LN2+M_LNPI);
    const double lnr_val = lg.val - (x-0.5)*lx + x - c;
    const double lnr_err = lg.err + 2.0 * GSL_DBL_EPSILON *((x+0.5)*fabs(lx) + c);
    const int stat_e  = gsl_sf_exp_err_e(lnr_val, lnr_err, result);
    return GSL_ERROR_SELECT_2(stat_lg, stat_e);
  }
  else if(x < 2.0) {
    const double t = 4.0/3.0*(x-0.5) - 1.0;
    return cheb_eval_e(&gstar_a_cs, t, result);
  }
  else if(x < 10.0) {
    const double t = 0.25*(x-2.0) - 1.0;
    gsl_sf_result c;
    cheb_eval_e(&gstar_b_cs, t, &c);
    result->val  = c.val/(x*x) + 1.0 + 1.0/(12.0*x);
    result->err  = c.err/(x*x);
    result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
    return GSL_SUCCESS;
  }
  else if(x < 1.0/GSL_ROOT4_DBL_EPSILON) {
    return gammastar_ser(x, result);
  }
  else if(x < 1.0/GSL_DBL_EPSILON) {
    /* Use Stirling formula for Gamma(x).
     */
    const double xi = 1.0/x;
    result->val = 1.0 + xi/12.0*(1.0 + xi/24.0*(1.0 - xi*(139.0/180.0 + 571.0/8640.0*xi)));
    result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
    return GSL_SUCCESS;
  }
  else {
    result->val = 1.0;
    result->err = 1.0/x;
    return GSL_SUCCESS;
  }
}


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

  if(x < 0.5) {
    gsl_sf_result lng;
    double sgn;
    int stat_lng = gsl_sf_lngamma_sgn_e(x, &lng, &sgn);
    if(stat_lng == GSL_EDOM) {
      result->val = 0.0;
      result->err = 0.0;
      return GSL_SUCCESS;
    }
    else if(stat_lng != GSL_SUCCESS) {
      result->val = 0.0;
      result->err = 0.0;
      return stat_lng;
    }
    else {
      return gsl_sf_exp_mult_err_e(-lng.val, lng.err, sgn, 0.0, result);
    }
  }
  else {
    gsl_sf_result g;
    int stat_g = gamma_xgthalf(x, &g);
    if(stat_g == GSL_EOVRFLW) {
      UNDERFLOW_ERROR(result);
    }
    else {
      result->val  = 1.0/g.val;
      result->err  = fabs(g.err/g.val) * fabs(result->val);
      result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
      CHECK_UNDERFLOW(result);
      return GSL_SUCCESS;
    }
  }
}


int
gsl_sf_lngamma_complex_e(double zr, double zi, gsl_sf_result * lnr, gsl_sf_result * arg)
{
  if(zr <= 0.5) {
    /* Transform to right half plane using reflection;
     * in fact we do a little better by stopping at 1/2.
     */
    double x = 1.0-zr;
    double y = -zi;
    gsl_sf_result a, b;
    gsl_sf_result lnsin_r, lnsin_i;

    int stat_l = lngamma_lanczos_complex(x, y, &a, &b);
    int stat_s = gsl_sf_complex_logsin_e(M_PI*zr, M_PI*zi, &lnsin_r, &lnsin_i);

    if(stat_s == GSL_SUCCESS) {
      int stat_r;
      lnr->val = M_LNPI - lnsin_r.val - a.val;
      lnr->err = lnsin_r.err + a.err + 2.0 * GSL_DBL_EPSILON * fabs(lnr->val);
      arg->val = -lnsin_i.val - b.val;
      arg->err = lnsin_i.err + b.err + 2.0 * GSL_DBL_EPSILON * fabs(arg->val);
      stat_r = gsl_sf_angle_restrict_symm_e(&(arg->val));
      return GSL_ERROR_SELECT_2(stat_r, stat_l);
    }
    else {
      DOMAIN_ERROR_2(lnr,arg);
    }
  }
  else {
    /* otherwise plain vanilla Lanczos */
    return lngamma_lanczos_complex(zr, zi, lnr, arg);
  }
}


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

  if(x < 0.0 || n < 0) {
    DOMAIN_ERROR(result);
  }
  else if(n == 0) {
    result->val = 1.0;
    result->err = 0.0;
    return GSL_SUCCESS;
  }
  else if(n == 1) {
    result->val = x;
    result->err = 0.0;
    return GSL_SUCCESS;
  }
  else if(x == 0.0) {
    result->val = 0.0;
    result->err = 0.0;
    return GSL_SUCCESS;
  }
  else {
    const double log2pi = M_LNPI + M_LN2;
    const double ln_test = n*(log(x)+1.0) + 1.0 - (n+0.5)*log(n+1.0) + 0.5*log2pi;

    if(ln_test < GSL_LOG_DBL_MIN+1.0) {
      UNDERFLOW_ERROR(result);
    }
    else if(ln_test > GSL_LOG_DBL_MAX-1.0) {
      OVERFLOW_ERROR(result);
    }
    else {
      double product = 1.0;
      int k;
      for(k=1; k<=n; k++) {
        product *= (x/k);
      }
      result->val = product;
      result->err = n * GSL_DBL_EPSILON * product;
      CHECK_UNDERFLOW(result);
      return GSL_SUCCESS;
    }    
  }
}


int gsl_sf_fact_e(const unsigned int n, gsl_sf_result * result)
{
  /* CHECK_POINTER(result) */

  if(n < 18) {
    result->val = fact_table[n].f;
    result->err = 0.0;
    return GSL_SUCCESS;
  }
  else if(n <= FACT_TABLE_MAX){
    result->val = fact_table[n].f;
    result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
    return GSL_SUCCESS;
  }
  else {
    OVERFLOW_ERROR(result);
  }
}


int gsl_sf_doublefact_e(const unsigned int n, gsl_sf_result * result)
{
  /* CHECK_POINTER(result) */

  if(n < 26) {
    result->val = doub_fact_table[n].f;
    result->err = 0.0;
    return GSL_SUCCESS;
  }
  else if(n <= DOUB_FACT_TABLE_MAX){
    result->val = doub_fact_table[n].f;
    result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
    return GSL_SUCCESS;
  }
  else {
    OVERFLOW_ERROR(result);
  }
}


int gsl_sf_lnfact_e(const unsigned int n, gsl_sf_result * result)
{
  /* CHECK_POINTER(result) */

  if(n <= FACT_TABLE_MAX){
    result->val = log(fact_table[n].f);
    result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
    return GSL_SUCCESS;
  }
  else {
    gsl_sf_lngamma_e(n+1.0, result);
    return GSL_SUCCESS;
  }
}


int gsl_sf_lndoublefact_e(const unsigned int n, gsl_sf_result * result)
{
  /* CHECK_POINTER(result) */

  if(n <= DOUB_FACT_TABLE_MAX){
    result->val = log(doub_fact_table[n].f);
    result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
    return GSL_SUCCESS;
  }
  else if(GSL_IS_ODD(n)) {
    gsl_sf_result lg;
    gsl_sf_lngamma_e(0.5*(n+2.0), &lg);
    result->val = 0.5*(n+1.0) * M_LN2 - 0.5*M_LNPI + lg.val;
    result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val) + lg.err;
    return GSL_SUCCESS;
  }
  else {
    gsl_sf_result lg;
    gsl_sf_lngamma_e(0.5*n+1.0, &lg);
    result->val = 0.5*n*M_LN2 + lg.val;
    result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val) + lg.err;
    return GSL_SUCCESS;
  }
}


int gsl_sf_lnchoose_e(unsigned int n, unsigned int m, gsl_sf_result * result)
{
  /* CHECK_POINTER(result) */

  if(m > n) {
    DOMAIN_ERROR(result);
  }
  else if(m == n || m == 0) {
    result->val = 0.0;
    result->err = 0.0;
    return GSL_SUCCESS;
  }
  else {
    gsl_sf_result nf;
    gsl_sf_result mf;
    gsl_sf_result nmmf;
    if(m*2 > n) m = n-m;
    gsl_sf_lnfact_e(n, &nf);
    gsl_sf_lnfact_e(m, &mf);
    gsl_sf_lnfact_e(n-m, &nmmf);
    result->val  = nf.val - mf.val - nmmf.val;
    result->err  = nf.err + mf.err + nmmf.err;
    result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
    return GSL_SUCCESS;
  }
}


int gsl_sf_choose_e(unsigned int n, unsigned int m, gsl_sf_result * result)
{
  if(m > n) {
    DOMAIN_ERROR(result);
  }
  else if(m == n || m == 0) {
    result->val = 1.0;
    result->err = 0.0;
    return GSL_SUCCESS;
  }
  else {
    double prod = 1.0;
    int k;
    for(k=n; k>=m+1; k--) {
      double tk = (double)k / (double)(k-m);
      if(tk > GSL_DBL_MAX/prod) {
        OVERFLOW_ERROR(result);
      }
      prod *= tk;
    }
    result->val = prod;
    result->err = 2.0 * GSL_DBL_EPSILON * prod * fabs(n-m);
    return GSL_SUCCESS;
  }
}


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

#include "eval.h"

double gsl_sf_fact(const unsigned int n)
{
  EVAL_RESULT(gsl_sf_fact_e(n, &result));
}

double gsl_sf_lnfact(const unsigned int n)
{
  EVAL_RESULT(gsl_sf_lnfact_e(n, &result));
}

double gsl_sf_doublefact(const unsigned int n)
{
  EVAL_RESULT(gsl_sf_doublefact_e(n, &result));
}

double gsl_sf_lndoublefact(const unsigned int n)
{
  EVAL_RESULT(gsl_sf_lndoublefact_e(n, &result));
}

double gsl_sf_lngamma(const double x)
{
  EVAL_RESULT(gsl_sf_lngamma_e(x, &result));
}

double gsl_sf_gamma(const double x)
{
  EVAL_RESULT(gsl_sf_gamma_e(x, &result));
}

double gsl_sf_gammastar(const double x)
{
  EVAL_RESULT(gsl_sf_gammastar_e(x, &result));
}

double gsl_sf_gammainv(const double x)
{
  EVAL_RESULT(gsl_sf_gammainv_e(x, &result));
}

double gsl_sf_taylorcoeff(const int n, const double x)
{
  EVAL_RESULT(gsl_sf_taylorcoeff_e(n, x, &result));
}

double gsl_sf_choose(unsigned int n, unsigned int m)
{
  EVAL_RESULT(gsl_sf_choose_e(n, m, &result));
}

double gsl_sf_lnchoose(unsigned int n, unsigned int m)
{
  EVAL_RESULT(gsl_sf_lnchoose_e(n, m, &result));
}