[Rd] pgamma discontinuity (PR#7307)

[Morten Welinder]

> Make that 30400 orders of magnitude (natural logs y'know)...

Right.  (/me raises hands showing 2.7 fingers.)

> What the devil are you calculating? The probability that a random
> configuration of atoms would make up the known universe?

Not quite.  Where you see a cdf for the gamma distribution I see the
incomplete gamma function.  Same function, different hat.  I am using
it to compute the Erlang B function ("Grade Of Service"), see

    http://www.dcss.mcmaster.ca/~qiao/publications/erlang/newerlang.html

And here is my code for the log version of this.  (Link's c==circuit;
link's rho==traffic; link's p is a typo for rho.)

-----------------------------------------------------------------------------
static gnm_float
calculate_loggos (gnm_float traffic, gnm_float circuits)
{
	double f;

	if (traffic < 0 || circuits < 1)
		return gnm_nan;
	if (traffic == 0)
		return gnm_ninf;

#ifdef CANCELLATION
	/* Calculated this way we get cancellation.  */
	f = circuits * loggnum (traffic) - lgamma1p (circuits) - traffic;
#else
	f = (circuits - traffic) +
		(1 - loggnum (sqrtgnum (2 * M_PIgnum))) -
		loggnum (circuits + 1) / 2.0 -
		logfbit (circuits) +
		circuits * (loggnum (traffic / (circuits + 1)));
#endif

	return f - pgamma (traffic, circuits + 1, 1, FALSE, TRUE);
}
-----------------------------------------------------------------------------

The two #ifdef branches calculate the same thing, but the bottom
version suffers a lot less from cancellation.  I might still need
to consider cancellation in the final subtraction.  (Read "double"
where the above says "gnm_float" and forget the "gnum" suffixes.)

In the traffic=1e6,circuits=1e5 case I quoted I could use the
second formula from the link above instead, but that won't work
when the two are close to each other.  Sadly I need it there too.

Googling suggests that the canonical reference for this problem is

    Temme N M (1987) On the computation of the incomplete gamma
    functions for large values of the parameters Algorithms for
    Approximation J C Mason and M G Cox (ed) Oxford University Press.

(This reference from nag's manual.)

Morten