Gamma Function by Stieltjes Continued Fraction

09012015, 10:31 PM
(This post was last modified: 09012015 10:40 PM by lcwright1964.)
Post: #1




Gamma Function by Stieltjes Continued Fraction
On his excellent page on Gamma/Factorial approximations, Peter Luschny recommends a simplification of the standard 3term Stieltjes continued fraction approximation of the Gamma function that is wellsuited 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 benefitI 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 JeanMichel'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 builtin 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 shiftanddivide 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+, JeanMichel cleverly uses storage register arithmetic with stack registers to save steps and minimize stack movemente.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 10digit 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:
Execution is selfevidentput 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 MCODE programming this is one I will try to adapt as I expect it will be fast and give 10 full digits when 13digit 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 

« Next Oldest  Next Newest »

User(s) browsing this thread: 1 Guest(s)