hyperg_1f1.c

开放gsl矩阵运算

 */
static
int
hyperg_1F1_CF1_p_ser(const double a, const double b, const double x, double * result)
{
  if(a == 0.0) {
    *result = 0.0;
    return GSL_SUCCESS;
  }
  else {
    const int maxiter = 5000;
    double sum  = 1.0;
    double pk   = 1.0;
    double rhok = 0.0;
    int k;
    for(k=1; k<maxiter; k++) {
      double ak = (a + k)*x/((b-x+k-1.0)*(b-x+k));
      rhok = -ak*(1.0 + rhok)/(1.0 + ak*(1.0+rhok));
      pk  *= rhok;
      sum += pk;
      if(fabs(pk/sum) < 2.0*GSL_DBL_EPSILON) break;
    }
    *result = a/(b-x) * sum;
    if(k == maxiter)
      GSL_ERROR ("error", GSL_EMAXITER);
    else
      return GSL_SUCCESS;
  }
}


/* 1F1(a+1,b,x)/1F1(a,b,x)
 *
 * I think this suffers from typical "anomalous convergence".
 * I could not find a region where it was truly useful.
 */
#if 0
static
int
hyperg_1F1_CF1(const double a, const double b, const double x, double * result)
{
  const double RECUR_BIG = GSL_SQRT_DBL_MAX;
  const int maxiter = 5000;
  int n = 1;
  double Anm2 = 1.0;
  double Bnm2 = 0.0;
  double Anm1 = 0.0;
  double Bnm1 = 1.0;
  double a1 = b - a - 1.0;
  double b1 = b - x - 2.0*(a+1.0);
  double An = b1*Anm1 + a1*Anm2;
  double Bn = b1*Bnm1 + a1*Bnm2;
  double an, bn;
  double fn = An/Bn;

  while(n < maxiter) {
    double old_fn;
    double del;
    n++;
    Anm2 = Anm1;
    Bnm2 = Bnm1;
    Anm1 = An;
    Bnm1 = Bn;
    an = (a + n - 1.0) * (b - a - n);
    bn = b - x - 2.0*(a+n);
    An = bn*Anm1 + an*Anm2;
    Bn = bn*Bnm1 + an*Bnm2;

    if(fabs(An) > RECUR_BIG || fabs(Bn) > RECUR_BIG) {
      An /= RECUR_BIG;
      Bn /= RECUR_BIG;
      Anm1 /= RECUR_BIG;
      Bnm1 /= RECUR_BIG;
      Anm2 /= RECUR_BIG;
      Bnm2 /= RECUR_BIG;
    }

    old_fn = fn;
    fn = An/Bn;
    del = old_fn/fn;
    
    if(fabs(del - 1.0) < 10.0*GSL_DBL_EPSILON) break;
  }

  *result = fn;
  if(n == maxiter)
    GSL_ERROR ("error", GSL_EMAXITER);
  else
    return GSL_SUCCESS;
}
#endif /* 0 */


/* 1F1(a,b+1,x)/1F1(a,b,x)
 *
 * This seemed to suffer from "anomalous convergence".
 * However, I have no theory for this recurrence.
 */
#if 0
static
int
hyperg_1F1_CF1_b(const double a, const double b, const double x, double * result)
{
  const double RECUR_BIG = GSL_SQRT_DBL_MAX;
  const int maxiter = 5000;
  int n = 1;
  double Anm2 = 1.0;
  double Bnm2 = 0.0;
  double Anm1 = 0.0;
  double Bnm1 = 1.0;
  double a1 = b + 1.0;
  double b1 = (b + 1.0) * (b - x);
  double An = b1*Anm1 + a1*Anm2;
  double Bn = b1*Bnm1 + a1*Bnm2;
  double an, bn;
  double fn = An/Bn;

  while(n < maxiter) {
    double old_fn;
    double del;
    n++;
    Anm2 = Anm1;
    Bnm2 = Bnm1;
    Anm1 = An;
    Bnm1 = Bn;
    an = (b + n) * (b + n - 1.0 - a) * x;
    bn = (b + n) * (b + n - 1.0 - x);
    An = bn*Anm1 + an*Anm2;
    Bn = bn*Bnm1 + an*Bnm2;

    if(fabs(An) > RECUR_BIG || fabs(Bn) > RECUR_BIG) {
      An /= RECUR_BIG;
      Bn /= RECUR_BIG;
      Anm1 /= RECUR_BIG;
      Bnm1 /= RECUR_BIG;
      Anm2 /= RECUR_BIG;
      Bnm2 /= RECUR_BIG;
    }

    old_fn = fn;
    fn = An/Bn;
    del = old_fn/fn;
    
    if(fabs(del - 1.0) < 10.0*GSL_DBL_EPSILON) break;
  }

  *result = fn;
  if(n == maxiter)
    GSL_ERROR ("error", GSL_EMAXITER);
  else
    return GSL_SUCCESS;
}
#endif /* 0 */


