This has turned up in the recent manual update, and of course I am eager how it will be implemented.

My suggestion is this: depending on the size of the x and y parameters (or the a and x parameters for those following along in A&S or the much-maligned Numerical Recipes), call gser() or gcf(), but passing zero as the gln parameter. If the result returned is what one wants, just return it. If what one wants is the complement (i.e. one wants the left side from the gcf result or the right side from the gser result), compute the complete gamma with a call to decNumberGamma() and do the subtraction and return that.

My concern is that if one does it the simplest way--computing the regularized incomplete gamma and the complete gamma and multiply the two together--one is computing the complete gamma a lot of times when one really doesn't need it. Moreover, computing an unnecessary LnGamma for gcf or gser and later a complete Gamma just to cancel it out seems redundant, and may compromise accuracy loss through rounding. If one hopes to speed up the inverse CDF by minimizing the need to repeatedly compute Gamma(0.5) = sqrt(Pi), this seems the way to do it.

Of course, my suggestion has subtraction error implications. In choosing between gser() or gcf() in any of these computations, the break point is x= y+1. If I use my user routines to compute RightSidedGamma(5,6) via the continued fraction, so that no computation of Gamma(5)=24 and no subtraction are required, I get a 34-digit result that is off just 1ULP. When I compute a left-sided incomplete gamma with (5,6), I get a 34-digit result of the LEFT-side that is again off just 1ULP. When I subtract this from Gamma(5) = 24, I effectively lose a digit as the right-sided gamma is off by 7ULP in the 34th digit. I am sure that more pathological examples can be found. The good news is that if all of this is done to 39 digits in the C code, this subtraction loss may not be detected in typical use. However, I do think the rounding error implications are worse than subtraction error. When I use the calculator's left and right regularized routines with the same arguments and multiply through by 24, I get a left-sided result (before any subtraction) that is off by 7ULP, and a right-sided result that is off by 53ULP. My humble user routines were inspired by the very same Numerical Recipes code yours were.

My personal fantasy is that gser() and gcf() do not accept a gln parameter and do not return a regularized result, deferring that decision to other routines downstream. Indeed, my own user code consists of such routines--IGS ("incomplete gamma series") and IGL ("incomplete gamma Lentz") that do just that. Therefore, when I use them in my own versions of the upper and lower cumulative normal distributions, things seem faster in many cases than with with built-in routine--I think because I regularize the gammas just by dividing through by sqrt(Pi). I don't need to compute Gamma(0.5) like the calculator's existing routines routines are obliged to do.

Just more musing from me.

Les

P.S., all of the above regarding the normal distribution applies to the use of the incomplete gammas to compute erf and erfc, which also regularize the relevant incomplete gammas by Gamma(0.5) = sqrt(Pi).

*Edited: 13 May 2012, 6:25 p.m. *