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 c1..c6 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.