Gamma Function by Stieltjes Continued Fraction
09-01-2015, 10:31 PM (This post was last modified: 09-01-2015 10:40 PM by lcwright1964.)
Post: #1
 lcwright1964 Member Posts: 66 Joined: Jan 2014
Gamma Function by Stieltjes Continued Fraction
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

Attached File(s)