hyperg_1f1.c

开放gsl矩阵运算

      double pair_ratio;
      double rat_0 = fabs(r_Mnm1.err/r_Mnm1.val);
      double rat_1 = fabs(r_Mn.err/r_Mn.val);
      for(n=eps+1.0; n<a-0.1; n++) {
        Mnp1 = ((b-n)*Mnm1 + (2*n-b+x)*Mn)/n;
        Mnm1 = Mn;
        Mn   = Mnp1;
        minim_pair = GSL_MIN_DBL(fabs(Mn) + fabs(Mnm1), minim_pair);
      }
      pair_ratio = start_pair/minim_pair;
      result->val  = Mn;
      result->err  = 2.0 * (rat_0 + rat_1 + GSL_DBL_EPSILON) * (fabs(a)+1.0) * fabs(Mn);
      result->err += 2.0 * (rat_0 + rat_1) * pair_ratio*pair_ratio * fabs(Mn);
      result->err += 2.0 * GSL_DBL_EPSILON * fabs(Mn);
      return GSL_ERROR_SELECT_2(stat_0, stat_1);
    }
  }
  else {
    /* x < 0
     * b < a
     */

    if(a <= 0.5*(b-x) || a >= -x) {
      /* Recurse down in b, from near the a=b line, b=a+eps,a+eps-1.
       */
      double N   = floor(a - b);
      double eps = 1.0 + N - a + b;
      gsl_sf_result r_Manp1;
      gsl_sf_result r_Man;
      int stat_0 = hyperg_1F1_beps_bgt0(-eps,    a+eps,     x, &r_Manp1);
      int stat_1 = hyperg_1F1_beps_bgt0(1.0-eps, a+eps-1.0, x, &r_Man);
      double Manp1 = r_Manp1.val;
      double Man   = r_Man.val;
      double Manm1;

      double n;
      double start_pair = fabs(Manp1) + fabs(Man);
      double minim_pair = GSL_DBL_MAX;
      double pair_ratio;
      double rat_0 = fabs(r_Manp1.err/r_Manp1.val);
      double rat_1 = fabs(r_Man.err/r_Man.val);
      for(n=a+eps-1.0; n>b+0.1; n -= 1.0) {
        Manm1 = (-n*(1-n-x)*Man - x*(n-a)*Manp1)/(n*(n-1.0));
        Manp1 = Man;
        Man = Manm1;
        minim_pair = GSL_MIN_DBL(fabs(Manp1) + fabs(Man), minim_pair);
      }

      /* FIXME: this is a nasty little hack; there is some
         (transient?) instability in this recurrence for some
	 values. I can tell when it happens, which is when
	 this pair_ratio is large. But I do not know how to
	 measure the error in terms of it. I guessed quadratic
	 below, but it is probably worse than that.
	 */
      pair_ratio = start_pair/minim_pair;
      result->val  = Man;
      result->err  = 2.0 * (rat_0 + rat_1 + GSL_DBL_EPSILON) * (fabs(b-a)+1.0) * fabs(Man);
      result->err *= pair_ratio*pair_ratio + 1.0;
      return GSL_ERROR_SELECT_2(stat_0, stat_1);
    }
    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.
       */
      double epsa = a - floor(a);
      double a0   = floor(0.5*(b-x)) + epsa;
      double N    = floor(a0 - b);
      double epsb = 1.0 + N - a0 + b;
      double Ma0b;
      double Ma0bp1;
      double Ma0p1b;
      int stat_a0;
      double Mnm1;
      double Mn;
      double Mnp1;
      double n;
      double err_rat;
      {
    	gsl_sf_result r_Ma0np1;
    	gsl_sf_result r_Ma0n;
        int stat_0 = hyperg_1F1_beps_bgt0(-epsb,    a0+epsb,     x, &r_Ma0np1);
        int stat_1 = hyperg_1F1_beps_bgt0(1.0-epsb, a0+epsb-1.0, x, &r_Ma0n);
    	double Ma0np1 = r_Ma0np1.val;
    	double Ma0n   = r_Ma0n.val;
    	double Ma0nm1;

	err_rat = fabs(r_Ma0np1.err/r_Ma0np1.val) + fabs(r_Ma0n.err/r_Ma0n.val);

    	for(n=a0+epsb-1.0; n>b+0.1; n -= 1.0) {
    	  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); /* right-down hook */
	stat_a0 = GSL_ERROR_SELECT_2(stat_0, stat_1);
      }

          
      /* Initialise the recurrence correctly BJG */

      if (a0 >= a - 0.1)
        { 
          Mn = Ma0b;
        }
      else if (a0 + 1>= a - 0.1)
        {
          Mn = Ma0p1b;
        }
      else
        {
          Mnm1 = Ma0b;
          Mn   = Ma0p1b;

          for(n=a0+1.0; n<a-0.1; n += 1.0) {
            Mnp1 = ((b-n)*Mnm1 + (2*n-b+x)*Mn)/n;
            Mnm1 = Mn;
            Mn   = Mnp1;
          }
        }

      result->val = Mn;
      result->err = (err_rat + GSL_DBL_EPSILON) * (fabs(b-a)+1.0) * fabs(Mn);
      return stat_a0;
    }
  }
}


