gamma.c

开放gsl矩阵运算

    psi_3.val = 0.0;
    psi_4.val = 0.0;
    psi_5.val = 0.0;
    psi_6.val = 0.0;
    gsl_sf_lnfact_e(N, &c0);
    gsl_sf_psi_int_e(N+1, &psi_0);
    gsl_sf_psi_1_int_e(N+1, &psi_1);
    if(aeps > 0.00001) gsl_sf_psi_n_e(2, N+1.0, &psi_2);
    if(aeps > 0.0002)  gsl_sf_psi_n_e(3, N+1.0, &psi_3);
    if(aeps > 0.001)   gsl_sf_psi_n_e(4, N+1.0, &psi_4);
    if(aeps > 0.005)   gsl_sf_psi_n_e(5, N+1.0, &psi_5);
    if(aeps > 0.01)    gsl_sf_psi_n_e(6, N+1.0, &psi_6);
    c1 = psi_0.val;
    c2 = psi_1.val/2.0;
    c3 = psi_2.val/6.0;
    c4 = psi_3.val/24.0;
    c5 = psi_4.val/120.0;
    c6 = psi_5.val/720.0;
    c7 = psi_6.val/5040.0;
    lng_ser = c0.val-eps*(c1-eps*(c2-eps*(c3-eps*(c4-eps*(c5-eps*(c6-eps*c7))))));

    /* calculate
     * g = ln(|eps gamma(-N+eps)|)
     *   = -ln(gamma(1+N-eps)) + ln(|eps Pi/sin(Pi(N+1+eps))|)
     */
    g = -lng_ser - log(sin_ser);

    lng->val = g - log(fabs(eps));
    lng->err = c0.err + 2.0 * GSL_DBL_EPSILON * (fabs(g) + fabs(lng->val));

    *sgn = ( GSL_IS_ODD(N) ? -1.0 : 1.0 ) * ( eps > 0.0 ? 1.0 : -1.0 );

    return GSL_SUCCESS;
  }
}


/* This gets bad near the negative half axis. However, this
 * region can be avoided by use of the reflection formula, as usual.
 * Only the first two terms of the series are kept.
 */
#if 0
static
int
lngamma_complex_stirling(const double zr, const double zi, double * lg_r, double * arg)
{
  double re_zinv,  im_zinv;
  double re_zinv2, im_zinv2;
  double re_zinv3, im_zinv3;
  double re_zhlnz, im_zhlnz;
  double r, lnr, theta;
  gsl_sf_complex_log_e(zr, zi, &lnr, &theta);  /* z = r e^{i theta} */
  r = exp(lnr);
  re_zinv =  (zr/r)/r;
  im_zinv = -(zi/r)/r;
  re_zinv2 = re_zinv*re_zinv - im_zinv*im_zinv;
  re_zinv2 = 2.0*re_zinv*im_zinv;
  re_zinv3 = re_zinv2*re_zinv - im_zinv2*im_zinv;
  re_zinv3 = re_zinv2*im_zinv + im_zinv2*re_zinv;
  re_zhlnz = (zr - 0.5)*lnr - zi*theta;
  im_zhlnz = zi*lnr + zr*theta;
  *lg_r = re_zhlnz - zr + 0.5*(M_LN2+M_LNPI) + re_zinv/12.0 - re_zinv3/360.0;
  *arg  = im_zhlnz - zi + 1.0/12.0*im_zinv - im_zinv3/360.0;
  return GSL_SUCCESS;
}
#endif /* 0 */


