hyperg_u.c

开放gsl矩阵运算

          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_for);
        }
        Ua1_for_val  = Ua;
	Ua1_for_err  = fabs(Ua) * fabs(r_Ua.err/r_Ua.val);
	Ua1_for_err += 2.0 * GSL_DBL_EPSILON * (fabs(a1-a0)+1.0) * fabs(Ua1_for_val);
      }

      /* Now do the matching to produce the final result.
       */
      if(Ua1_bck_val == 0.0) {
        result->val = 0.0;
	result->err = 0.0;
	result->e10 = 0;
        GSL_ERROR ("error", GSL_EZERODIV);
      }
      else if(Ua1_for_val == 0.0) {
        /* Should never happen. */
        UNDERFLOW_ERROR_E10(result);
      }
      else {
        double lns = (scale_count_for - scale_count_bck)*log(scale_factor);
	double ln_for_val = log(fabs(Ua1_for_val));
	double ln_for_err = GSL_DBL_EPSILON + fabs(Ua1_for_err/Ua1_for_val);
	double ln_bck_val = log(fabs(Ua1_bck_val));
	double ln_bck_err = GSL_DBL_EPSILON + fabs(Ua1_bck_err/Ua1_bck_val);
        double lnr_val = lm_for.val + ln_for_val - ln_bck_val + lns;
	double lnr_err = lm_for.err + ln_for_err + ln_bck_err
	  + 2.0 * GSL_DBL_EPSILON * (fabs(lm_for.val) + fabs(ln_for_val) + fabs(ln_bck_val) + fabs(lns));
	double sgn = GSL_SIGN(Ua1_for_val) * GSL_SIGN(Ua1_bck_val);
        int stat_e = gsl_sf_exp_err_e10_e(lnr_val, lnr_err, result);
	result->val *= sgn;
        return GSL_ERROR_SELECT_3(stat_e, stat_bck, stat_for);
      }
    }
  }
}


/* Handle b >= 1 for generic a,b values.
 */
