hyperg_1f1.c

math library from gnu

    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;
        }
        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);
      }

      /* Initialise the recurrence correctly BJG */

      if (a0 >= a)
        { 
          Mn = Ma0b;
        }
      else if (a0 + 1>= a)
        {
          Mn = Ma0p1b;
        }
      else
        {
          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,