Gamma Function Using Spouge's Method

[Dieter]

(08-12-2015 05:39 PM)Dieter Wrote:  Addendum: This also means that a value as high as a=12,5 does not make much sense here. Try Namir's program with a=7 (replace 12.5 with 7 and 1.012 with 1.006). My tests returned results at least as accurate as with a=12,5 (mostly even better) while the program runs about twice as fast.

FTR, here is my implementation for the '41, based on a "hard-wired" a=7.
Like the factorial function on various other HPs, it actually calculates Γ(x+1).

Code:

01 LBL "GAM+1"
02 STO 00
03 6
04 +
05 STO 04
06 5,5
07 STO 02
08 INT
09 FACT
10 STO 03
11 SIGN
12 STO 01
13 CLx
14 LBL 01
15 RCL 01
16 RCL 02
17 y^x
18 RCL 01
19 e^x
20 *
21 RCL 03
22 /
23 RCL 04
24 /
25 +
26 RCL 02
27 INT
28 X=0?
29 GTO 02
30 CHS
31 ST/ 03
32 SIGN
33 ST- 01
34 ST+ 02
35 ST+ 04
36 RDN
37 GTO 01
38 LBL 02
39 RDN
40 ABS
41 Pi
42 ST+ X
43 SQRT
44 +
45 RCL 00
46 7
47 +
48 STO 03
49 RCL 00
50 ,5
51 +
52 y^x
53 *
54 RCL 03
55 e^x
56 /
57 END

Since it doesn't seem to get much better than 8 valid digits in the final result, a=7 is sufficient here. Evaluated exactly, the largest relative error is approx. 2,3 E–9. Actually the error may be slightly higher. But it should be less than in the original a=12,5 version.

Some examples:

Code:

 3   XEQ"GAM+1"  =>  5,999999994
 5   XEQ"GAM+1"  =>  120,0000000
 7   XEQ"GAM+1"  =>  5040,000009
8,35 XEQ"GAM+1"  =>  85899,08333      (exact: ...363)
17,3 XEQ"GAM+1"  =>  8,416039696 E+14 (exact: ...738)
50,5 XEQ"GAM+1"  =>  2,166683759 E+65 (exact: ...771)

The current version throws an OUT OF RANGE error for x > 54,9. Maybe this can be improved somehow.

Edit: the following modification may be a solution. However, the resulting accuracy is slightly less than in the first version.

Code:

...
45 RCL 00
46 7
47 +
48 STO 03
49 LN
50 RCL 00
51 ,5
52 +
53 *
54 RCL 03
55 -
56 e^x
57 *
58 END

This way x=69,9575744 returns 9,999997034 E+99 (exact: ...7560).

Dieter