hyperg_u.c

开放gsl矩阵运算

        const double xi  = istrt + i;
        const double xi1 = istrt + i - 1;
	const double tmp = (a-1.0)*(xn+2.0*xi-1.0) + xi*(xi-beps);
	const double b0_multiplier = (a+xi1-beps)*x/((xn+xi1)*(xi-beps));
	const double c0_multiplier_1 = (a+xi1)*x/((b+xi1)*xi);
	const double c0_multiplier_2 = tmp / (xi*(b+xi1)*(a+xi1-beps));
        b0_val *= b0_multiplier;
	b0_err += fabs(b0_multiplier) * b0_err + fabs(b0_val) * 8.0 * 2.0 * GSL_DBL_EPSILON;
        c0_val  = c0_multiplier_1 * c0_val - c0_multiplier_2 * b0_val;
	c0_err  =  fabs(c0_multiplier_1) * c0_err
	         + fabs(c0_multiplier_2) * b0_err
		 + fabs(c0_val) * 8.0 * 2.0 * GSL_DBL_EPSILON
		 + fabs(b0_val * c0_multiplier_2) * 16.0 * 2.0 * GSL_DBL_EPSILON;
        t_val  = c0_val + xeps1_val*b0_val;
	t_err  = c0_err + fabs(xeps1_val)*b0_err;
	t_err += fabs(b0_val*lnx) * dexprl.err;
	t_err += fabs(b0_val)*xeps1_err;
	dchu_val += t_val;
	dchu_err += t_err;
        if(fabs(t_val) < EPS*fabs(dchu_val)) break;
      }

      result->val  = dchu_val;
      result->err  = 2.0 * dchu_err;
      result->err += 2.0 * fabs(t_val);
      result->err += 4.0 * GSL_DBL_EPSILON * (i+2.0) * fabs(dchu_val);
      result->err *= 2.0; /* FIXME: fudge factor */

      if(i >= 2000) {
        GSL_ERROR ("error", GSL_EMAXITER);
      }
      else {
        return stat_all;
      }
    }
    else {
      /*
       C  X**(-BEPS) IS VERY DIFFERENT FROM 1.0, SO THE
       C  STRAIGHTFORWARD FORMULATION IS STABLE.
       */
      int i;
      double dchu_val;
      double dchu_err;
      double t_val;
      double t_err;
      gsl_sf_result dgamrbxi;
      int stat_dgamrbxi = gsl_sf_gammainv_e(b+xi, &dgamrbxi);
      double a0_val = factor_val * pochai.val * dgamrbxi.val * gamri1.val / beps;
      double a0_err =  fabs(factor_val * pochai.val * dgamrbxi.val / beps) * gamri1.err
                     + fabs(factor_val * pochai.val * gamri1.val / beps) * dgamrbxi.err
		     + fabs(factor_val * dgamrbxi.val * gamri1.val / beps) * pochai.err
		     + fabs(pochai.val * dgamrbxi.val * gamri1.val / beps) * factor_err
                     + 2.0 * GSL_DBL_EPSILON * fabs(a0_val);
      stat_all = GSL_ERROR_SELECT_2(stat_all, stat_dgamrbxi);

      b0_val = xeps * b0_val / beps;
      b0_err = fabs(xeps / beps) * b0_err + 4.0 * GSL_DBL_EPSILON * fabs(b0_val);
      dchu_val = sum.val + a0_val - b0_val;
      dchu_err = sum.err + a0_err + b0_err
        + 2.0 * GSL_DBL_EPSILON * (fabs(sum.val) + fabs(a0_val) + fabs(b0_val));

      for(i=1; i<2000; i++) {
        double xi = istrt + i;
        double xi1 = istrt + i - 1;
	double a0_multiplier = (a+xi1)*x/((b+xi1)*xi);
	double b0_multiplier = (a+xi1-beps)*x/((aintb+xi1)*(xi-beps));
        a0_val *= a0_multiplier;
	a0_err += fabs(a0_multiplier) * a0_err;
        b0_val *= b0_multiplier;
	b0_err += fabs(b0_multiplier) * b0_err;
        t_val = a0_val - b0_val;
	t_err = a0_err + b0_err;
        dchu_val += t_val;
	dchu_err += t_err;
        if(fabs(t_val) < EPS*fabs(dchu_val)) break;
      }

      result->val  = dchu_val;
      result->err  = 2.0 * dchu_err;
      result->err += 2.0 * fabs(t_val);
      result->err += 4.0 * GSL_DBL_EPSILON * (i+2.0) * fabs(dchu_val);
      result->err *= 2.0; /* FIXME: fudge factor */

      if(i >= 2000) {
        GSL_ERROR ("error", GSL_EMAXITER);
      }
      else {
        return stat_all;
      }
    }
  }
}