inline
static
int
lngamma_1_pade(const double eps, gsl_sf_result * result)
{
  /* Use (2,2) Pade for Log[Gamma[1+eps]]/eps
   * plus a correction series.
   */
  const double n1 = -1.0017419282349508699871138440;
  const double n2 =  1.7364839209922879823280541733;
  const double d1 =  1.2433006018858751556055436011;
  const double d2 =  5.0456274100274010152489597514;
  const double num = (eps + n1) * (eps + n2);
  const double den = (eps + d1) * (eps + d2);
  const double pade = 2.0816265188662692474880210318 * num / den;
  const double c0 =  0.004785324257581753;
  const double c1 = -0.01192457083645441;
  const double c2 =  0.01931961413960498;
  const double c3 = -0.02594027398725020;
  const double c4 =  0.03141928755021455;
  const double eps5 = eps*eps*eps*eps*eps;
  const double corr = eps5 * (c0 + eps*(c1 + eps*(c2 + eps*(c3 + c4*eps))));
  result->val = eps * (pade + corr);
  result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
  return GSL_SUCCESS;
}

inline
static
int
lngamma_2_pade(const double eps, gsl_sf_result * result)
{
  /* Use (2,2) Pade for Log[Gamma[2+eps]]/eps
   * plus a correction series.
   */
  const double n1 = 1.000895834786669227164446568;
  const double n2 = 4.209376735287755081642901277;
  const double d1 = 2.618851904903217274682578255;
  const double d2 = 10.85766559900983515322922936;
  const double num = (eps + n1) * (eps + n2);
  const double den = (eps + d1) * (eps + d2);
  const double pade = 2.85337998765781918463568869 * num/den;
  const double c0 =  0.0001139406357036744;
  const double c1 = -0.0001365435269792533;
  const double c2 =  0.0001067287169183665;
  const double c3 = -0.0000693271800931282;
  const double c4 =  0.0000407220927867950;
  const double eps5 = eps*eps*eps*eps*eps;
  const double corr = eps5 * (c0 + eps*(c1 + eps*(c2 + eps*(c3 + c4*eps))));
  result->val = eps * (pade + corr);
  result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
  return GSL_SUCCESS;
}


/* series for gammastar(x)
 * double-precision for x > 10.0
 */
static
int
gammastar_ser(const double x, gsl_sf_result * result)
{
  /* Use the Stirling series for the correction to Log(Gamma(x)),
   * which is better behaved and easier to compute than the
   * regular Stirling series for Gamma(x). 
   */
  const double y = 1.0/(x*x);
  const double c0 =  1.0/12.0;
  const double c1 = -1.0/360.0;
  const double c2 =  1.0/1260.0;
  const double c3 = -1.0/1680.0;
  const double c4 =  1.0/1188.0;
  const double c5 = -691.0/360360.0;
  const double c6 =  1.0/156.0;
  const double c7 = -3617.0/122400.0;
  const double ser = c0 + y*(c1 + y*(c2 + y*(c3 + y*(c4 + y*(c5 + y*(c6 + y*c7))))));
  result->val = exp(ser/x);
  result->err = 2.0 * GSL_DBL_EPSILON * result->val * GSL_MAX_DBL(1.0, ser/x);
  return GSL_SUCCESS;
}


/* Chebyshev expansion for log(gamma(x)/gamma(8))
 * 5 < x < 10
 * -1 < t < 1
 */
static double gamma_5_10_data[24] = {
 -1.5285594096661578881275075214,
  4.8259152300595906319768555035,
  0.2277712320977614992970601978,
 -0.0138867665685617873604917300,
  0.0012704876495201082588139723,
 -0.0001393841240254993658962470,
  0.0000169709242992322702260663,
 -2.2108528820210580075775889168e-06,
  3.0196602854202309805163918716e-07,
 -4.2705675000079118380587357358e-08,
  6.2026423818051402794663551945e-09,
 -9.1993973208880910416311405656e-10,
  1.3875551258028145778301211638e-10,
 -2.1218861491906788718519522978e-11,
  3.2821736040381439555133562600e-12,
 -5.1260001009953791220611135264e-13,
  8.0713532554874636696982146610e-14,
 -1.2798522376569209083811628061e-14,
  2.0417711600852502310258808643e-15,
 -3.2745239502992355776882614137e-16,
  5.2759418422036579482120897453e-17,
 -8.5354147151695233960425725513e-18,
  1.3858639703888078291599886143e-18,
 -2.2574398807738626571560124396e-19
};
static const cheb_series gamma_5_10_cs = {
  gamma_5_10_data,
  23,
  -1, 1,
  11
};