/* Assumes b != integer
 * Assumes a != integer when x > 0
 * Assumes b-a != neg integer when x < 0
 */
static
int
hyperg_1F1_ab_neg(const double a, const double b, const double x,
                  gsl_sf_result * result)
{
  const double bma = b - a;
  const double abs_x = fabs(x);
  const double abs_a = fabs(a);
  const double abs_b = fabs(b);
  const double size_a = GSL_MAX(abs_a, 1.0);
  const double size_b = GSL_MAX(abs_b, 1.0);
  const int bma_integer = ( bma - floor(bma+0.5) < _1F1_INT_THRESHOLD );

  if(   (abs_a < 10.0 && abs_b < 10.0 && abs_x < 5.0)
     || (b > 0.8*GSL_MAX(fabs(a),1.0)*fabs(x))
    ) {
    return gsl_sf_hyperg_1F1_series_e(a, b, x, result);
  }
  else if(   x > 0.0
          && size_b > size_a
          && size_a*log(M_E*x/size_b) < GSL_LOG_DBL_EPSILON+7.0
    ) {
    /* Series terms are positive definite up until
     * there is a sign change. But by then the
     * terms are small due to the last condition.
     */
    return gsl_sf_hyperg_1F1_series_e(a, b, x, result);
  }
  else if(   (abs_x < 5.0 && fabs(bma) < 10.0 && abs_b < 10.0)
          || (b > 0.8*GSL_MAX_DBL(fabs(bma),1.0)*abs_x)
    ) {
    /* Use Kummer transformation to render series safe.
     */
    gsl_sf_result Kummer_1F1;
    int stat_K = gsl_sf_hyperg_1F1_series_e(bma, b, -x, &Kummer_1F1);
    int stat_e = gsl_sf_exp_mult_err_e(x, GSL_DBL_EPSILON * fabs(x),
                                      Kummer_1F1.val, Kummer_1F1.err,
                                      result);
    return GSL_ERROR_SELECT_2(stat_e, stat_K);
  }
  else if(   x < -30.0
          && GSL_MAX_DBL(fabs(a),1.0)*GSL_MAX_DBL(fabs(1.0+a-b),1.0) < 0.99*fabs(x)
    ) {
    /* Large negative x asymptotic.
     * Note that we do not check if b-a is a negative integer.
     */
    return hyperg_1F1_asymp_negx(a, b, x, result);
  }
  else if(   x > 100.0
          && GSL_MAX_DBL(fabs(bma),1.0)*GSL_MAX_DBL(fabs(1.0-a),1.0) < 0.99*fabs(x)
    ) {
    /* Large positive x asymptotic.
     * Note that we do not check if a is a negative integer.
     */
    return hyperg_1F1_asymp_posx(a, b, x, result);
  }
  else if(x > 0.0 && !(bma_integer && bma > 0.0)) {
    return hyperg_1F1_U(a, b, x, result);
  }
  else {
    /* FIXME:  if all else fails, try the series... BJG */
    if (x < 0.0) {
      /* Apply Kummer Transformation */
      int status = gsl_sf_hyperg_1F1_series_e(b-a, b, -x, result);
      double K_factor = exp(x);
      result->val *= K_factor;
      result->err *= K_factor;
      return status;
    } else {
      int status = gsl_sf_hyperg_1F1_series_e(a, b, x, result);
      return status;
    }

    /* Sadness... */
    /* result->val = 0.0; */
    /* result->err = 0.0; */
    /* GSL_ERROR ("error", GSL_EUNIMPL); */
  }
}


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