/* Assumes b > 0 and x > 0.
 */
static
int
hyperg_U_small_ab(const double a, const double b, const double x, gsl_sf_result * result)
{
  if(a == -1.0) {
    /* U(-1,c+1,x) = Laguerre[c,0,x] = -b + x
     */
    result->val  = -b + x;
    result->err  = 2.0 * GSL_DBL_EPSILON * (fabs(b) + fabs(x));
    result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
    return GSL_SUCCESS;
  }
  else if(a == 0.0) {
    /* U(0,c+1,x) = Laguerre[c,0,x] = 1
     */
    result->val = 1.0;
    result->err = 0.0;
    return GSL_SUCCESS;
  }
  else if(ASYMP_EVAL_OK(a,b,x)) {
    double p = pow(x, -a);
    gsl_sf_result asymp;
    int stat_asymp = hyperg_zaU_asymp(a, b, x, &asymp);
    result->val  = asymp.val * p;
    result->err  = asymp.err * p;
    result->err += fabs(asymp.val) * GSL_DBL_EPSILON * fabs(a) * p;
    result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
    return stat_asymp;
  }
  else {
    return hyperg_U_series(a, b, x, result);
  }
}


/* Assumes b > 0 and x > 0.
 */
static
int
hyperg_U_small_a_bgt0(const double a, const double b, const double x,
                      gsl_sf_result * result,
		      double * ln_multiplier
		      )
{
  if(a == 0.0) {
    result->val = 1.0;
    result->err = 1.0;
    *ln_multiplier = 0.0;
    return GSL_SUCCESS;
  }
  else if(   (b > 5000.0 && x < 0.90 * fabs(b))
          || (b >  500.0 && x < 0.50 * fabs(b))
    ) {
    int stat = gsl_sf_hyperg_U_large_b_e(a, b, x, result, ln_multiplier);
    if(stat == GSL_EOVRFLW)
      return GSL_SUCCESS;
    else
      return stat;
  }
  else if(b > 15.0) {
    /* Recurse up from b near 1.
     */
    double eps = b - floor(b);
    double b0  = 1.0 + eps;
    gsl_sf_result r_Ubm1;
    gsl_sf_result r_Ub;
    int stat_0 = hyperg_U_small_ab(a, b0,     x, &r_Ubm1);
    int stat_1 = hyperg_U_small_ab(a, b0+1.0, x, &r_Ub);
    double Ubm1 = r_Ubm1.val;
    double Ub   = r_Ub.val;
    double Ubp1;
    double bp;

    for(bp = b0+1.0; bp<b-0.1; bp += 1.0) {
      Ubp1 = ((1.0+a-bp)*Ubm1 + (bp+x-1.0)*Ub)/x;
      Ubm1 = Ub;
      Ub   = Ubp1;
    }
    result->val  = Ub;
    result->err  = (fabs(r_Ubm1.err/r_Ubm1.val) + fabs(r_Ub.err/r_Ub.val)) * fabs(Ub);
    result->err += 2.0 * GSL_DBL_EPSILON * (fabs(b-b0)+1.0) * fabs(Ub);
    *ln_multiplier = 0.0;
    return GSL_ERROR_SELECT_2(stat_0, stat_1);
  }
  else {
    *ln_multiplier = 0.0;
    return hyperg_U_small_ab(a, b, x, result);
  }
}