/* gamma(x) for x >= 1/2
 * assumes x >= 1/2
 */
static
int
gamma_xgthalf(const double x, gsl_sf_result * result)
{
  /* CHECK_POINTER(result) */

  if(x == 0.5) {
    result->val = 1.77245385090551602729817;
    result->err = GSL_DBL_EPSILON * result->val;
    return GSL_SUCCESS;
  }
  else if(fabs(x - 1.0) < 0.01) {
    /* Use series for Gamma[1+eps] - 1/(1+eps).
     */
    const double eps = x - 1.0;
    const double c1 =  0.4227843350984671394;
    const double c2 = -0.01094400467202744461;
    const double c3 =  0.09252092391911371098;
    const double c4 = -0.018271913165599812664;
    const double c5 =  0.018004931096854797895;
    const double c6 = -0.006850885378723806846;
    const double c7 =  0.003998239557568466030;
    result->val = 1.0/x + eps*(c1+eps*(c2+eps*(c3+eps*(c4+eps*(c5+eps*(c6+eps*c7))))));
    result->err = GSL_DBL_EPSILON;
    return GSL_SUCCESS;
  }
  else if(fabs(x - 2.0) < 0.01) {
    /* Use series for Gamma[1 + eps].
     */
    const double eps = x - 2.0;
    const double c1 =  0.4227843350984671394;
    const double c2 =  0.4118403304264396948;
    const double c3 =  0.08157691924708626638;
    const double c4 =  0.07424901075351389832;
    const double c5 = -0.00026698206874501476832;
    const double c6 =  0.011154045718130991049;
    const double c7 = -0.002852645821155340816;
    const double c8 =  0.0021039333406973880085;
    result->val = 1.0 + eps*(c1+eps*(c2+eps*(c3+eps*(c4+eps*(c5+eps*(c6+eps*(c7+eps*c8)))))));
    result->err = GSL_DBL_EPSILON;
    return GSL_SUCCESS;
  }
  else if(x < 5.0) {
    /* Exponentiating the logarithm is fine, as
     * long as the exponential is not so large
     * that it greatly amplifies the error.
     */
    gsl_sf_result lg;
    lngamma_lanczos(x, &lg);
    result->val = exp(lg.val);
    result->err = result->val * (lg.err + 2.0 * GSL_DBL_EPSILON);
    return GSL_SUCCESS;
  }
  else if(x < 10.0) {
    /* This is a sticky area. The logarithm
     * is too large and the gammastar series
     * is not good.
     */
    const double gamma_8 = 5040.0;
    const double t = (2.0*x - 15.0)/5.0;
    gsl_sf_result c;
    cheb_eval_e(&gamma_5_10_cs, t, &c);
    result->val  = exp(c.val) * gamma_8;
    result->err  = result->val * c.err;
    result->err += 2.0 * GSL_DBL_EPSILON * result->val;
    return GSL_SUCCESS;
  }
  else if(x < GSL_SF_GAMMA_XMAX) {
    /* We do not want to exponentiate the logarithm
     * if x is large because of the inevitable
     * inflation of the error. So we carefully
     * use pow() and exp() with exact quantities.
     */
    double p = pow(x, 0.5*x);
    double e = exp(-x);
    double q = (p * e) * p;
    double pre = M_SQRT2 * M_SQRTPI * q/sqrt(x);
    gsl_sf_result gstar;
    int stat_gs = gammastar_ser(x, &gstar);
    result->val = pre * gstar.val;
    result->err = (x + 2.5) * GSL_DBL_EPSILON * result->val;
    return stat_gs;
  }
  else {
    OVERFLOW_ERROR(result);
  }
}


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


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

  if(fabs(x - 1.0) < 0.01) {
    /* Note that we must amplify the errors
     * from the Pade evaluations because of
     * the way we must pass the argument, i.e.
     * writing (1-x) is a loss of precision
     * when x is near 1.
     */
    int stat = lngamma_1_pade(x - 1.0, result);
    result->err *= 1.0/(GSL_DBL_EPSILON + fabs(x - 1.0));
    return stat;
  }
  else if(fabs(x - 2.0) < 0.01) {
    int stat = lngamma_2_pade(x - 2.0, result);
    result->err *= 1.0/(GSL_DBL_EPSILON + fabs(x - 2.0));
    return stat;
  }
  else if(x >= 0.5) {
    return lngamma_lanczos(x, result);
  }
  else if(x == 0.0) {
    DOMAIN_ERROR(result);
  }
  else if(fabs(x) < 0.02) {
    double sgn;
    return lngamma_sgn_0(x, result, &sgn);
  }
  else if(x > -0.5/(GSL_DBL_EPSILON*M_PI)) {
    /* Try to extract a fractional
     * part from x.
     */
    double z  = 1.0 - x;
    double s  = sin(M_PI*z);
    double as = fabs(s);
    if(s == 0.0) {
      DOMAIN_ERROR(result);
    }
    else if(as < M_PI*0.015) {
      /* x is near a negative integer, -N */
      if(x < INT_MIN + 2.0) {
        result->val = 0.0;
        result->err = 0.0;
        GSL_ERROR ("error", GSL_EROUND);
      }
      else {
        int N = -(int)(x - 0.5);
        double eps = x + N;
        double sgn;
        return lngamma_sgn_sing(N, eps, result, &sgn);
      }
    }
    else {
      gsl_sf_result lg_z;
      lngamma_lanczos(z, &lg_z);
      result->val = M_LNPI - (log(as) + lg_z.val);
      result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val) + lg_z.err;
      return GSL_SUCCESS;
    }
  }
  else {
    /* |x| was too large to extract any fractional part */
    result->val = 0.0;
    result->err = 0.0;
    GSL_ERROR ("error", GSL_EROUND);
  }
}