static
int
hyperg_U_bge1(const double a, const double b, const double x,
              gsl_sf_result_e10 * result)
{
  const double rinta = floor(a+0.5);
  const int a_neg_integer = (a < 0.0 && fabs(a - rinta) < INT_THRESHOLD);

  if(a == 0.0) {
    result->val = 1.0;
    result->err = 0.0;
    result->e10 = 0;
    return GSL_SUCCESS;
  }
  else if(a_neg_integer && fabs(rinta) < INT_MAX) {
    /* U(-n,b,x) = (-1)^n n! Laguerre[n,b-1,x]
     */
    const int n = -(int)rinta;
    const double sgn = (GSL_IS_ODD(n) ? -1.0 : 1.0);
    gsl_sf_result lnfact;
    gsl_sf_result L;
    const int stat_L = gsl_sf_laguerre_n_e(n, b-1.0, x, &L);
    gsl_sf_lnfact_e(n, &lnfact);
    {
      const int stat_e = gsl_sf_exp_mult_err_e10_e(lnfact.val, lnfact.err,
                                                      sgn*L.val, L.err,
                                                      result);
      return GSL_ERROR_SELECT_2(stat_e, stat_L);
    }
  }
  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(fabs(a) <= 1.0) {
    gsl_sf_result rU;
    double ln_multiplier;
    int stat_U = hyperg_U_small_a_bgt0(a, b, x, &rU, &ln_multiplier);
    int stat_e = gsl_sf_exp_mult_err_e10_e(ln_multiplier, 2.0*GSL_DBL_EPSILON*fabs(ln_multiplier), rU.val, rU.err, result);
    return GSL_ERROR_SELECT_2(stat_U, stat_e);
  }
  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.0) {
    /* Recurse backward on a and then upward on b.
     */
    const double scale_factor = GSL_SQRT_DBL_MAX;
    const double a0 = a - floor(a) - 1.0;
    const double b0 = b - floor(b) + 1.0;
    int scale_count = 0;
    double lm_0, lm_1;
    double lm_max;
    gsl_sf_result r_Uap1;
    gsl_sf_result r_Ua;
    int stat_0 = hyperg_U_small_a_bgt0(a0+1.0, b0, x, &r_Uap1, &lm_0);
    int stat_1 = hyperg_U_small_a_bgt0(a0,     b0, x, &r_Ua,   &lm_1);
    int stat_e;
    double Uap1 = r_Uap1.val;
    double Ua   = r_Ua.val;
    double Uam1;
    double ap;
    lm_max = GSL_MAX(lm_0, lm_1);
    Uap1 *= exp(lm_0-lm_max);
    Ua   *= exp(lm_1-lm_max);

    /* Downward recursion on a.
     */
    for(ap=a0; ap>a+0.1; ap -= 1.0) {
      Uam1 = ap*(b0-ap-1.0)*Uap1 + (x+2.0*ap-b0)*Ua;
      Uap1 = Ua;
      Ua   = Uam1;
      RESCALE_2(Ua,Uap1,scale_factor,scale_count);
    }

    if(b < 2.0) {
      /* b == b0, so no recursion necessary
       */
      const double lnscale = log(scale_factor);
      gsl_sf_result lnm;
      gsl_sf_result y;
      lnm.val = lm_max + scale_count * lnscale;
      lnm.err = 2.0 * GSL_DBL_EPSILON * (fabs(lm_max) + scale_count * fabs(lnscale));
      y.val  = Ua;
      y.err  = fabs(r_Uap1.err/r_Uap1.val) * fabs(Ua);
      y.err += fabs(r_Ua.err/r_Ua.val) * fabs(Ua);
      y.err += 2.0 * GSL_DBL_EPSILON * (fabs(a-a0) + 1.0) * fabs(Ua);
      y.err *= fabs(lm_0-lm_max) + 1.0;
      y.err *= fabs(lm_1-lm_max) + 1.0;
      stat_e = gsl_sf_exp_mult_err_e10_e(lnm.val, lnm.err, y.val, y.err, result);
    }
    else {
      /* Upward recursion on b.
       */
      const double err_mult = fabs(b-b0) + fabs(a-a0) + 1.0;
      const double lnscale = log(scale_factor);
      gsl_sf_result lnm;
      gsl_sf_result y;

      double Ubm1 = Ua; 				/* U(a,b0)   */
      double Ub   = (a*(b0-a-1.0)*Uap1 + (a+x)*Ua)/x;	/* U(a,b0+1) */
      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;
        RESCALE_2(Ub,Ubm1,scale_factor,scale_count);
      }

      lnm.val = lm_max + scale_count * lnscale;
      lnm.err = 2.0 * GSL_DBL_EPSILON * (fabs(lm_max) + fabs(scale_count * lnscale));
      y.val = Ub;
      y.err  = 2.0 * err_mult * fabs(r_Uap1.err/r_Uap1.val) * fabs(Ub);
      y.err += 2.0 * err_mult * fabs(r_Ua.err/r_Ua.val) * fabs(Ub);
      y.err += 2.0 * GSL_DBL_EPSILON * err_mult * fabs(Ub);
      y.err *= fabs(lm_0-lm_max) + 1.0;
      y.err *= fabs(lm_1-lm_max) + 1.0;
      stat_e = gsl_sf_exp_mult_err_e10_e(lnm.val, lnm.err, y.val, y.err, result);
    }
    return GSL_ERROR_SELECT_3(stat_e, stat_0, stat_1);
  }
  else if(b >= 2*a + x) {
    /* Recurse forward from a near zero.
     * Note that we cannot cross the singularity at
     * the line b=a+1, because the only way we could
     * be in that little wedge is if a < 1. But we
     * have already dealt with the small a case.
     */
    int scale_count = 0;
    const double a0 = a - floor(a);
    const double scale_factor = GSL_SQRT_DBL_MAX;
    double lnscale;
    double lm_0, lm_1, lm_max;
    gsl_sf_result r_Uam1;
    gsl_sf_result r_Ua;
    int stat_0 = hyperg_U_small_a_bgt0(a0-1.0, b, x, &r_Uam1, &lm_0);
    int stat_1 = hyperg_U_small_a_bgt0(a0,     b, x, &r_Ua,   &lm_1);
    int stat_e;
    gsl_sf_result lnm;
    gsl_sf_result y;
    double Uam1 = r_Uam1.val;
    double Ua   = r_Ua.val;
    double Uap1;
    double ap;
    lm_max = GSL_MAX(lm_0, lm_1);
    Uam1 *= exp(lm_0-lm_max);
    Ua   *= exp(lm_1-lm_max);

    for(ap=a0; ap<a-0.1; ap += 1.0) {
      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_max + scale_count * lnscale;
    lnm.err = 2.0 * GSL_DBL_EPSILON * (fabs(lm_max) + fabs(scale_count * lnscale));
    y.val  = Ua;
    y.err  = fabs(r_Uam1.err/r_Uam1.val) * fabs(Ua);
    y.err += fabs(r_Ua.err/r_Ua.val) * fabs(Ua);
    y.err += 2.0 * GSL_DBL_EPSILON * (fabs(a-a0) + 1.0) * fabs(Ua);
    y.err *= fabs(lm_0-lm_max) + 1.0;
    y.err *= fabs(lm_1-lm_max) + 1.0;
    stat_e = gsl_sf_exp_mult_err_e10_e(lnm.val, lnm.err, y.val, y.err, result);
    return GSL_ERROR_SELECT_3(stat_e, stat_0, stat_1);
  }
  else {
    if(b <= x) {
      /* Recurse backward to a near zero.
       */
      const double a0 = a - floor(a);
      const double scale_factor = GSL_SQRT_DBL_MAX;
      int scale_count = 0;
      gsl_sf_result lnm;
      gsl_sf_result y;
      double lnscale;
      double lm_0;
      double Uap1;
      double Ua;
      double Uam1;
      gsl_sf_result U0;
      double ap;
      double ru;
      double r;
      int CF1_count;
      int stat_CF1 = hyperg_U_CF1(a, b, 0, x, &ru, &CF1_count);
      int stat_U0;
      int stat_e;
      r = ru/a;
      Ua   = GSL_SQRT_DBL_MIN;
      Uap1 = r * Ua;
      for(ap=a; ap>a0+0.1; ap -= 1.0) {
        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);
      }

      stat_U0 = hyperg_U_small_a_bgt0(a0, b, x, &U0, &lm_0);

      lnscale = log(scale_factor);
      lnm.val = lm_0 - scale_count * lnscale;
      lnm.err = 2.0 * GSL_DBL_EPSILON * (fabs(lm_0) + fabs(scale_count * lnscale));
      y.val  = GSL_SQRT_DBL_MIN*(U0.val/Ua);
      y.err  = GSL_SQRT_DBL_MIN*(U0.err/fabs(Ua));
      y.err += 2.0 * GSL_DBL_EPSILON * (fabs(a0-a) + CF1_count + 1.0) * fabs(y.val);
      stat_e = gsl_sf_exp_mult_err_e10_e(lnm.val, lnm.err, y.val, y.err, result);
      return GSL_ERROR_SELECT_3(stat_e, stat_U0, stat_CF1);
    }
    else {
      /* Recurse backward to near the b=2a+x line, then
       * forward from a near zero to get the normalization.
       */
      int scale_count_for = 0;
      int scale_count_bck = 0;
      const double scale_factor = GSL_SQRT_DBL_MAX;
      const double eps = a - floor(a);
      const double a0 = ( eps == 0.0 ? 1.0 : eps );
      const double a1 = a0 + ceil(0.5*(b-x) - a0);
      gsl_sf_result lnm;
      gsl_sf_result y;
      double lm_for;
      double lnscale;
      double Ua1_bck;
      double Ua1_for;
      int stat_for;
      int stat_bck;
      int stat_e;
      int CF1_count;

      {
        /* Recurse back to determine U(a1,b), sans normalization.
         */
        double Uap1;
        double Ua;
        double Uam1;
        double ap;
        double ru;
        double r;
        int stat_CF1 = hyperg_U_CF1(a, b, 0, x, &ru, &CF1_count);
        r = ru/a;
        Ua   = GSL_SQRT_DBL_MIN;
        Uap1 = r * Ua;
        for(ap=a; ap>a1+0.1; ap -= 1.0) {
          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 = Ua;
        stat_bck = stat_CF1;
      }
      {
        /* Recurse forward to determine U(a1,b) with
         * absolute normalization.
         */
	gsl_sf_result r_Uam1;
	gsl_sf_result r_Ua;
	double lm_0, lm_1;
        int stat_0 = hyperg_U_small_a_bgt0(a0-1.0, b, x, &r_Uam1, &lm_0);
        int stat_1 = hyperg_U_small_a_bgt0(a0,     b, x, &r_Ua,   &lm_1);
        double Uam1 = r_Uam1.val;
        double Ua   = r_Ua.val;
        double Uap1;
        double ap;

	lm_for = GSL_MAX(lm_0, lm_1);
	Uam1 *= exp(lm_0 - lm_for);
	Ua   *= exp(lm_1 - lm_for);

        for(ap=a0; ap<a1-0.1; ap += 1.0) {
          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_for);
        }
        Ua1_for = Ua;
	stat_for = GSL_ERROR_SELECT_2(stat_0, stat_1);
      }

      lnscale = log(scale_factor);
      lnm.val = lm_for + (scale_count_for - scale_count_bck)*lnscale;
      lnm.err = 2.0 * GSL_DBL_EPSILON * (fabs(lm_for) + fabs(scale_count_for - scale_count_bck)*fabs(lnscale));
      y.val = GSL_SQRT_DBL_MIN*Ua1_for/Ua1_bck;
      y.err = 2.0 * GSL_DBL_EPSILON * (fabs(a-a0) + CF1_count + 1.0) * fabs(y.val);
      stat_e = gsl_sf_exp_mult_err_e10_e(lnm.val, lnm.err, y.val, y.err, result);
      return GSL_ERROR_SELECT_3(stat_e, stat_bck, stat_for);
    }
  }
}


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