/* We use this to keep track of large
 * dynamic ranges in the recursions.
 * This can be important because sometimes
 * we want to calculate a very large and
 * a very small number and the answer is
 * the product, of order 1. This happens,
 * for instance, when we apply a Kummer
 * transform to make b positive and
 * both x and b are large.
 */
#define RESCALE_2(u0,u1,factor,count)      \
do {                                       \
  double au0 = fabs(u0);                   \
  if(au0 > factor) {                       \
    u0 /= factor;                          \
    u1 /= factor;                          \
    count++;                               \
  }                                        \
  else if(au0 < 1.0/factor) {              \
    u0 *= factor;                          \
    u1 *= factor;                          \
    count--;                               \
  }                                        \
} while (0)


/* Specialization to b >= 1, for integer parameters.
 * Assumes x > 0.
 */
static
int
hyperg_U_int_bge1(const int a, const int b, const double x,
                  gsl_sf_result_e10 * result)
{
  if(a == 0) {
    result->val = 1.0;
    result->err = 0.0;
    result->e10 = 0;
    return GSL_SUCCESS;
  }
  else if(a == -1) {
    result->val  = -b + x;
    result->err  = 2.0 * GSL_DBL_EPSILON * (fabs(b) + fabs(x));
    result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
    result->e10  = 0;
    return GSL_SUCCESS;
  }
  else if(b == a + 1) {
    /* U(a,a+1,x) = x^(-a)
     */
    return gsl_sf_exp_e10_e(-a*log(x), result);
  }
  else if(ASYMP_EVAL_OK(a,b,x)) {
    const double ln_pre_val = -a*log(x);
    const double ln_pre_err = 2.0 * GSL_DBL_EPSILON * fabs(ln_pre_val);
    gsl_sf_result asymp;
    int stat_asymp = hyperg_zaU_asymp(a, b, x, &asymp);
    int stat_e = gsl_sf_exp_mult_err_e10_e(ln_pre_val, ln_pre_err,
                                              asymp.val, asymp.err,
                                              result);
    return GSL_ERROR_SELECT_2(stat_e, stat_asymp);
  }
  else if(SERIES_EVAL_OK(a,b,x)) {
    gsl_sf_result ser;
    const int stat_ser = hyperg_U_series(a, b, x, &ser);
    result->val = ser.val;
    result->err = ser.err;
    result->e10 = 0;
    return stat_ser;
  }
  else if(a < 0) {
    /* Recurse backward from a = -1,0.
     */
    int scale_count = 0;
    const double scale_factor = GSL_SQRT_DBL_MAX;
    gsl_sf_result lnm;
    gsl_sf_result y;
    double lnscale;
    double Uap1 = 1.0;     /* U(0,b,x)  */
    double Ua   = -b + x;  /* U(-1,b,x) */
    double Uam1;
    int ap;

    for(ap=-1; ap>a; ap--) {
      Uam1 = ap*(b-ap-1.0)*Uap1 + (x+2.0*ap-b)*Ua;
      Uap1 = Ua;
      Ua   = Uam1;
      RESCALE_2(Ua,Uap1,scale_factor,scale_count);
    }

    lnscale = log(scale_factor);
    lnm.val = scale_count*lnscale;
    lnm.err = 2.0 * GSL_DBL_EPSILON * fabs(lnm.val);
    y.val = Ua;
    y.err = 4.0 * GSL_DBL_EPSILON * (fabs(a)+1.0) * fabs(Ua);
    return gsl_sf_exp_mult_err_e10_e(lnm.val, lnm.err, y.val, y.err, result);
  }
  else if(b >= 2.0*a + x) {
    /* Recurse forward from a = 0,1.
     */
    int scale_count = 0;
    const double scale_factor = GSL_SQRT_DBL_MAX;
    gsl_sf_result r_Ua;
    gsl_sf_result lnm;
    gsl_sf_result y;
    double lnscale;
    double lm;
    int stat_1 = hyperg_U_small_a_bgt0(1.0, b, x, &r_Ua, &lm);  /* U(1,b,x) */
    int stat_e;
    double Uam1 = 1.0;  /* U(0,b,x) */
    double Ua   = r_Ua.val;
    double Uap1;
    int ap;

    Uam1 *= exp(-lm);

    for(ap=1; ap<a; ap++) {
      Uap1 = -(Uam1 + (b-2.0*ap-x)*Ua)/(ap*(1.0+ap-b));
      Uam1 = Ua;
      Ua   = Uap1;
      RESCALE_2(Ua,Uam1,scale_factor,scale_count);
    }

    lnscale = log(scale_factor);
    lnm.val = lm + scale_count * lnscale;
    lnm.err = 2.0 * GSL_DBL_EPSILON * (fabs(lm) + fabs(scale_count*lnscale));
    y.val  = Ua;
    y.err  = fabs(r_Ua.err/r_Ua.val) * fabs(Ua);
    y.err += 2.0 * GSL_DBL_EPSILON * (fabs(a) + 1.0) * fabs(Ua);
    stat_e = gsl_sf_exp_mult_err_e10_e(lnm.val, lnm.err, y.val, y.err, result);
    return GSL_ERROR_SELECT_2(stat_e, stat_1);
  }
  else {
    if(b <= x) {
      /* Recurse backward either to the b=a+1 line
       * or to a=0, whichever we hit.
       */
      const double scale_factor = GSL_SQRT_DBL_MAX;
      int scale_count = 0;
      int stat_CF1;
      double ru;
      int CF1_count;
      int a_target;
      double lnU_target;
      double Ua;
      double Uap1;
      double Uam1;
      int ap;

      if(b < a + 1) {
        a_target = b-1;
	lnU_target = -a_target*log(x);
      }
      else {
        a_target = 0;
	lnU_target = 0.0;
      }

      stat_CF1 = hyperg_U_CF1(a, b, 0, x, &ru, &CF1_count);

      Ua   = 1.0;
      Uap1 = ru/a * Ua;
      for(ap=a; ap>a_target; ap--) {
        Uam1 = -((b-2.0*ap-x)*Ua + ap*(1.0+ap-b)*Uap1);
        Uap1 = Ua;
        Ua   = Uam1;
        RESCALE_2(Ua,Uap1,scale_factor,scale_count);
      }

      if(Ua == 0.0) {
        result->val = 0.0;
        result->err = 0.0;
	result->e10 = 0;
        GSL_ERROR ("error", GSL_EZERODIV);
      }
      else {
        double lnscl = -scale_count*log(scale_factor);
        double lnpre_val = lnU_target + lnscl;
        double lnpre_err = 2.0 * GSL_DBL_EPSILON * (fabs(lnU_target) + fabs(lnscl));
        double oUa_err   = 2.0 * (fabs(a_target-a) + CF1_count + 1.0) * GSL_DBL_EPSILON * fabs(1.0/Ua);
        int stat_e = gsl_sf_exp_mult_err_e10_e(lnpre_val, lnpre_err,
        				          1.0/Ua, oUa_err,
        				          result);
        return GSL_ERROR_SELECT_2(stat_e, stat_CF1);
      }
    }
    else {
      /* Recurse backward to near the b=2a+x line, then
       * determine normalization by either direct evaluation
       * or by a forward recursion. The direct evaluation
       * is needed when x is small (which is precisely
       * when it is easy to do).
       */
      const double scale_factor = GSL_SQRT_DBL_MAX;
      int scale_count_for = 0;
      int scale_count_bck = 0;
      int a0 = 1;
      int a1 = a0 + ceil(0.5*(b-x) - a0);
      double Ua1_bck_val;
      double Ua1_bck_err;
      double Ua1_for_val;
      double Ua1_for_err;
      int stat_for;
      int stat_bck;
      gsl_sf_result lm_for;

      {
        /* Recurse back to determine U(a1,b), sans normalization.
         */
        double ru;
	int CF1_count;
        int stat_CF1 = hyperg_U_CF1(a, b, 0, x, &ru, &CF1_count);
        double Ua   = 1.0;
        double Uap1 = ru/a * Ua;
        double Uam1;
        int ap;
        for(ap=a; ap>a1; ap--) {
          Uam1 = -((b-2.0*ap-x)*Ua + ap*(1.0+ap-b)*Uap1);
          Uap1 = Ua;
          Ua   = Uam1;
	  RESCALE_2(Ua,Uap1,scale_factor,scale_count_bck);
        }
        Ua1_bck_val = Ua;
	Ua1_bck_err = 2.0 * GSL_DBL_EPSILON * (fabs(a1-a)+CF1_count+1.0) * fabs(Ua);
        stat_bck = stat_CF1;
      }

      if(b == 2*a1 && a1 > 1) {
        /* This can happen when x is small, which is
	 * precisely when we need to be careful with
	 * this evaluation.
	 */
	hyperg_lnU_beq2a((double)a1, x, &lm_for);
	Ua1_for_val = 1.0;
	Ua1_for_err = 0.0;
        stat_for = GSL_SUCCESS;
      }
      else if(b == 2*a1 - 1 && a1 > 1) {
        /* Similar to the above. Happens when x is small.
	 * Use
	 *   U(a,2a-1) = (x U(a,2a) - U(a-1,2(a-1))) / (2a - 2)
	 */
	gsl_sf_result lnU00, lnU12;
	gsl_sf_result U00, U12;
	hyperg_lnU_beq2a(a1-1.0, x, &lnU00);
	hyperg_lnU_beq2a(a1,	 x, &lnU12);
	if(lnU00.val > lnU12.val) {
	  lm_for.val = lnU00.val;
	  lm_for.err = lnU00.err;
	  U00.val = 1.0;
	  U00.err = 0.0;
	  gsl_sf_exp_err_e(lnU12.val - lm_for.val, lnU12.err + lm_for.err, &U12);
	}
	else {
	  lm_for.val = lnU12.val;
	  lm_for.err = lnU12.err;
	  U12.val = 1.0;
	  U12.err = 0.0;
	  gsl_sf_exp_err_e(lnU00.val - lm_for.val, lnU00.err + lm_for.err, &U00);
	}
	Ua1_for_val  = (x * U12.val - U00.val) / (2.0*a1 - 2.0);
	Ua1_for_err  = (fabs(x)*U12.err + U00.err) / fabs(2.0*a1 - 2.0);
	Ua1_for_err += 2.0 * GSL_DBL_EPSILON * fabs(Ua1_for_val);
	stat_for = GSL_SUCCESS;
      }
      else {
        /* Recurse forward to determine U(a1,b) with
         * absolute normalization.
         */
	gsl_sf_result r_Ua;
        double Uam1 = 1.0;  /* U(a0-1,b,x) = U(0,b,x) */
        double Ua;
        double Uap1;
        int ap;
	double lm_for_local;
        stat_for = hyperg_U_small_a_bgt0(a0, b, x, &r_Ua, &lm_for_local); /* U(1,b,x) */
	Ua = r_Ua.val;
        Uam1 *= exp(-lm_for_local);
	lm_for.val = lm_for_local;
	lm_for.err = 0.0;

        for(ap=a0; ap<a1; ap++) {