/* 1F1(a,b,x)
 * |a| <= 1, b > 0
 */
static
int
hyperg_1F1_small_a_bgt0(const double a, const double b, const double x, gsl_sf_result * result)
{
  const double bma = b-a;
  const double oma = 1.0-a;
  const double ap1mb = 1.0+a-b;
  const double abs_bma = fabs(bma);
  const double abs_oma = fabs(oma);
  const double abs_ap1mb = fabs(ap1mb);

  const double ax = fabs(x);

  if(a == 0.0) {
    result->val = 1.0;
    result->err = 0.0;
    return GSL_SUCCESS;
  }
  else if(a == 1.0) {
    return hyperg_1F1_1(b, x, result);
  }
  else if(a == -1.0) {
    result->val  = 1.0 + a/b * x;
    result->err  = GSL_DBL_EPSILON * (1.0 + fabs(a/b * x));
    result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
    return GSL_SUCCESS;
  }
  else if(b >= 1.4*ax) {
    return gsl_sf_hyperg_1F1_series_e(a, b, x, result);
  }
  else if(x > 0.0) {
    if(x > 100.0 && abs_bma*abs_oma < 0.5*x) {
      return hyperg_1F1_asymp_posx(a, b, x, result);
    }
    else if(b < 5.0e+06) {
      /* Recurse backward on b from
       * a suitably high point.
       */
      const double b_del = ceil(1.4*x-b) + 1.0;
      double bp = b + b_del;
      gsl_sf_result r_Mbp1;
      gsl_sf_result r_Mb;
      double Mbp1;
      double Mb;
      double Mbm1;
      int stat_0 = gsl_sf_hyperg_1F1_series_e(a, bp+1.0, x, &r_Mbp1);
      int stat_1 = gsl_sf_hyperg_1F1_series_e(a, bp,     x, &r_Mb);
      const double err_rat = fabs(r_Mbp1.err/r_Mbp1.val) + fabs(r_Mb.err/r_Mb.val);
      Mbp1 = r_Mbp1.val;
      Mb   = r_Mb.val;
      while(bp > b+0.1) {
        /* Do backward recursion. */
        Mbm1 = ((x+bp-1.0)*Mb - x*(bp-a)/bp*Mbp1)/(bp-1.0);
        bp -= 1.0;
	Mbp1 = Mb;
	Mb   = Mbm1;
      }
      result->val  = Mb;
      result->err  = err_rat * (fabs(b_del)+1.0) * fabs(Mb);
      result->err += 2.0 * GSL_DBL_EPSILON * fabs(Mb);
      return GSL_ERROR_SELECT_2(stat_0, stat_1);
    }
    else {
      return hyperg_1F1_large2bm4a(a, b, x, result);
    }
  }
  else {
    /* x < 0 and b not large compared to |x|
     */
    if(ax < 10.0 && b < 10.0) {
      return gsl_sf_hyperg_1F1_series_e(a, b, x, result);
    }
    else if(ax >= 100.0 && GSL_MAX(abs_ap1mb,1.0) < 0.99*ax) {
      return hyperg_1F1_asymp_negx(a, b, x, result);
    }
    else {
      return hyperg_1F1_luke(a, b, x, result);
    }
  }
}


/* 1F1(b+eps,b,x)
 * |eps|<=1, b > 0
 */
