hyperg_1f1.c

开放gsl矩阵运算

        result->val  = Mn;
	result->err  = (x + 1.0) * GSL_DBL_EPSILON * fabs(Mn);
	result->err *= fabs(a-b)+1.0;
        return GSL_SUCCESS;
      }
      else {
        OVERFLOW_ERROR(result);
      }
    }
    else {
      /* b > a
       * b < 2a + x 
       * b <= x (otherwise we would have finished above)
       *
       * Gautschi anomalous convergence region. However, we can
       * recurse forward all the way from a=0,1 because we are
       * always underneath the b=2a+x line.
       */
      gsl_sf_result r_Mn;
      double Mnm1 = 1.0;    /* 1F1(0,b,x) */
      double Mn;	    /* 1F1(1,b,x)  */
      double Mnp1;
      int n;
      gsl_sf_exprel_n_e(b-1, x, &r_Mn);
      Mn = r_Mn.val;
      for(n=1; n<a; n++) {
        Mnp1 = ((b-n)*Mnm1 + (2*n-b+x)*Mn)/n;
     	Mnm1 = Mn;
     	Mn   = Mnp1;
      }
      result->val  = Mn;
      result->err  = fabs(Mn) * (1.0 + fabs(a)) * fabs(r_Mn.err / r_Mn.val);
      result->err += 2.0 * GSL_DBL_EPSILON * fabs(Mn);
      return GSL_SUCCESS;
    }
  }
  else {
    /* x < 0
     * b < a (otherwise we would have tripped one of the above)
     */

    if(a <= 0.5*(b-x) || a >= -x) {
      /* Gautschi continued fraction is in the anomalous region,
       * so we must find another way. We recurse down in b,
       * from the a=b line.
       */
      double ex = exp(x);
      double Manp1 = ex;
      double Man   = ex * (1.0 + x/(a-1.0));
      double Manm1;
      int n;
      for(n=a-1; n>b; n--) {
        Manm1 = (-n*(1-n-x)*Man - x*(n-a)*Manp1)/(n*(n-1.0));
        Manp1 = Man;
        Man = Manm1;
      }
      result->val  = Man;
      result->err  = (fabs(x) + 1.0) * GSL_DBL_EPSILON * fabs(Man);
      result->err *= fabs(b-a)+1.0;
      return GSL_SUCCESS;
    }
    else {
      /* Pick a0 such that b ~= 2a0 + x, then
       * recurse down in b from a0,a0 to determine
       * the values near the line b=2a+x. Then recurse
       * forward on a from a0.
       */
      int a0 = ceil(0.5*(b-x));
      double Ma0b;    /* M(a0,b)   */
      double Ma0bp1;  /* M(a0,b+1) */
      double Ma0p1b;  /* M(a0+1,b) */
      double Mnm1;
      double Mn;
      double Mnp1;
      int n;
      {
        double ex = exp(x);
    	double Ma0np1 = ex;
    	double Ma0n   = ex * (1.0 + x/(a0-1.0));
    	double Ma0nm1;
    	for(n=a0-1; n>b; n--) {
    	  Ma0nm1 = (-n*(1-n-x)*Ma0n - x*(n-a0)*Ma0np1)/(n*(n-1.0));
    	  Ma0np1 = Ma0n;
    	  Ma0n = Ma0nm1;
    	}
	Ma0bp1 = Ma0np1;
    	Ma0b   = Ma0n;
	Ma0p1b = (b*(a0+x)*Ma0b + x*(a0-b)*Ma0bp1)/(a0*b);
      }

      Mnm1 = Ma0b;
      Mn   = Ma0p1b;
      for(n=a0+1; n<a; n++) {
    	Mnp1 = ((b-n)*Mnm1 + (2*n-b+x)*Mn)/n;
    	Mnm1 = Mn;
    	Mn   = Mnp1;
      }
      result->val  = Mn;
      result->err  = (fabs(x) + 1.0) * GSL_DBL_EPSILON *  fabs(Mn);
      result->err *= fabs(b-a)+1.0;
      return GSL_SUCCESS;
    }
  }
}


/* Evaluate a <= 0, a integer, cases directly. (Polynomial; Horner)
 * When the terms are all positive, this
 * must work. We will assume this here.
 */
