Gamma Function by Stieltjes Continued Fraction

[lcwright1964]

On his excellent page on Gamma/Factorial approximations, Peter Luschny recommends a simplification of the standard 3-term Stieltjes continued fraction approximation of the Gamma function that is well-suited to calculator programming. Not only does the simplification render all constants as single digits, it just so happens to improve the performance by roughly an extra digit of accuracy over the original version. (Peter doesn't discuss this added benefit--I discovered it when examining the approximation myself in high precision environments on Free42 and in Mathematica.)

The structure of the program is inspired largely by JM Baillard's high quality Stirling series approximation, GAM+. I have in particular adopted Jean-Michel's treatment of how the leading factors get treated (i.e., what gets exponentiated and what is multiplied straight through), and his device of splitting one term into two factors so as to avoid overflow when the argument exceeds roughly 56. To save steps I have eliminated his treatment of integer arguments with the calculator's built-in FACT function.

The work is all done on the stack, except for the synthetic register M, which is used for the accumulation of the Pochhammer shift-and-divide product when the argument is 7 or less. If the program is entered from the keyboard using any regular register for this is fine.

In GAM+, Jean-Michel cleverly uses storage register arithmetic with stack registers to save steps and minimize stack movement--e.g. ST+ X accomplishes the same as ENTER + or 2 *, and ST- Y alters the contents of the Y register while leaving the X register unchanged. I have not used these devices, given that I wish to do a straight port to the HP67/97 forthwith and on those machines that is not possible. If the user wishes to use these concise commands the program can be shortened by 4 or 5 steps at least.

In high precision environments where adequate guard digits are not an issue (e.g., quadruple precision on Free42 Decimal), the approximation gives 10 digits easily (a minimum of about 10.4 exact decimal digits [edd] to be precise), and usually 11 or more across the argument range 0 to 70. Of course, rounding error in the 10-digit HP41 environment compromises this. I am finding in general that the program gives at least 8 good digits, usually 9 (with the 10th off by at most a couple ULP) and occasionally all 10. In my experience this is about as good as GAM+, only mine is smaller and seems to run at least a little faster.

Without further ado, here is the listing:

Code:


 01 LBL "GMST"
 02 7
 03 STO M
 04 ST/ M
 05 X<>Y
 06 LBL 01
 07 X>Y?
 08 GTO 02
 09 ST* M
 10 1
 11 +
 12 GTO 01
 13 LBL 02
 14 ENTER
 15 ENTER
 16 ENTER
 17 4
 18 *
 19 1/X
 20 +
 21 5
 22 *
 23 1/X
 24 X<>Y
 25 6
 26 *
 27 +
 28 ENTER
 29 +
 30 1/X
 31 X<>Y
 32 -
 33 E^X
 34 X<>Y
 35 R^
 36 0.5
 37 -
 38 LASTX
 39 *
 40 Y^X
 41 ENTER
 42 RDN
 43 *
 44 R^
 45 *
 46 PI
 47 ENTER
 48 +
 49 SQRT
 50 *
 51 0
 52 X<> M
 53 /
 54 END

Execution is self-evident--put the argument in the X register and XEQ [alpha]GMST[alpha].

I plan to port this directly and forthwith to the 67/97. I am keen on suggestions to save steps, because I would like to be able to fit the routine into the 49 steps of the HP33C/E or the HP25. And if I ever learn M-CODE programming this is one I will try to adapt as I expect it will be fast and give 10 full digits when 13-digit internal precision is used.

RAW and LIF files are attached in the ZIP file (I can't seem to upload them unzipped). Barcode PDF also attached.

Looking forward to feedback.

Les