Gamma Function Using Spouge's Method

[lcwright1964]

(08-23-2015 09:01 PM)Dieter Wrote:  Now maybe you or someone else may setup a "really" optimized approximation including the final 3,84... coefficient. ;-)

Not quite sure what you did to tweak things. I started with the 3.84... (which is actually 4.34 - .5) parameter as provided by Pugh's research--this is in an appendix to the thesis. The table values Pugh produced gives the values of the parameter that provides the best performance of the Lanczos formula with a given number of terms. However, Pugh examined not only real arguments, but allowable complex arguments too, and his reported max relative errors reflect this. So it could be that there is something with the same parameter and the same number of terms that does even better if one confines oneself to not only the positive real line, but the very specific subdomain of about 0 < z < 70. I don't know what you did to optimize my function as presented, but I would like to know how you went about it and what you got.

The theory behind where Pugh got his optimized parameters is discussed in the thesis and much of it sails over my head. Like I said above, I started with the parameter for my desired number of terms, run things through a Maple worksheet that ports Viktor Toth's C++ implementation of Godfrey's matrix approach to computing the coefficients (Pugh references this in the thesis), compute everything to an arbitrary number of digits, and in this case just present 10 digits as that is what would be used in FOCAL.

If you are curious about more digits with this number of terms, here is 16 digit precision to match the precision of Pugh parameter in full:

(0.3264892412943410e-1 + 1.188688981251149/z - 1.068410309377202/(z+1) + .1887119901913100/(z+2) - 0.2745021708565992e-2/(z+3)) * ((z + 3.840881908632582)/exp(1))^(z-1/2)

If curious you could try your optimization method on my n = 5 result, reproduced here with 16-digit coefficients:

(0.9446965718573752e-2 + 1.367185253712206/z - 1.839535625334429/(z+1) + .6838535775172448/(z+2) - 0.6542905498247415e-1/(z+3) + 0.6572071053708592e-3/(z+4)) * ((z + 5.081000210222555)/exp(1))^(z-1/2)

Finally, Angel uses an n = 6 version as calculated by Lanczos himself, but rearranged by Viktor Toth. Here is my n = 6 version using Pugh's suggested parameter, again expressed with 16-digit coefficients for your optimizing pleasure:

(0.2849626875039468e-2 + 1.519901663012889/z - 2.706750200734360/(z+1) + 1.554955033608236/(z+2) - .3214751818926898/(z+3) + 0.1886890596091082e-1/(z+4) - 0.1321092292846819e-3/(z+5)) * ((z + 6.279505747515941)/exp(1))^(z-1/2)

Glen Pugh continued to refine his work in subsequent years, and in 2006 shared with me an extended and corrected table of his parameters, up to n = 100. It doesn't include all of the columns of the original thesis table, but it does have the relevant ones. It is attached.

I am keen to learn about your strategy of optimizing these formulae for use on a specific real domain. If you are squeezing out an extra digit of accuracy, maybe we could generate something that is shorter for Angel's MCODE yet does just as well. If you have managed to produce an improved n = 4 formula with 10-digit constants that gives 10+ edd in the 0 to 70 domain, then perhaps Angel might have something to consider that offers real improvement over the excellent results his existing code already gets with an n = 6 formula using an arbitrary parameter that we know from Pugh's research is not optimal.

Let me know how you are doing this. You really have me curious! If you can start by sending me your optimized version of, say, my n = 4 result, that would really be fun to play with in Mathematica, particularly if you can produce one with more digits based on what I have given above.

Les