int
gsl_sf_hyperg_U_int_e10_e(const int a, const int b, const double x,
                             gsl_sf_result_e10 * result)
{
  /* CHECK_POINTER(result) */

  if(x <= 0.0) {
    DOMAIN_ERROR_E10(result);
  }
  else {
    if(b >= 1) {
      return hyperg_U_int_bge1(a, b, x, result);
    }
    else {
      /* Use the reflection formula
       * U(a,b,x) = x^(1-b) U(1+a-b,2-b,x)
       */
      gsl_sf_result_e10 U;
      double ln_x = log(x);
      int ap = 1 + a - b;
      int bp = 2 - b;
      int stat_e;
      int stat_U = hyperg_U_int_bge1(ap, bp, x, &U);
      double ln_pre_val = (1.0-b)*ln_x;
      double ln_pre_err = 2.0 * GSL_DBL_EPSILON * (fabs(b)+1.0) * fabs(ln_x);
      ln_pre_err += 2.0 * GSL_DBL_EPSILON * fabs(1.0-b); /* error in log(x) */
      stat_e = gsl_sf_exp_mult_err_e10_e(ln_pre_val + U.e10*M_LN10, ln_pre_err,
                                            U.val, U.err,
                                            result);
      return GSL_ERROR_SELECT_2(stat_e, stat_U);
    }
  }
}


