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.

[ Return to the Message Index ]

Go back to the main exhibit hall