static
int
hyperg_1F1_a_negint_poly(const int a, const double b, const double x, gsl_sf_result * result)
{
  if(a == 0) {
    result->val = 1.0;
    result->err = 1.0;
    return GSL_SUCCESS;
  }
  else {
    int N = -a;
    double poly = 1.0;
    int k;
    for(k=N-1; k>=0; k--) {
      double t = (a+k)/(b+k) * (x/(k+1));
      double r = t + 1.0/poly;
      if(r > 0.9*GSL_DBL_MAX/poly) {
        OVERFLOW_ERROR(result);
      }
      else {
        poly *= r;  /* P_n = 1 + t_n P_{n-1} */
      }
    }
    result->val = poly;
    result->err = 2.0 * (sqrt(N) + 1.0) * GSL_DBL_EPSILON * fabs(poly);
    return GSL_SUCCESS;
  }
}


/* Evaluate negative integer a case by relation
 * to Laguerre polynomials. This is more general than
 * the direct polynomial evaluation, but is safe
 * for all values of x.
 *
 * 1F1(-n,b,x) = n!/(b)_n Laguerre[n,b-1,x]
 *	       = n B(b,n) Laguerre[n,b-1,x]
 *
 * assumes b is not a negative integer
 */
static
int
hyperg_1F1_a_negint_lag(const int a, const double b, const double x, gsl_sf_result * result)
{
  const int n = -a;

  gsl_sf_result lag;
  const int stat_l = gsl_sf_laguerre_n_e(n, b-1.0, x, &lag);
  if(b < 0.0) {
    gsl_sf_result lnfact;
    gsl_sf_result lng1;
    gsl_sf_result lng2;
    double s1, s2;
    const int stat_f  = gsl_sf_lnfact_e(n, &lnfact);
    const int stat_g1 = gsl_sf_lngamma_sgn_e(b + n, &lng1, &s1);
    const int stat_g2 = gsl_sf_lngamma_sgn_e(b, &lng2, &s2);
    const double lnpre_val = lnfact.val - (lng1.val - lng2.val);
    const double lnpre_err = lnfact.err + lng1.err + lng2.err
      + 2.0 * GSL_DBL_EPSILON * fabs(lnpre_val);
    const int stat_e = gsl_sf_exp_mult_err_e(lnpre_val, lnpre_err,
                                                s1*s2*lag.val, lag.err,
                                                result);
    return GSL_ERROR_SELECT_5(stat_e, stat_l, stat_g1, stat_g2, stat_f);
  }
  else {
    gsl_sf_result lnbeta;
    gsl_sf_lnbeta_e(b, n, &lnbeta);
    if(fabs(lnbeta.val) < 0.1) {
      /* As we have noted, when B(x,y) is near 1,
       * evaluating log(B(x,y)) is not accurate.
       * Instead we evaluate B(x,y) directly.
       */
      const double ln_term_val = log(1.25*n);
      const double ln_term_err = 2.0 * GSL_DBL_EPSILON * ln_term_val;
      gsl_sf_result beta;
      int stat_b = gsl_sf_beta_e(b, n, &beta);
      int stat_e = gsl_sf_exp_mult_err_e(ln_term_val, ln_term_err,
                                            lag.val, lag.err,
                                            result);
      result->val *= beta.val/1.25;
      result->err *= beta.val/1.25;
      return GSL_ERROR_SELECT_3(stat_e, stat_l, stat_b);
    }
    else {
      /* B(x,y) was not near 1, so it is safe to use
       * the logarithmic values.
       */
      const double ln_n = log(n);
      const double ln_term_val = lnbeta.val + ln_n;
      const double ln_term_err = lnbeta.err + 2.0 * GSL_DBL_EPSILON * fabs(ln_n);
      int stat_e = gsl_sf_exp_mult_err_e(ln_term_val, ln_term_err,
                                            lag.val, lag.err,
                                            result);
      return GSL_ERROR_SELECT_2(stat_e, stat_l);
    }
  }
}


/* Handle negative integer a case for x > 0 and
 * generic b.
 *
 * Combine [Abramowitz+Stegun, 13.6.9 + 13.6.27]
 * M(-n,b,x) = (-1)^n / (b)_n U(-n,b,x) = n! / (b)_n Laguerre^(b-1)_n(x)
 */
