Gamma Function Using Spouge's Method

08242015, 07:32 PM
(This post was last modified: 08242015 07:39 PM by Dieter.)
Post: #24




RE: Gamma Function Using Spouge's Methjod
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.
(08242015 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 wellknown 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. (08242015 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+10,5) There are six unknowns: a0...a4 and c, the constant in the second factor. This nonlinear twist makes things a bit more complicated, since I have no tools for solving nonlinear 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 12digit result: Code: a0 = 3,2648 92414 09 E02 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. (08242015 09:23 AM)lcwright1964 Wrote: If curious you could try your optimization method on my n = 5 result, reproduced here with 16digit coefficients: I did all this in a simple Excel spreadsheet with standard 15digit accuracy. My Gamma/Factorial implementation (the reference against which the error is calculated) has 1213 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. (08242015 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. What I can say is, that n=4 and 12digit coefficients allow for a relative error less than ±4 E–11 (if evaluated exactly), so the 10digit 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). (08242015 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 

« Next Oldest  Next Newest »

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