static
int
hyperg_1F1_beps_bgt0(const double eps, const double b, const double x, gsl_sf_result * result)
{
  if(b > fabs(x) && fabs(eps) < GSL_SQRT_DBL_EPSILON) {
    /* If b-a is very small and x/b is not too large we can
     * use this explicit approximation.
     *
     * 1F1(b+eps,b,x) = exp(ax/b) (1 - eps x^2 (v2 + v3 x + ...) + ...)
     *
     *   v2 = a/(2b^2(b+1))
     *   v3 = a(b-2a)/(3b^3(b+1)(b+2))
     *   ...
     *
     * See [Luke, Mathematical Functions and Their Approximations, p.292]
     *
     * This cannot be used for b near a negative integer or zero.
     * Also, if x/b is large the deviation from exp(x) behaviour grows.
     */
    double a = b + eps;
    gsl_sf_result exab;
    int stat_e = gsl_sf_exp_e(a*x/b, &exab);
    double v2 = a/(2.0*b*b*(b+1.0));
    double v3 = a*(b-2.0*a)/(3.0*b*b*b*(b+1.0)*(b+2.0));
    double v  = v2 + v3 * x;
    double f  = (1.0 - eps*x*x*v);
    result->val  = exab.val * f;
    result->err  = exab.err * fabs(f);
    result->err += fabs(exab.val) * GSL_DBL_EPSILON * (1.0 + fabs(eps*x*x*v));
    result->err += 4.0 * GSL_DBL_EPSILON * fabs(result->val);
    return stat_e;
  }
  else {
    /* Otherwise use a Kummer transformation to reduce
     * it to the small a case.
     */
    gsl_sf_result Kummer_1F1;
    int stat_K = hyperg_1F1_small_a_bgt0(-eps, b, -x, &Kummer_1F1);
    if(Kummer_1F1.val != 0.0) {
      int stat_e = gsl_sf_exp_mult_err_e(x, 2.0*GSL_DBL_EPSILON*fabs(x),
                                            Kummer_1F1.val, Kummer_1F1.err,
                                            result);
      return GSL_ERROR_SELECT_2(stat_e, stat_K);
    }
    else {
      result->val = 0.0;
      result->err = 0.0;
      return stat_K;
    }
  }
}


/* 1F1(a,2a,x) = Gamma(a + 1/2) E(x) (|x|/4)^(-a+1/2) scaled_I(a-1/2,|x|/2)
 *
 * E(x) = exp(x) x > 0
 *      = 1      x < 0
 *
 * a >= 1/2
 */
static
int
hyperg_1F1_beq2a_pos(const double a, const double x, gsl_sf_result * result)
{
  if(x == 0.0) {
    result->val = 1.0;
    result->err = 0.0;
    return GSL_SUCCESS;
  }
  else {
    gsl_sf_result I;
    int stat_I = gsl_sf_bessel_Inu_scaled_e(a-0.5, 0.5*fabs(x), &I);
    gsl_sf_result lg;
    int stat_g = gsl_sf_lngamma_e(a + 0.5, &lg);
    double ln_term   = (0.5-a)*log(0.25*fabs(x));
    double lnpre_val = lg.val + GSL_MAX_DBL(x,0.0) + ln_term;
    double lnpre_err = lg.err + GSL_DBL_EPSILON * (fabs(ln_term) + fabs(x));
    int stat_e = gsl_sf_exp_mult_err_e(lnpre_val, lnpre_err,
    					  I.val, I.err,
    					  result);
    return GSL_ERROR_SELECT_3(stat_e, stat_g, stat_I);
  }
}


/* Determine middle parts of diagonal recursion along b=2a
 * from two endpoints, i.e.
 *
 * given:  M(a,b)      and  M(a+1,b+2)
 * get:    M(a+1,b+1)  and  M(a,b+1)
 */
#if 0
inline
static
int
hyperg_1F1_diag_step(const double a, const double b, const double x,
                     const double Mab, const double Map1bp2,
                     double * Map1bp1, double * Mabp1)
{
  if(a == b) {
    *Map1bp1 = Mab;
    *Mabp1   = Mab - x/(b+1.0) * Map1bp2;
  }
  else {
    *Map1bp1 = Mab - x * (a-b)/(b*(b+1.0)) * Map1bp2;
    *Mabp1   = (a * *Map1bp1 - b * Mab)/(a-b);
  }
  return GSL_SUCCESS;
}
#endif /* 0 */


/* Determine endpoint of diagonal recursion.
 *
 * given:  M(a,b)    and  M(a+1,b+2)
 * get:    M(a+1,b)  and  M(a+1,b+1)
 */
#if 0
inline
static
int
hyperg_1F1_diag_end_step(const double a, const double b, const double x,
                         const double Mab, const double Map1bp2,
                         double * Map1b, double * Map1bp1)
{
  *Map1bp1 = Mab - x * (a-b)/(b*(b+1.0)) * Map1bp2;
  *Map1b   = Mab + x/b * *Map1bp1;
  return GSL_SUCCESS;
}
#endif /* 0 */