#if 0
static
int
hyperg_1F1_a_negint_U(const int a, const double b, const double x, gsl_sf_result * result)
{
  const int n = -a;
  const double sgn = ( GSL_IS_ODD(n) ? -1.0 : 1.0 );
  double sgpoch;
  gsl_sf_result lnpoch;
  gsl_sf_result U;
  const int stat_p = gsl_sf_lnpoch_sgn_e(b, n, &lnpoch, &sgpoch);
  const int stat_U = gsl_sf_hyperg_U_e(-n, b, x, &U);
  const int stat_e = gsl_sf_exp_mult_err_e(-lnpoch.val, lnpoch.err,
  					      sgn * sgpoch * U.val, U.err,
  					      result);
  return GSL_ERROR_SELECT_3(stat_e, stat_U, stat_p);
}
#endif


/* Assumes a <= -1,  b <= -1, and b <= a.
 */
static
int
hyperg_1F1_ab_negint(const int a, const int b, const double x, gsl_sf_result * result)
{
  if(x == 0.0) {
    result->val = 1.0;
    result->err = 0.0;
    return GSL_SUCCESS;
  }
  else if(x > 0.0) {
    return hyperg_1F1_a_negint_poly(a, b, x, result);
  }
  else {
    /* Apply a Kummer transformation to make x > 0 so
     * we can evaluate the polynomial safely. Of course,
     * this assumes b <= a, which must be true for
     * a<0 and b<0, since otherwise the thing is undefined.
     */
    gsl_sf_result K;
    int stat_K = hyperg_1F1_a_negint_poly(b-a, b, -x, &K);
    int stat_e = gsl_sf_exp_mult_err_e(x, 2.0 * GSL_DBL_EPSILON * fabs(x),
    					  K.val, K.err,
    					  result);
    return GSL_ERROR_SELECT_2(stat_e, stat_K);
  }
}


/* [Abramowitz+Stegun, 13.1.3]
 *
 * M(a,b,x) = Gamma(1+a-b)/Gamma(2-b) x^(1-b) *
 *            { Gamma(b)/Gamma(a) M(1+a-b,2-b,x) - (b-1) U(1+a-b,2-b,x) }
 *
 * b not an integer >= 2
 * a-b not a negative integer
 */
static
int
hyperg_1F1_U(const double a, const double b, const double x, gsl_sf_result * result)
{
  const double bp = 2.0 - b;
  const double ap = a - b + 1.0;

  gsl_sf_result lg_ap, lg_bp;
  double sg_ap;
  int stat_lg0 = gsl_sf_lngamma_sgn_e(ap, &lg_ap, &sg_ap);
  int stat_lg1 = gsl_sf_lngamma_e(bp, &lg_bp);
  int stat_lg2 = GSL_ERROR_SELECT_2(stat_lg0, stat_lg1);
  double t1 = (bp-1.0) * log(x);
  double lnpre_val = lg_ap.val - lg_bp.val + t1;
  double lnpre_err = lg_ap.err + lg_bp.err + 2.0 * GSL_DBL_EPSILON * fabs(t1);

  gsl_sf_result lg_2mbp, lg_1papmbp;
  double sg_2mbp, sg_1papmbp;
  int stat_lg3 = gsl_sf_lngamma_sgn_e(2.0-bp,    &lg_2mbp,    &sg_2mbp);
  int stat_lg4 = gsl_sf_lngamma_sgn_e(1.0+ap-bp, &lg_1papmbp, &sg_1papmbp);
  int stat_lg5 = GSL_ERROR_SELECT_2(stat_lg3, stat_lg4);
  double lnc1_val = lg_2mbp.val - lg_1papmbp.val;
  double lnc1_err = lg_2mbp.err + lg_1papmbp.err
                    + GSL_DBL_EPSILON * (fabs(lg_2mbp.val) + fabs(lg_1papmbp.val));

  gsl_sf_result M;
  gsl_sf_result_e10 U;
  int stat_F = gsl_sf_hyperg_1F1_e(ap, bp, x, &M);
  int stat_U = gsl_sf_hyperg_U_e10_e(ap, bp, x, &U);
  int stat_FU = GSL_ERROR_SELECT_2(stat_F, stat_U);

  gsl_sf_result_e10 term_M;
  int stat_e0 = gsl_sf_exp_mult_err_e10_e(lnc1_val, lnc1_err,
                                             sg_2mbp*sg_1papmbp*M.val, M.err,
                                             &term_M);

  const double ombp = 1.0 - bp;
  const double Uee_val = U.e10*M_LN10;
  const double Uee_err = 2.0 * GSL_DBL_EPSILON * fabs(Uee_val);
  const double Mee_val = term_M.e10*M_LN10;
  const double Mee_err = 2.0 * GSL_DBL_EPSILON * fabs(Mee_val);
  int stat_e1;

  /* Do a little dance with the exponential prefactors
   * to avoid overflows in intermediate results.
   */
  if(Uee_val > Mee_val) {
    const double factorM_val = exp(Mee_val-Uee_val);
    const double factorM_err = 2.0 * GSL_DBL_EPSILON * (fabs(Mee_val-Uee_val)+1.0) * factorM_val;
    const double inner_val = term_M.val*factorM_val - ombp*U.val;
    const double inner_err =
        term_M.err*factorM_val + fabs(ombp) * U.err
      + fabs(term_M.val) * factorM_err
      + GSL_DBL_EPSILON * (fabs(term_M.val*factorM_val) + fabs(ombp*U.val));
    stat_e1 = gsl_sf_exp_mult_err_e(lnpre_val+Uee_val, lnpre_err+Uee_err,
                                       sg_ap*inner_val, inner_err,
                                       result);
  }
  else {
    const double factorU_val = exp(Uee_val - Mee_val);
    const double factorU_err = 2.0 * GSL_DBL_EPSILON * (fabs(Mee_val-Uee_val)+1.0) * factorU_val;
    const double inner_val = term_M.val - ombp*factorU_val*U.val;
    const double inner_err =
        term_M.err + fabs(ombp*factorU_val*U.err)
      + fabs(ombp*factorU_err*U.val)
      + GSL_DBL_EPSILON * (fabs(term_M.val) + fabs(ombp*factorU_val*U.val));
    stat_e1 = gsl_sf_exp_mult_err_e(lnpre_val+Mee_val, lnpre_err+Mee_err,
                                       sg_ap*inner_val, inner_err,
                                       result);
  }

  return GSL_ERROR_SELECT_5(stat_e1, stat_e0, stat_FU, stat_lg5, stat_lg2);
}


