Gamma Function Using Spouge's Method

[Dieter]

Les, first of all thank you very much for sharing your insights regarding the Gamma function. There is a lot of interesting information in your post, and yes: the most crucial point is the limited working precision of the calculators we use.

(08-24-2015 09:23 AM)lcwright1964 Wrote:  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.

Right, that's the basic idea. I have been using such approximations (mostly rational ones) for other functions where the domain is limited by the calculator's working range. Take for instance the Normal quantile function. There is the well-known Hastings approximation which seems work for any argment greater than zero. But since the '41 only works down to 1E–99, the same (2;3) rational approximation may be used with a different set of coefficients optimized for exactly this working range. Et voilà: the resulting error is much lower and we can squeeze out another valid digit.

(08-24-2015 09:23 AM)lcwright1964 Wrote:  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.

OK, here we go.

First of all, my approximation determines the factorial instead of Gamma, so the argument is off by one. But the basic idea is the same for both cases. Simply take a look at the Lanczos approximation:

x! = [a0 + a1/(x+1) + a2/(x+2) + a3/(x+3) + a4/(x+4)] * ((x+1+c)/e)^(x+1-0,5)

There are six unknowns: a0...a4 and c, the constant in the second factor. This non-linear twist makes things a bit more complicated, since I have no tools for solving non-linear equation systems. So I simply decided to set c to the value you already metioned, i.e. c = 3,84088190863. This way the remaining coefficients can be determined quite easily – you can set up a simple linear equation system:

[a0 + a1 * 1/(x+1) + a2 * 1/(x+2) + a3 * 1/(x+3) + a4 * 1/(x+4)] = x! / ((x+4,84088190863)/e)^(x+0,5)

Do this for five different x, spread over the range 0...70, and you get a linear system with five equations in five unknowns. Which can be solved easily and gives five coefficients a0...a4.

Now plot an error graph and see where the error has its maxima and minima. Adjust your five nodes so that the error peaks even out and start over again. After some iterations the nodes were at 0, near 1.54 and 5.43, at 16.6 and 53. This way I got the following 12-digit result:

Code:

a0 =  3,2648 92414 09 E-02
a1 =  1,1886 88980 04
a2 = -1,0684 10297 87
a3 =  1,8871 19655 59 E-01
a4 = -2,7450 07080 40 E-03

With c = 3,84088190863 the relative error for 0!...70! oscillates nicely between plus and minus 3,9 E–11. Which means 10,4 edd.

(08-24-2015 09:23 AM)lcwright1964 Wrote:  If curious you could try your optimization method on my n = 5 result, reproduced here with 16-digit coefficients:

I did all this in a simple Excel spreadsheet with standard 15-digit accuracy. My Gamma/Factorial implementation (the reference against which the error is calculated) has 12-13 valid digits, so all this is a bit of working at the limit. That's why I do not trust my results for the n=5 case, although the resulting coefficients are very close to the ones you posted. But at least I can say that another term will reduce the relative error in 0...70 to not much more than about ±10^–12.

(08-24-2015 09:23 AM)lcwright1964 Wrote:  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.

What I can say is, that n=4 and 12-digit coefficients allow for a relative error less than ±4 E–11 (if evaluated exactly), so the 10-digit result should exact or maybe ±1 ULP off. Which seems to be comparable to the current Sandmath GAMMA implementation. One more term would probably yield 11 digits, or even twelve (if evaluated with sufficient precision).

(08-24-2015 09:23 AM)lcwright1964 Wrote:  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.

Well, it's not magic. Just simple (really simple) math. ;-)

Dieter