[Rd] pgamma discontinuity (PR#7307)

[terra at gnome.org]

For the record, this is what I am going to use in Gnumeric for now.

Morten



/*
 * Problems fixed.
 *
 * 1. x = -inf, scale = +inf: result now 0, was nan.
 *
 * 2. No discontinuity for very small alph.
 *      Old pgamma(1,1e-8,1,FALSE,TRUE) --> -19.9376 [ok]
 *      Old pgamma(1,0.9999999999e-8,1,FALSE,TRUE) --> -5556027.755 [bad]
 *    That's more than <dr_evil>two million</dr_evil> orders of magnitude caused
 *    by the pnorm approximation being used in an area it was not intended for.
 *
 * 3. No discontinuity at alph=100000.
 *      Old pgamma(1e6,100000,1,FALSE,TRUE) --> -669750.36332851 [ok]
 *      Old pgamma(1e6,100001,1,FALSE,TRUE) --> -599731.36192347 [bad]
 *
 * 4. Improved precision in calculating log of series sum using log1p.
 *
 * 5. Avoid using series and continued fraction too close to the critical
 *    line.
 *
 * 6. Improved precision of argument in pnorm approximation.
 *
 * 7. Use a power of 2 for rescaling in continued fraction so rescaling does
 *    not introduce any rounding errors.
 */





/*
 * Abramowitz and Stegun 6.5.29
 */
static double
pgamma_series_sum (double x, double alph)
{
     double sum = 0, c = 1, a = alph;

     do {
	  a += 1.;
	  c *= x / a;
	  sum += c;
     } while (c > DBL_EPSILON * sum);

     return sum;
}

/*
 * Abramowitz and Stegun 6.5.31
 */
static double
pgamma_cont_frac (double x, double alph)
{
     double a, b, n, an, pn1, pn2, pn3, pn4, pn5, pn6, sum, osum;
     double xlarge = pow (2, 128);

     a = 1. - alph;
     b = a + x + 1.;
     pn1 = 1.;
     pn2 = x;
     pn3 = x + 1.;
     pn4 = x * b;
     sum = pn3 / pn4;
     for (n = 1; ; n++) {
	  a += 1.;/* =   n+1 -alph */
	  b += 2.;/* = 2(n+1)-alph+x */
	  an = a * n;
	  pn5 = b * pn3 - an * pn1;
	  pn6 = b * pn4 - an * pn2;
	  if (fabs(pn6) > 0.) {
	       osum = sum;
	       sum = pn5 / pn6;
	       if (fabs(osum - sum) <= DBL_EPSILON * fmin2(1., sum))
		    break;
	  }
	  pn1 = pn3;
	  pn2 = pn4;
	  pn3 = pn5;
	  pn4 = pn6;
	  if (fabs(pn5) >= xlarge) {
	       /* re-scale the terms in continued fraction if they are large */
#ifdef DEBUG_p
	       REprintf(" [r] ");
#endif
	       pn1 /= xlarge;
	       pn2 /= xlarge;
	       pn3 /= xlarge;
	       pn4 /= xlarge;
	  }
     }

     return sum;
}


double pgamma(double x, double alph, double scale, int lower_tail, int log_p)
{
     double arg;

#ifdef IEEE_754
     if (ISNAN(x) || ISNAN(alph) || ISNAN(scale))
	  return x + alph + scale;
#endif
     if(alph <= 0. || scale <= 0.)
	  ML_ERR_return_NAN;
     if (x <= 0.)
	  return R_DT_0;
     x /= scale;
#ifdef IEEE_754
     if (ISNAN(x)) /* eg. original x = scale = +Inf */
	  return x;
#endif

     if (x < 0.99 * (alph + 1000) && x < 10 + alph && x < 10 * alph) {
	  /*
	   * The first condition makes sure that terms in the sum get at
	   * least 1% smaller, no later than after 1000 terms.
	   *
	   * The second condition makes sure that the terms start getting
	   * smaller (even if just a tiny bit) after no more than 10 terms.
	   *
	   * The third condition makes sure that the terms don't get huge
	   * before they start getting smaller.
	   */
	  double C1mls2p = /* 1 - log (sqrt (2 * M_PI)) */
	       0.081061466795327258219670263594382360138602526362216587182846;
	  double logsum = log1p (pgamma_series_sum (x, alph));
	  /*
	   * The first two terms here would cause cancellation for extremely
	   * large and near-equal x and alph.  The first condition above
	   * prevents that.
	   */
	  arg = (alph - x) + alph * log (x / (alph + 1)) + C1mls2p -
	       log1p (alph) / 2 - stirlerr (alph + 1) + logsum; 
     } else if (alph < x && alph - 100 < 0.99 * x) {
	  /*
	   * The first condition guarantees that we will start making
	   * a tiny bit of progress right now.
	   *
	   * The second condition guarantees that we will make decent
	   * progress after no more than 100 loops.
	   */
	  double lfrac = log (pgamma_cont_frac (x, alph));
	  arg = lfrac + alph * log(x) - x - lgammafn(alph);
	  lower_tail = !lower_tail;
     } else {
	  /*
	   * Approximation using pnorm.  We only get here for fairly large
	   * x and alph that are within 1% of each other.
	   */
	  double r3m1 = expm1 (log (x / alph) / 3);
	  return pnorm (sqrt (alph) * 3 * (r3m1 + 1 / (9 * alph)),
			0, 1, lower_tail, log_p);
     }

     if (lower_tail)
	  return log_p ? arg : exp(arg);
     else
	  return log_p ? R_Log1_Exp(arg) : -expm1(arg);
}