/* Handle case of generic positive a, b.
 * Assumes b-a is not a negative integer.
 */
static
int
hyperg_1F1_ab_pos(const double a, const double b,
                  const double x,
                  gsl_sf_result * result)
{
  const double ax = fabs(x);

  if(   ( b < 10.0 && a < 10.0 && ax < 5.0 )
     || ( b > a*ax )
     || ( b > a && ax < 5.0 )
    ) {
    return gsl_sf_hyperg_1F1_series_e(a, b, x, result);
  }
  else if(   x < -100.0
          && GSL_MAX_DBL(fabs(a),1.0)*GSL_MAX_DBL(fabs(1.0+a-b),1.0) < 0.7*fabs(x)
    ) {
    /* Large negative x asymptotic.
     */
    return hyperg_1F1_asymp_negx(a, b, x, result);
  }
  else if(   x > 100.0
          && GSL_MAX_DBL(fabs(b-a),1.0)*GSL_MAX_DBL(fabs(1.0-a),1.0) < 0.7*fabs(x)
    ) {
    /* Large positive x asymptotic.
     */
    return hyperg_1F1_asymp_posx(a, b, x, result);
  }
  else if(fabs(b-a) <= 1.0) {
    /* Directly handle b near a.
     */
    return hyperg_1F1_beps_bgt0(a-b, b, x, result);  /* a = b + eps */
  }

  else if(b > a && b >= 2*a + x) {
    /* Use the Gautschi CF series, then
     * recurse backward to a near 0 for normalization.
     * This will work for either sign of x.
     */ 
    double rap;
    int stat_CF1 = hyperg_1F1_CF1_p_ser(a, b, x, &rap);
    double ra = 1.0 + x/a * rap;

    double Ma	= GSL_SQRT_DBL_MIN;
    double Map1 = ra * Ma;
    double Mnp1 = Map1;
    double Mn	= Ma;
    double Mnm1;
    gsl_sf_result Mn_true;
    int stat_Mt;
    double n;
    for(n=a; n>0.5; n -= 1.0) {
      Mnm1 = (n * Mnp1 - (2.0*n-b+x) * Mn) / (b-n);
      Mnp1 = Mn;
      Mn   = Mnm1;
    }

    stat_Mt = hyperg_1F1_small_a_bgt0(n, b, x, &Mn_true);

    result->val  = (Ma/Mn) * Mn_true.val;
    result->err  = fabs(Ma/Mn) * Mn_true.err;
    result->err += 2.0 * GSL_DBL_EPSILON * (fabs(a)+1.0) * fabs(result->val);
    return GSL_ERROR_SELECT_2(stat_Mt, stat_CF1);
  }
  else if(b > a && b < 2*a + x && b > x) {
    /* Use the Gautschi series representation of
     * the continued fraction. Then recurse forward
     * to near the a=b line for normalization. This will
     * work for either sign of x, although we do need
     * to check for b > x, which is relevant when x is positive.
     */
    gsl_sf_result Mn_true;
    int stat_Mt;
    double rap;
    int stat_CF1 = hyperg_1F1_CF1_p_ser(a, b, x, &rap);
    double ra = 1.0 + x/a * rap;
    double Ma	= GSL_SQRT_DBL_MIN;
    double Mnm1 = Ma;
    double Mn	= ra * Mnm1;
    double Mnp1;
    double n;
    for(n=a+1.0; n<b-0.5; n += 1.0) {
      Mnp1 = ((b-n)*Mnm1 + (2*n-b+x)*Mn)/n;
      Mnm1 = Mn;
      Mn   = Mnp1;
    }
    stat_Mt = hyperg_1F1_beps_bgt0(n-b, b, x, &Mn_true);
    result->val  = Ma/Mn * Mn_true.val;
    result->err  = fabs(Ma/Mn) * Mn_true.err;
    result->err += 2.0 * GSL_DBL_EPSILON * (fabs(b-a)+1.0) * fabs(result->val);
    return GSL_ERROR_SELECT_2(stat_Mt, stat_CF1);
  }
  else if(x >= 0.0) {

    if(b < a) {
      /* Forward recursion on a from a=b+eps-1,b+eps.
       */
      double N   = floor(a-b);
      double eps = a - b - N;
      gsl_sf_result r_M0;
      gsl_sf_result r_M1;
      int stat_0 = hyperg_1F1_beps_bgt0(eps-1.0, b, x, &r_M0);
      int stat_1 = hyperg_1F1_beps_bgt0(eps,     b, x, &r_M1);
      double M0 = r_M0.val;
      double M1 = r_M1.val;

      double Mam1 = M0;
      double Ma   = M1;
      double Map1;
      double ap;
      double start_pair = fabs(M0) + fabs(M1);
      double minim_pair = GSL_DBL_MAX;
      double pair_ratio;
      double rat_0 = fabs(r_M0.err/r_M0.val);
      double rat_1 = fabs(r_M1.err/r_M1.val);
      for(ap=b+eps; ap<a-0.1; ap += 1.0) {
        Map1 = ((b-ap)*Mam1 + (2.0*ap-b+x)*Ma)/ap;
        Mam1 = Ma;
        Ma   = Map1;
        minim_pair = GSL_MIN_DBL(fabs(Mam1) + fabs(Ma), minim_pair);
      }
      pair_ratio = start_pair/minim_pair;
      result->val  = Ma;
      result->err  = 2.0 * (rat_0 + rat_1 + GSL_DBL_EPSILON) * (fabs(b-a)+1.0) * fabs(Ma);
      result->err += 2.0 * (rat_0 + rat_1) * pair_ratio*pair_ratio * fabs(Ma);
      result->err += 2.0 * GSL_DBL_EPSILON * fabs(Ma);
      return GSL_ERROR_SELECT_2(stat_0, stat_1);
    }
    else {
      /* b > a
       * b < 2a + x 
       * b <= x
       *
       * Recurse forward on a from a=eps,eps+1.
       */
      double eps = a - floor(a);
      gsl_sf_result r_Mnm1;
      gsl_sf_result r_Mn;
      int stat_0 = hyperg_1F1_small_a_bgt0(eps,     b, x, &r_Mnm1);
      int stat_1 = hyperg_1F1_small_a_bgt0(eps+1.0, b, x, &r_Mn);
      double Mnm1 = r_Mnm1.val;
      double Mn   = r_Mn.val;
      double Mnp1;

      double n;
      double start_pair = fabs(Mn) + fabs(Mnm1);
      double minim_pair = GSL_DBL_MAX;