int
gsl_sf_hyperg_1F1_int_e(const int a, const int b, const double x, gsl_sf_result * result)
{
  /* CHECK_POINTER(result) */

  if(x == 0.0) {
    result->val = 1.0;
    result->err = 0.0;
    return GSL_SUCCESS;
  }
  else if(a == b) {
    return gsl_sf_exp_e(x, result);
  }
  else if(b == 0) {
    DOMAIN_ERROR(result);
  }
  else if(a == 0) {
    result->val = 1.0;
    result->err = 0.0;
    return GSL_SUCCESS;
  }
  else if(b < 0 && (a < b || a > 0)) {
    /* Standard domain error due to singularity. */
    DOMAIN_ERROR(result);
  }
  else if(x > 100.0  && GSL_MAX_DBL(1.0,fabs(b-a))*GSL_MAX_DBL(1.0,fabs(1-a)) < 0.5 * x) {
    /* x -> +Inf asymptotic */
    return hyperg_1F1_asymp_posx(a, b, x, result);
  }
  else if(x < -100.0 && GSL_MAX_DBL(1.0,fabs(a))*GSL_MAX_DBL(1.0,fabs(1+a-b)) < 0.5 * fabs(x)) {
    /* x -> -Inf asymptotic */
    return hyperg_1F1_asymp_negx(a, b, x, result);
  }
  else if(a < 0 && b < 0) {
    return hyperg_1F1_ab_negint(a, b, x, result);
  }
  else if(a < 0 && b > 0) {
    /* Use Kummer to reduce it to the positive integer case.
     * Note that b > a, strictly, since we already trapped b = a.
     */
    gsl_sf_result Kummer_1F1;
    int stat_K = hyperg_1F1_ab_posint(b-a, b, -x, &Kummer_1F1);
    int stat_e = gsl_sf_exp_mult_err_e(x, GSL_DBL_EPSILON * fabs(x),
                                      Kummer_1F1.val, Kummer_1F1.err,
                                      result); 
    return GSL_ERROR_SELECT_2(stat_e, stat_K);
  }
  else {
    /* a > 0 and b > 0 */
    return hyperg_1F1_ab_posint(a, b, x, result);
  }
}