int gsl_sf_lngamma_sgn_e(double x, gsl_sf_result * result_lg, double * sgn)
{
  if(fabs(x - 1.0) < 0.01) {
    *sgn = 1.0;
    return lngamma_1_pade(x-1.0, result_lg);
  }
  else if(fabs(x - 2.0) < 0.01) {
    *sgn = 1.0;
    return lngamma_2_pade(x-2.0, result_lg);
  }
  else if(x >= 0.5) {
    *sgn = 1.0;
    return lngamma_lanczos(x, result_lg);
  }
  else if(x == 0.0) {
    *sgn = 0.0;
    DOMAIN_ERROR(result_lg);
  }
  else if(fabs(x) < 0.02) {
    return lngamma_sgn_0(x, result_lg, sgn);
  }
  else if(x > -0.5/(GSL_DBL_EPSILON*M_PI)) {
    double z = 1.0 - x;
    double s = sin(M_PI*x);
    double as = fabs(s);
    if(s == 0.0) {
      *sgn = 0.0;
      DOMAIN_ERROR(result_lg);
    }
    else if(as < M_PI*0.015) {
      /* x is near a negative integer, -N */
      if(x < INT_MIN + 2.0) {
        result_lg->val = 0.0;
        result_lg->err = 0.0;
	*sgn = 0.0;
	GSL_ERROR ("error", GSL_EROUND);
      }
      else {
        int N = -(int)(x - 0.5);
	double eps = x + N;
        return lngamma_sgn_sing(N, eps, result_lg, sgn);
      }
    }
    else {
      gsl_sf_result lg_z;
      lngamma_lanczos(z, &lg_z);
      *sgn = (s > 0.0 ? 1.0 : -1.0);
      result_lg->val = M_LNPI - (log(as) + lg_z.val);
      result_lg->err = 2.0 * GSL_DBL_EPSILON * fabs(result_lg->val) + lg_z.err;
      return GSL_SUCCESS;
    }
  }
  else {
    /* |x| was too large to extract any fractional part */
    result_lg->val = 0.0;
    result_lg->err = 0.0;
    *sgn = 0.0;
    GSL_ERROR ("error", GSL_EROUND);
  }
}