int
gsl_sf_hyperg_U_e10_e(const double a, const double b, const double x,
                         gsl_sf_result_e10 * result)
{
  const double rinta = floor(a + 0.5);
  const double rintb = floor(b + 0.5);
  const int a_integer = ( fabs(a - rinta) < INT_THRESHOLD );
  const int b_integer = ( fabs(b - rintb) < INT_THRESHOLD );

  /* CHECK_POINTER(result) */

  if(x <= 0.0) {
    DOMAIN_ERROR_E10(result);
  }
  else if(a == 0.0) {
    result->val = 1.0;
    result->err = 0.0;
    result->e10 = 0;
    return GSL_SUCCESS;
  }
  else if(a_integer && b_integer) {
    return gsl_sf_hyperg_U_int_e10_e(rinta, rintb, x, result);
  }
  else {
    if(b >= 1.0) {
      /* Use b >= 1 function.
       */
      return hyperg_U_bge1(a, b, x, result);
    }
    else {
      /* Use the reflection formula
       * U(a,b,x) = x^(1-b) U(1+a-b,2-b,x)
       */
      const double lnx = log(x);
      const double ln_pre_val = (1.0-b)*lnx;
      const double ln_pre_err = fabs(lnx) * 2.0 * GSL_DBL_EPSILON * (1.0 + fabs(b));
      const double ap = 1.0 + a - b;
      const double bp = 2.0 - b;
      gsl_sf_result_e10 U;
      int stat_U = hyperg_U_bge1(ap, bp, x, &U);
      int stat_e = gsl_sf_exp_mult_err_e10_e(ln_pre_val + U.e10*M_LN10, ln_pre_err,
                                            U.val, U.err,
                                            result);
      return GSL_ERROR_SELECT_2(stat_e, stat_U);
    }
  }
}


int
gsl_sf_hyperg_U_int_e(const int a, const int b, const double x, gsl_sf_result * result)
{
  gsl_sf_result_e10 re;
  int stat_U = gsl_sf_hyperg_U_int_e10_e(a, b, x, &re);
  int stat_c = gsl_sf_result_smash_e(&re, result);
  return GSL_ERROR_SELECT_2(stat_c, stat_U);
}


int
gsl_sf_hyperg_U_e(const double a, const double b, const double x, gsl_sf_result * result)
{
  gsl_sf_result_e10 re;
  int stat_U = gsl_sf_hyperg_U_e10_e(a, b, x, &re);
  int stat_c = gsl_sf_result_smash_e(&re, result);
  return GSL_ERROR_SELECT_2(stat_c, stat_U);
}


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

#include "eval.h"

double gsl_sf_hyperg_U_int(const int a, const int b, const double x)
{
  EVAL_RESULT(gsl_sf_hyperg_U_int_e(a, b, x, &result));
}

double gsl_sf_hyperg_U(const double a, const double b, const double x)
{
  EVAL_RESULT(gsl_sf_hyperg_U_e(a, b, x, &result));
}