Inside Zen Gamma

If you're a floating point designer, it's a good idea to keep a close eye on the mathematics community. Although we're still using lots of algorithms that were discovered centuries ago, that doesn't mean that there aren't any new developments.

If you don't know what a Gamma function is, I can't blame you. It's not something you're using every day, although it is based on a concept that we floating point designers are only too familiar with: factorials. A factorial uses the following notation:

5!

Which means:

5! = 1 * 2 * 3 * 4 * 5

A factorial can also be expressed as a "Pi" function. You may have heard of the "Sigma" function, which means: the sum of "everything that follows it". The "Pi" function means: the product of "everything that follows it". So this is equivalent:

     __
n! = ||(n)

But a factorial only deals with natural numbers, where the Gamma function is also able to handle real numbers - and even complex numbers. The relation between the factorial and the Gamma function is expressed in this identity:

n! = Gamma (n+1)

But where calculating a factorial is trivial, calculating the Gamma function is whole different ballgame. Many are still using Stirling's approximation which was discovered in the 18th century.

But the story doesn't end there. In 1964 Cornelius Lanczos published a significantly better approximation, but it was much harder to use. Still, it gained some popularity when it was used in the popular "Numerical Recipes in C" - BTW, a volume which FOSS developers should avoid like the plague, because it is entirely comprised of non-free software.

In 2009 Cristinel Mortici wrote an article titled "A substantial improvement of the Stirling formula" in which he presented an equally simple formula which was more accurate, delivering about 7 good digits. But no matter how good the new formula is, it has the same flaw as the original equation: it becomes inaccurate close to zero.

There is a way around that. Simply use this identity:

         G (x + n)
G (x) =  ---------
         (x*(x+1)*(x+2)*..(x+n-1))

In laymans terms: if we calculate "Gamma(10)" instead of "Gamma(2)", we're getting into reasonably safe territory. Since Zen float has a precision of around 9 digits, that ain't half bad. We only have to apply Mortici's approximation and we're home free:

n! ~ sqrt(2*pi*n)(n/e + 1/(12*e*n))^n

So, now for the results. You'll see that the relative error is slowly decreasing, but even for "0.5", the result is usable if you're not too picky:

x	Calc(y)	Real(y)	Abs. Error	Rel. Error
0.5	1.7724507	1.7724538509	3.15E-06	1.78E-06
1	0.999998524	1	1.48E-06	1.48E-06
1.5	0.88622582	0.8862269255	1.11E-06	1.25E-06
2	0.999998908	1	1.09E-06	1.09E-06
2.5	1.32933903	1.3293403882	1.36E-06	1.02E-06
3	1.99999824	2	1.76E-06	8.80E-07
3.5	3.32334831	3.3233509704	2.66E-06	8.01E-07
4	5.99999559	6	4.41E-06	7.35E-07
4.5	11.6317201	11.6317283966	8.30E-06	7.13E-07
5	23.9999845	24	1.55E-05	6.46E-07
5.5	52.3427478	52.3427777846	3.00E-05	5.73E-07
6	119.9999381	120	6.19E-05	5.16E-07
6.5	287.885118	287.885277815	1.60E-04	5.55E-07
7	719.999568	720	4.32E-04	6.00E-07
7.5	1871.25312	1871.2543057978	1.19E-03	6.34E-07
8	5039.99719	5040	2.81E-03	5.58E-07
8.5	14034.39896	14034.4072934834	8.33E-03	5.94E-07
9	40319.9786	40320	2.14E-02	5.31E-07
9.5	119292.3913	119292.461994609	7.07E-02	5.93E-07
10	362879.79	362880	2.10E-01	5.79E-07

Mortici continues to make improvements, so may be in the future we'll use another approximation of his hand. Just wait and see!