GAMMA and GAMMLN: Fast, Accurate Lanczos Approximation

[eigen]

This code computes the Gamma and log-Gamma functions, Γ(x) and log|Γ(x)|, over all x with full or nearly full 10 digit accuracy. It uses Lanczos's appoximation, which is amazing: very high accuracy Γ(x) with a couple of LN calls and a handful of arithmetic operations. The Lanczos coefficients for (ϒ = 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 Γ(x) and log|Γ(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:

Code: Select all

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 c1..c6 are stored in the HP41 as:

Code: Select all

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.)

Code: Select all

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

[floppy_stuttgart]

Is the .raw file of this gamma function anywhere available? (I will upload in "V41 R9F" then transfer to my HP41 on my desk).

Furthermore, is there anywhere an inverse of the gamma function in ROMs/FOCAL? (area between x 0.00 and 1.461632144968.. which second number is the minimum in positive x)
https://mathoverflow.net/questions/1282 ... a-function
https://math.stackexchange.com/question ... an-inverse
https://en.wikipedia.org/wiki/Gamma_function
If not, I will do it

Update 27Aug2021 :
- a "gamma" was found in sandmath44
- "reverse gamma" done with the help of "solving" from math module (in sandmath44 too, however few issues with it)