Incomplete Gamma Function

[Namir]

This program implements an incomplete Gamma function base don an HP-25 program by Peter Henrici.

Memory Map

R0 = X
R1 = A
R2 = L(n-2)
R3 = L(n-1)
R4 = s
R5 = gamma(n)
R6 = n
R7 = tolerance

Listing

Code:

1 LBL a
2 STO 7
3 R/S
4 LBL A
5 STO 1
6 X<>Y
7 STO 0
8 CLX
9 STO 2
10 STO 4
11 1
12 STO 3
13 STO 6
14 RCL 1
15 RCL 0
16 Y^X
17 STO 5
18 RCL 1
19 EXP
20 STO/ 5
21 LBL 0
22 RCL 5
23 RCL 6
24 RCL 0
25 -
26 STO* 5
27 1
28 -
29 RCL 3
30 RCL 2
31 -
32 *
33 RCL 6
34 RCL 1
35 +
36 RCL 3
37 STO 2
38 *
39 +
40 RCL 6
41 /
42 STO 3
43 RCL 2
44 *
45 /
46 STO+ 4
47 PAUSE
48 ABS
49 RCL 7
50 X<>Y
51 X<=Y?
52 GTO 1
53 1
54 STO+ 6
55 RCL 6
56 STO/ 5
57 GTO 0
58 LBL 1
59 RCL 4
60 R/S
61 GTO A

Usage


1. Enter the tolerance value and press [f][A]. The program automatically resumes with step 2.
2. Enter the values for A and X, and then press [A].
3. The program pauses to display the improvements in the calculated incomplete gamma function. When these improvements are less than the tolerance value, the program displays the result. If you press [R/S] after viewing the value of the incomplete gamma, the program flow resumes with step 2.

[bshoring]

Namir,
I wonder if you can enlighten me on this program.

I am no expert on Gamma nor a mathematician, but I have run other programs that compute both Gamma and Incomplete Gamma.

I have run this program as well as the one on page 223 of Henrici's book. Using the examples he gives on page 229, I get the same results with both programs IF I use a=0.

I am confused on:
1. What should one use for "a"? (I see a reference to a < 0 on p. 223)

2. Am I correct that the argument you want the result for is x?

3. With other incomplete gamma programs I have run, if I use a=0.5 and x=>35, I get a result close to 1.77245, comparable to an actual gamma function. I guess I am just trying to figure out how to relate the Henrici Incomplete Gamma to others.

Do you think you can shed any light on this for me? I would just like to be able to understand what I am getting.

Thanks,
Bob

[Dieter]

(03-26-2015 10:41 PM)bshoring Wrote:  Do you think you can shed any light on this for me? I would just like to be able to understand what I am getting.

There are different Gamma functions.

1. The "normal" Gamma function Γ(a). For integer a this is equal to (a–1)!. So Γ(5) = 4! = 24.
This, let's say "complete" Gamma function Γ(a) can be expressed as an integral from 0 to infinity.

2. The incomplete Gamma functions γ(a, x) and Γ(a, x). These are the integrals from 0 to x resp. from x to infinity. So they sum up to Γ(a).
Example: γ(5, 3) + Γ(5, 3) = 4,4337 + 19,5663 = 24 = Γ(5).
Since Γ(a) is the integral from 0 to infinity, and Γ(a, x) is the integral from x to infinity, it's clear that Γ(a, 0) = Γ(a).

3. The regularized Gamma functions P(a, x) and Q(a, x). These are simply γ(a, x) resp. Γ(a, x) divided by Γ(a). So they sum up to 1.
Example: P(5, 3) + Q(5, 3) = 4,4337/24 + 19,5663/24 = 0,18474 + 0,81526 = 1.

Dieter

[bshoring]

OK, I'm starting to see it. It appears to me the Henrici program gives us the Γ(a, x). Given (5,3), it yields 19,5663, which, if added to y(5,3), or 4,4337, gives us 24.

Then, would the above program be called the upper incomplete gamma?

Thanks,
Bob

(03-26-2015 11:10 PM)Dieter Wrote:  

(03-26-2015 10:41 PM)bshoring Wrote:  Do you think you can shed any light on this for me? I would just like to be able to understand what I am getting.

There are different Gamma functions.

1. The "normal" Gamma function Γ(a). For integer a this is equal to (a–1)!. So Γ(5) = 4! = 24.
This, let's say "complete" Gamma function Γ(a) can be expressed as an integral from 0 to infinity.

2. The incomplete Gamma functions γ(a, x) and Γ(a, x). These are the integrals from 0 to x resp. from x to infinity. So they sum up to Γ(a).
Example: γ(5, 3) + Γ(5, 3) = 4,4337 + 19,5663 = 24 = Γ(5).
Since Γ(a) is the integral from 0 to infinity, and Γ(a, x) is the integral from x to infinity, it's clear that Γ(a, 0) = Γ(a).

3. The regularized Gamma functions P(a, x) and Q(a, x). These are simply γ(a, x) resp. Γ(a, x) divided by Γ(a). So they sum up to 1.
Example: P(5, 3) + Q(5, 3) = 4,4337/24 + 19,5663/24 = 0,18474 + 0,81526 = 1.

Dieter

[Dieter]

(03-27-2015 07:26 PM)bshoring Wrote:  OK, I'm starting to see it. It appears to me the Henrici program gives us the Γ(a, x). Given (5,3), it yields 19,5663, which, if added to y(5,3), or 4,4337, gives us 24.

Yes, it seems to return Γ(a, x) if the arguments are entered a [ENTER] x.

(03-27-2015 07:26 PM)bshoring Wrote:  Then, would the above program be called the upper incomplete gamma?

That's how you may call it.

Dieter