[WP34S] Eagerly Anticipating Non-Regularized Incomplete Gammas


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.


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.


I just saw the recently posted code. You anticipated my rambling :) You did exactly what I suggested before I suggested it!

Look forward to the next build.


Edited: 13 May 2012, 6:36 p.m.


Les, might I suggest looking at the source code before raising all these (predominately needless) concerns?

I'm not calculating the LnGamma for the non-regularised versions. I only even calculate it in the regularised case. I'm not naively dividing to scale. If the non-regularised answer requires inversion, I calculate Gamma then and only then. This introduces a risk of overflow and maybe cancellation since gamma gets large rather fast but this is unavoidable at this point.

The normal distribution will be able to use the non-regularised gamma and divide by the constant scaling factor. No problems here. This will result in a speed up half the time.

I am not going to expose the continued fraction and series functions separately -- they are way too arcane even for me.

The only thing we won't get is sharing the regularisation factor when it is needed for the inversion or when doing a chi2 inverse. I can live with these inefficiencies.

Now, I should give you subversion access so you can fix the xrom code :-)

- Pauli


I did!

I downloaded only the emulator build to test what looked even to me to be a bug in the C-code.

Lines 2183 and 2187 should call gser and gcf with plga as the last parameter. The existing code continues to work great with &lga as that parameter for the regularized version, since it is computed up at line 2179, where it is also equated to plga that is subsequently not used.

However, when plga is set to NULL in the non-regularized case, that is not passed to gcf() and gser() and the errors are returned.

But I bet you caught this already!


Edited: 13 May 2012, 7:09 p.m.


Yeah, fixed this. Not committed or built yet.

Still not getting the correct answers from the normal cdf though. Must be something else wrong :-(

- Pauli


Found the other bug and a new firmware/emulator will be built soon.

- Pauli


Using a 34 digits constant for gamma(0.5) will be less accurate than calculating gamma(0.5) in 39 digit precision, so we'll be introducing another small source of error into the normal cdf. I doubt it will make much difference though.

- Pauli

Possibly Related Threads...
Thread Author Replies Views Last Post
  [WP34S] WP34S firmware on the AT91SAM7L-STK dev kit? jerome ibanes 1 883 10-04-2012, 04:59 PM
Last Post: Paul Dale
  [WP34S] Improving the Incomplete Gammas Les Wright 19 3,683 05-16-2012, 05:15 PM
Last Post: Marcus von Cube, Germany
  [wp34s] Incomplete Gamma on the wp34s Les Wright 18 3,688 12-06-2011, 11:07 AM
Last Post: Namir
  [wp34s] Romberg Integration on the wp34s Les Wright 2 1,126 12-04-2011, 10:49 PM
Last Post: Les Wright
  Incomplete gamma function et al. for 35s Les Wright 8 1,686 08-27-2007, 10:21 PM
Last Post: Karl Schneider

Forum Jump: