HP Articles Forum
[Return to the Index ]
[ Previous | Next ]
Incomplete Gamma Function for wp34s
Posted by Les Wright on 5 Dec 2011, 2:40 a.m.
Most readers will know that the "complete" gamma function is related to the factorial function by the following: Gamma(a) = (a-1)!. For real numbers, this function is computed by integrating t^(a-1)*exp(-t) from t = 0 to t = +infinity.
The so-called incomplete gamma functions are what one gets when one computes this integration with one of the limits of integration changed to something different. For example, the left-sided incomplete gamma, which here is computed by IGS, is given when one takes the integration from t = 0 to t = some x. Analogously the right-sided incomplete gamma function, which here is computed by IGL, takes the integration from t = some x to t = +infinity. Obviously, the sum of IGS and IGL for the same underlying "a" gives the "complete" gamma for that "a".
The left-sided incomplete gamma is best computed by a series approximation when x < a, so the routine is call IGS, where the "S" stands for "series". The right-sided incomplete gamma is best computed by a continued fraction expansion using the modified Lentz method when x > a, so is called IGL, where the "L" is for Lentz.
The code below was adapted from the routines given in the Numerical Recipes books:
001 LBL'IGS' 002 STO 00 003 x[<->] Y 004 STO 01 005 STO 02 006 1/x 007 STO 03 008 STO 04 009 LBL 00 010 STO 05 011 INC 02 012 RCL 00 013 STO[times] 03 014 RCL 02 015 STO/ 03 016 RCL 03 017 STO+ 04 018 RCL 05 019 RCL 04 020 x[!=]? Y 021 GTO 00 022 GTO 04 023 LBL'IGL' 024 STO 00 025 x[<->] Y 026 STO 01 027 - 028 1 029 + 030 STO 02 031 1 032 EEX 033 2 034 5 035 0 036 STO 05 037 1/x 038 STO 06 039 + 040 1/x 041 STO 03 042 STO 04 043 CLx 044 STO 08 045 LBL 09 046 INC 08 047 INC 02 048 INC 02 049 RCL 01 050 RCL- 08 051 RCL[times] 08 052 STO 07 053 RCL[times] 03 054 RCL+ 02 055 RCL+ 06 056 STO 03 057 RCL 07 058 RCL/ 05 059 RCL+ 02 060 RCL+ 06 061 STO 05 062 RCL 03 063 STO/ 03 064 STO/ 03 065 / 066 STO[times] 04 067 x[!=]1? 068 GTO 09 069 RCL 04 070 LBL 04 071 RCL 00 072 RCL 01 073 y[^x] 074 [times] 075 RCL 00 076 +/- 077 e[^x] 078 [times] 079 RTN
Use is simple: put the arguments on the stack as a [ENTER] x and XEQ either IGS or IGL. Given the high internal accuracy of the wp34s, the displayed digits (and then some) are almost always completely accurate.
Example: Compute the left and right incomplete gammas for a = 5 and x =5.
5 [ENTER] 5 {XEQ] [alpha] IGS returns 13.4281611584 in ALL 00.
5 [ENTER] 5 {XEQ] [alpha] IGL returns 10.5718388416.
It can be easily seen that these outputs sum to 24, which is the "complete" gamma of 5, or the factorial of 5-1 = 4.
I am aware that this routine may be of only academic interest to most. The incomplete gamma functions are of practical use in that they form the basis of the error function, the cumulative standard normal distribution, and the cumulative chi-square distribution. All of these applications have their own excellent routines already built into the wp34s. I originally ported the Numerical Recipes program for use on the HP41C, HP42S, HP35s, and older calcs that did not have these statistical functions built in. Nonetheless, I hope this code interests someone.
Edited: 5 Dec 2011, 6:04 a.m.