int
gsl_sf_hyperg_1F1_e(const double a, const double b, const double x,
                       gsl_sf_result * result
                       )
{
  const double bma = b - a;
  const double rinta = floor(a + 0.5);
  const double rintb = floor(b + 0.5);
  const double rintbma = floor(bma + 0.5);
  const int a_integer   = ( fabs(a-rinta) < _1F1_INT_THRESHOLD && rinta > INT_MIN && rinta < INT_MAX );
  const int b_integer   = ( fabs(b-rintb) < _1F1_INT_THRESHOLD && rintb > INT_MIN && rintb < INT_MAX );
  const int bma_integer = ( fabs(bma-rintbma) < _1F1_INT_THRESHOLD && rintbma > INT_MIN && rintbma < INT_MAX );
  const int b_neg_integer   = ( b < -0.1 && b_integer );
  const int a_neg_integer   = ( a < -0.1 && a_integer );
  const int bma_neg_integer = ( bma < -0.1 &&  bma_integer );

  /* CHECK_POINTER(result) */

  if(x == 0.0) {
    /* Testing for this before testing a and b
     * is somewhat arbitrary. The result is that
     * we have 1F1(a,0,0) = 1.
     */
    result->val = 1.0;
    result->err = 0.0;
    return GSL_SUCCESS;
  }
  else if(b == 0.0) {
    DOMAIN_ERROR(result);
  }
  else if(a == 0.0) {
    result->val = 1.0;
    result->err = 0.0;
    return GSL_SUCCESS;
  }
  else if(a == b) {
    /* case: a==b; exp(x)
     * It's good to test exact equality now.
     * We also test approximate equality later.
     */
    return gsl_sf_exp_e(x, result);
  }
  else if(fabs(b) < _1F1_INT_THRESHOLD) {
    /* Note that neither a nor b is zero, since
     * we eliminated that with the above tests.
     */
    if(fabs(a) < _1F1_INT_THRESHOLD) {
      /* a and b near zero: 1 + a/b (exp(x)-1)
       */
      gsl_sf_result exm1;
      int stat_e = gsl_sf_expm1_e(x, &exm1);
      double sa = ( a > 0.0 ? 1.0 : -1.0 );
      double sb = ( b > 0.0 ? 1.0 : -1.0 );
      double lnab = log(fabs(a/b)); /* safe */
      gsl_sf_result hx;
      int stat_hx = gsl_sf_exp_mult_err_e(lnab, GSL_DBL_EPSILON * fabs(lnab),
        				     sa * sb * exm1.val, exm1.err,
        				     &hx);
      result->val = (hx.val == GSL_DBL_MAX ? hx.val : 1.0 + hx.val);  /* FIXME: excessive paranoia ? what is DBL_MAX+1 ?*/
      result->err = hx.err;
      return GSL_ERROR_SELECT_2(stat_hx, stat_e);
    }
    else {
      /* b near zero and a not near zero
       */
      const double m_arg = 1.0/(0.5*b);
      gsl_sf_result F_renorm;
      int stat_F = hyperg_1F1_renorm_b0(a, x, &F_renorm);
      int stat_m = gsl_sf_multiply_err_e(m_arg, 2.0 * GSL_DBL_EPSILON * m_arg,
                                            0.5*F_renorm.val, 0.5*F_renorm.err,
				            result);
      return GSL_ERROR_SELECT_2(stat_m, stat_F);
    }
  }
  else if(a_integer && b_integer) {
    /* Check for reduction to the integer case.
     * Relies on the arbitrary "near an integer" test.
     */
    return gsl_sf_hyperg_1F1_int_e((int)rinta, (int)rintb, x, result);
  }
  else if(b_neg_integer && !(a_neg_integer && a > b)) {
    /* Standard domain error due to
     * uncancelled singularity.
     */
    DOMAIN_ERROR(result);
  }
  else if(a_neg_integer) {
    return hyperg_1F1_a_negint_lag((int)rinta, b, x, result);
  }
  else if(b > 0.0) {
    if(-1.0 <= a && a <= 1.0) {
      /* Handle small a explicitly.
       */
      return hyperg_1F1_small_a_bgt0(a, b, x, result);
    }
    else if(bma_neg_integer) {
      /* Catch this now, to avoid problems in the
       * generic evaluation code.
       */
      gsl_sf_result Kummer_1F1;
      int stat_K = hyperg_1F1_a_negint_lag((int)rintbma, b, -x, &Kummer_1F1);
      int stat_e = gsl_sf_exp_mult_err_e(x, GSL_DBL_EPSILON * fabs(x),
        				    Kummer_1F1.val, Kummer_1F1.err,
        				    result);
      return GSL_ERROR_SELECT_2(stat_e, stat_K);
    }
    else if(a < 0.0) {
      /* Use Kummer to reduce it to the generic positive case.
       * Note that b > a, strictly, since we already trapped b = a.
       * Also b-(b-a)=a, and a is not a negative integer here,
       * so the generic evaluation is safe.
       */
      gsl_sf_result Kummer_1F1;
      int stat_K = hyperg_1F1_ab_pos(b-a, b, -x, &Kummer_1F1);
      int stat_e = gsl_sf_exp_mult_err_e(x, GSL_DBL_EPSILON * fabs(x),
        				    Kummer_1F1.val, Kummer_1F1.err,
        				    result);
      return GSL_ERROR_SELECT_2(stat_e, stat_K);
    }
    else {
      /* a > 0.0 */
      return hyperg_1F1_ab_pos(a, b, x, result);
    }
  }
  else {
    /* b < 0.0 */

    if(bma_neg_integer && x < 0.0) {
      /* Handle this now to prevent problems
       * in the generic evaluation.
       */
      gsl_sf_result K;
      int stat_K;
      int stat_e;
      if(a < 0.0) {
        /* Kummer transformed version of safe polynomial.
	 * The condition a < 0 is equivalent to b < b-a,
	 * which is the condition required for the series
	 * to be positive definite here.
         */
        stat_K = hyperg_1F1_a_negint_poly((int)rintbma, b, -x, &K);
      }
      else {
        /* Generic eval for negative integer a. */
	stat_K = hyperg_1F1_a_negint_lag((int)rintbma, b, -x, &K);
      }
      stat_e = gsl_sf_exp_mult_err_e(x, GSL_DBL_EPSILON * fabs(x),
                                        K.val, K.err,
                                        result);
      return GSL_ERROR_SELECT_2(stat_e, stat_K);
    }
    else if(a > 0.0) {
      /* Use Kummer to reduce it to the generic negative case.
       */
      gsl_sf_result K;
      int stat_K = hyperg_1F1_ab_neg(b-a, b, -x, &K);
      int stat_e = gsl_sf_exp_mult_err_e(x, GSL_DBL_EPSILON * fabs(x),
        				    K.val, K.err,
        				    result);
      return GSL_ERROR_SELECT_2(stat_e, stat_K);
    }
    else {
      return hyperg_1F1_ab_neg(a, b, x, result);
    }
  }
}


  
#if 0  
    /* Luke in the canonical case.
   */
  if(x < 0.0 && !a_neg_integer && !bma_neg_integer) {
    double prec;
    return hyperg_1F1_luke(a, b, x, result, &prec);
  }


  /* Luke with Kummer transformation.
   */
  if(x > 0.0 && !a_neg_integer && !bma_neg_integer) {
    double prec;
    double Kummer_1F1;
    double ex;
    int stat_F = hyperg_1F1_luke(b-a, b, -x, &Kummer_1F1, &prec);
    int stat_e = gsl_sf_exp_e(x, &ex);
    if(stat_F == GSL_SUCCESS && stat_e == GSL_SUCCESS) {
      double lnr = log(fabs(Kummer_1F1)) + x;
      if(lnr < GSL_LOG_DBL_MAX) {
        *result = ex * Kummer_1F1;
	return GSL_SUCCESS;
      }
      else {
        *result = GSL_POSINF; 
	GSL_ERROR ("overflow", GSL_EOVRFLW);
      }
    }
    else if(stat_F != GSL_SUCCESS) {
      *result = 0.0;
      return stat_F;
    }
    else {
      *result = 0.0;
      return stat_e;
    }
  }
#endif



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

#include "eval.h"

double gsl_sf_hyperg_1F1_int(const int m, const int n, double x)
{
  EVAL_RESULT(gsl_sf_hyperg_1F1_int_e(m, n, x, &result));
}

double gsl_sf_hyperg_1F1(double a, double b, double x)
{
  EVAL_RESULT(gsl_sf_hyperg_1F1_e(a, b, x, &result));
}