HP Articles Forum

[Return to the Index ]
[ Previous | Next ]

GAMMA and GAMMLN: Fast, Accurate Lanczos Approximation for Gamma Functions

Posted by Steven Thomas Smith on 1 Sept 2009, 12:49 p.m.

This code computes the Gamma and log-Gamma functions, Gamma(x) and log|Gamma(x)|, over all x with full or nearly full 10 digit accuracy. It uses Lanczos's appoximation, which is amazing: very high accuracy Gamma(x) with a couple of LN calls and a handful of arithmetic operations. The Lanczos coefficients for (gamma = 5, N = 6) from Numerical Recipes in C, 2d ed. are used. The Lanczos coefficients are stored 3 at a time in the synthetic MNO status registers, and simple flag setting is used for the fast two-pass loop. The excellent SANDMATH'-5 packages's compiled M-code is only about 40% faster. (This may suggest that SANDMATH might consider using a fast implementation of Lanczos's approximation, which would provide both compiled Gamma(x) and log|Gamma(x)| within SANDMATH.) Lanczos effectively turns GAMMA and GAMMLN into fast system calls like LN and SIN.

Test case:

.5 XEQ GAMMA            -> 1.772453851 == (-1/2)! == SQRT(PI)      [in 3.1 s on i41CX+'s ridiculously slow "Default Speed"]
.5 XEQ SANDMATH-5 GAMMA -> 1.772453850 ==  SQRT(PI) - 1e-9         [in 1.9 s]

The synthetic text below may be generated using the PPC ROM's "LB" load bytes function with these hex codes:

Line 57: FE 95 39 52 39 38 59 94 01 20 86 50 97 49 97
Line 58: F8 7F 91 23 17 39 57 20 00
Line 86: FE 02 40 14 09 82 40 01 98 65 05 32 03 30 01
Line 87: F8 7F 07 61 80 09 17 30 01
Line 93: F7 09 18 93 85 33 29 99

The explanation (see Synthetic Programming Made Easy, p. 40f.) for these byte codes are that the 6 Lanczos coefficients c₁..c₆ are stored in the HP41 as:

76.18009172947146      = +7.618009173 E1   = 07 61 80 09 17 30 01
-86.50532032941677     = -8.650532033 E1   = 98 65 05 32 03 30 01
24.01409824083091      = +2.401409824 E1   = 02 40 14 09 82 40 01
-1.231739572450155     = -1.231739572 E0   = 91 23 17 39 57 20 00
1.208650973866179e-3   = +1.208650974 E-3  = 01 20 86 50 97 49 97
-5.395239384953e-6     = -5.395239385 E-6  = 95 39 52 39 38 59 94

This code is generated by Antonio Lagana's i41CX+ iPhone emulator, with editing for the synthetic text lines. (Antonio is upgrading i41CX+ in v. 3.2 to do synthetic text viewing automatically.)

01 LBL "GAMMA"
02 CF 05
03 CF 07
04 X<0?
05 SF 07
06 ENTER
07 INT
08 2
09 /
10 FRC
11 X!=0?
12 CF 07
13 RDN
14 ENTER
15 FRC
16 CF 06
17 X=0?
18 SF 06
19 RDN
20 FC?C 06
21 GTO 01
22 E
23 -
24 FACT
25 RTN
26 LBL "GAMMLN"
27 SF 05
28 LBL 01
29 CLA
30 E
31 X<>Y
32 GTO 03
33 LBL 02
34 ENTER
35 ABS
36 LN
37 ST- N
38 RDN
39 E
40 +
41 LBL 03
42 X<Y?
43 GTO 02
44 ENTER
45 ENTER
46 4.5
47 +
48 ST- N
49 LN
50 RCL Y
51 .5
52 -
53 *
54 ST+ N
55 CLX
56 RCL N
57 "\95\39\52\39\38\59\94\01\20\86\50\97\49\97"
58 >"\91\23\17\39\57\20\00"  ; hex, -5.395239385 E-6->O, +1.208650974 E-3->N, -1.231739572 E0-M
59 .
60 RCL Z
61 5
62 +
63 SF 06
64 LBL 04
65 X<> O
66 RCL O
67 /
68 +
69 E
70 ST- O
71 X<> N
72 RCL O
73 /
74 +
75 RCL N
76 ST- O
77 X<> M
78 RCL O
79 /
80 +
81 FC? 06
82 GTO 05
83 RCL O
84 RCL M
85 -
86 "\95\39\52\39\38\59\94\01\20\86\50\97\49\97"
87 >"07\61\80\09\17\30\01"  ; hex, +2.401409824 E1->O, -8.650532033 E1->N, +7.618009173 E1->M
88 FS?C 06
89 GTO 04
90 LBL 05
91 LN1+X
92 +
93 "09\18\93\85\33\29\99"  ; hex, LN(SQRT(2*PI)) = 9.189385332 E-1->M
94 RCL M
95 +
96 FS?C 05
97 RTN
98 E^X
99 FS?C 07
100 CHS
101 END

Edited: 25 Sept 2009, 12:00 p.m.

[ Return to the Message Index ]

Go back to the main exhibit hall