/* Handle the case of a and b both positive integers.
 * Assumes a > 0 and b > 0.
 */
static
int
hyperg_1F1_ab_posint(const int a, const int b, const double x, gsl_sf_result * result)
{
  double ax = fabs(x);

  if(a == b) {
    return gsl_sf_exp_e(x, result);             /* 1F1(a,a,x) */
  }
  else if(a == 1) {
    return gsl_sf_exprel_n_e(b-1, x, result);   /* 1F1(1,b,x) */
  }
  else if(b == a + 1) {
    gsl_sf_result K;
    int stat_K = gsl_sf_exprel_n_e(a, -x, &K);  /* 1F1(1,1+a,-x) */
    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);
  }
  else if(a == b + 1) {
    gsl_sf_result ex;
    int stat_e = gsl_sf_exp_e(x, &ex);
    result->val  = ex.val * (1.0 + x/b);
    result->err  = ex.err * (1.0 + x/b);
    result->err += ex.val * GSL_DBL_EPSILON * (1.0 + fabs(x/b));
    result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
    return stat_e;
  }
  else if(a == b + 2) {
    gsl_sf_result ex;
    int stat_e = gsl_sf_exp_e(x, &ex);
    double poly  = (1.0 + x/b*(2.0 + x/(b+1.0)));
    result->val  = ex.val * poly;
    result->err  = ex.err * fabs(poly);
    result->err += ex.val * GSL_DBL_EPSILON * (1.0 + fabs(x/b) * (2.0 + fabs(x/(b+1.0))));
    result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
    return stat_e;
  }
  else if(b == 2*a) {
    return hyperg_1F1_beq2a_pos(a, x, result);  /* 1F1(a,2a,x) */
  }
  else if(   ( b < 10 && a < 10 && ax < 5.0 )
          || ( b > a*ax )
	  || ( b > a && ax < 5.0 )
    ) {
    return gsl_sf_hyperg_1F1_series_e(a, b, x, result);
  }
  else if(b > a && b >= 2*a + x) {
    /* Use the Gautschi CF series, then
     * recurse backward to a=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;
    int n;
    for(n=a; n>0; n--) {
      Mnm1 = (n * Mnp1 - (2*n-b+x) * Mn) / (b-n);
      Mnp1 = Mn;
      Mn   = Mnm1;
    }
    result->val = Ma/Mn;
    result->err = 2.0 * GSL_DBL_EPSILON * (fabs(a) + 1.0) * fabs(Ma/Mn);
    return 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 the a=b line for normalization. This will
     * work for either sign of x, although we do need
     * to check for b > x, for when x is positive.
     */
    double rap;
    int stat_CF1 = hyperg_1F1_CF1_p_ser(a, b, x, &rap);
    double ra = 1.0 + x/a * rap;
    gsl_sf_result ex;
    int stat_ex;

    double Ma   = GSL_SQRT_DBL_MIN;
    double Map1 = ra * Ma;
    double Mnm1 = Ma;
    double Mn	= Map1;
    double Mnp1;
    int n;
    for(n=a+1; n<b; n++) {
      Mnp1 = ((b-n)*Mnm1 + (2*n-b+x)*Mn)/n;
      Mnm1 = Mn;
      Mn   = Mnp1;
    }

    stat_ex = gsl_sf_exp_e(x, &ex);  /* 1F1(b,b,x) */
    result->val  = ex.val * Ma/Mn;
    result->err  = ex.err * fabs(Ma/Mn);
    result->err += 4.0 * GSL_DBL_EPSILON * (fabs(b-a)+1.0) * fabs(result->val);
    return GSL_ERROR_SELECT_2(stat_ex, stat_CF1);
  }
  else if(x >= 0.0) {

    if(b < a) {
      /* The point b,b is below the b=2a+x line.
       * Forward recursion on a from b,b+1 is possible.
       * Note that a > b + 1 as well, since we already tried a = b + 1.
       */
      if(x + log(fabs(x/b)) < GSL_LOG_DBL_MAX-2.0) {
        double ex = exp(x);
        int n;
        double Mnm1 = ex;		  /* 1F1(b,b,x)   */
    	double Mn   = ex * (1.0 + x/b);   /* 1F1(b+1,b,x) */
    	double Mnp1;
        for(n=b+1; n<a; n++) {
          Mnp1 = ((b-n)*Mnm1 + (2*n-b+x)*Mn)/n;
          Mnm1 = Mn;
          Mn   = Mnp1;
        }