Gamma function has room for improvement

[verybadcat]

The gamma function produces wrong values for integer inputs. For example, gamma(4) returns not 6, but 5.99999999999999, and gamma(5) returns not 24, but 24.00000000000011. Similar errors happen for other integers. I realize that some error is inevitable as doubles actually have discrete values. But there are still two problems:

1. For integer inputs, it really should return the correct integer output.

2. The error in these particular cases are respectively 14x and 32x the granularity of double near their output values. That seems like more than it would need to be.

I'm new to this project but may use it heavily. Therefore, I'm willing to volunteer myself for a fix. But it seems like the kind of thing that should be discussed a little first.

[cdrnet]

For such small positive integers we can actually leverage _factorialCache which contains precomputed factorials for integers up to 170. Although we should be careful to keep it continuous.

But if there is a simple way to tune the Gamma function to full precision without sacrificing to much performance, yes, we should consider this as well.

[verybadcat]

Simple -- probably. Easy -- probably not. For example, opensource.apple.com has an implementation that looks quite simple, but it appears to have gone through several generations of tweaking, and it says copyright Apple, all rights reserved. http://www.opensource.apple.com/source/Libm/Libm-315/Source/ARM/tgamma.c

It's not much good to just correct the integer values. It should also behave smoothly, have the correct derivative, and so on.

[verybadcat]

It seems that Apple took their algorithm from netlib originally, and it has since been used by people all over the world. A Google search on the number 6.64561438202405440627855 is enough to verify that. In light of that, I think we are safe in using it. I see that John D. Cook has an implementation in C# which he says is public domain. I'm going to run some accuracy tests on it.

http://www.netlib.org/specfun/gamma http://www.johndcook.com/Gamma.cs

[cdrnet]

Great, thanks!

[verybadcat]

OK, it looks to me like it would be worth using. I made some changes to the public domain code, which appears to have been written when the precision of floating-point types is less than it is today. My changes are here, and you are more than welcome to use them if you want to.

https://github.com/verybadcat/mathnet-numerics

[evolvedmicrobe]

So there seem to be a lot of gammas out there, did you just copy this one from the johndcook site essentially? That implementation doesn't seem to be widely used so I am a bit scared of it so to speak and it doesn't seem to be quite like the apple one (they seem to overlap in certain ways but not others).

Another option might be to use the R option, they seem to have a lot of correction coefficients in there code, but it is not quite the same. Note that the R code detects integer values less than 50 directly and just does a factorial calculation with integers for better precision in this case, and the current gamma function could also be modified for this.

[verybadcat]

I took the John D. Cook one, then tested it at a bunch of values where I thought there was a significant risk of being wrong. It was pretty far off for some of them, so I made a lot of changes, e.g. considerably shrinking the region near zero where a certain power series is used. He has three ranges; my version actually only uses two of them, though the code for the third one is still in there -- I didn't want to remove it completely myself, before anyone else had a chance to take a look.

By the end, I wasn't finding any places where it was far wrong. Unfortunately, I didn't keep records while testing, but I believe everything was within 50x the precision of double, and frequently much closer than that.

I have no opinion one way or the other about the R implementation. My main concerns are that it should produce values that are close to correct, are exact at integers, and smooth everywhere. For example, it would not be good to take an implementation that is wrong at the integers, then correct it at the integers only, as the resulting function would not be smooth.

[verybadcat]

I just downloaded Math.Net numerics, and asked it to calculate SpecialFunctions.Gamma(5). The correct answer is 24, but it returned 24.000000000000114.

Can this be fixed? (My proposed changes would fix this. I don't really care how we fix it; if there is a better way to fix it, I would be fine with that.)

[cdrnet]

I finally had a look at how the new algorithm performs numerically. It is indeed mostly 1-2 digits more accurate than the existing implementation and behaves very nicely around positive integers.

I noticed two potential problems so far:

• Accuracy drops up to 5 digits around 10e-5 (but is much better again at 10e-10).
• GammaLn no longer supports 0.0 and negative numbers (where Gamma is non-negative).

[verybadcat]

I've got a fix for the accuracy issue, by simply using more terms of the power series and moving the lower breakpoint. I'll post once I've had a look at the other issue.

[verybadcat]

I've just posted a version that addresses both of your issues. I'm returning NaN for GammaLn in cases where Gamma is negative. But it's unclear if that's the right choice; one could instead return the real part.

https://github.com/mathnet/mathnet-numerics/pull/204/commits