# Algorithm to compute Gamma function

The question is simple. I would like to implement the Gamma function in my calculator written in C; however, I have not been able to find an easy way to programmatically compute an approximation to arbitrary precision.

Is there a good algorithm to compute approximations of the Gamma function?

Thanks!

-
Do you want an algorithm for its complex domain or just for real numbers? – Gerben Jan 27 '11 at 15:39
complex would be better – houbysoft Jan 27 '11 at 15:47
For reals it in is the GNU C library (haven't checked if it is mandated by the standard). – vonbrand Mar 1 '13 at 16:07
@vonbrand: not to arbitrary precision, though. – houbysoft Jan 8 '14 at 10:19

Looks like the Lanczos approximation will suit my needs : http://en.wikipedia.org/wiki/Lanczos_approximation

-
This note should be of interest. – J. M. Jul 27 '11 at 9:49
-
Down thumb for doing nothing but referring people to the Wikipedia. – CogitoErgoCogitoSum Jun 11 '15 at 1:03

Someone asked a similar question yesterday. I thought of replacing $e^{-t}$ by a series. $$\Gamma (z) = \int_{0}^{\infty} t^{z-1} e^{-t} dt \approx \sum_{j=0}^{a} \frac{(-1)^j b^{j+z}}{(j + z) j !} . \text{Choose } a > b ,$$ but as J. M. points out, I should have checked this a bit better. Take great care in the choice of $a, b$.

-
Have you tried implementing this? :) – J. M. Oct 9 '11 at 10:17
Not with more than one $z$. It can be very inaccurate depending on the choice of $a, b$. Now that I've tried a few more $z$, it looks like $a = b z^2$ is about right, but $b$ must also be chosen not too far from $z$. Example: $z = 6.6$, $b=2 \times z$, $a = b z^2$. This is only a guess! – Samuel Hambleton Oct 9 '11 at 11:32

It is now part of the C++11 standard library.

http://en.cppreference.com/w/cpp/numeric/math/tgamma

-
What you link to is the C++ interface to the standard C library. – vonbrand Mar 1 '13 at 16:09
This is not to arbitrary precision, though. – houbysoft Jan 8 '14 at 10:19

Try Nemes' approximation:

$$\ln ( \Gamma( x ) ) = \frac12 \ln( 2 \pi ) + \left( x - \frac12 \right) \ln( x ) - x + \frac x2 \ln\left( x \sinh\left( \frac1x \right) + \frac 1 {810 x^6} \right)$$

The last term, $\dfrac1 { 810 x^6}$ is an error-checking term and may be eliminated from your calculations.

-
Welcome to Maths.SE! I've edited your question using MathJax, our maths renderer. Check the source (by clicking edit) to see how it works. For further information about writing maths at this site see e.g. here, here, here and here. – Lord_Farin Mar 6 '15 at 21:25
Error-checking? No. Its a term used to increase precision. It has nothing to do with error checking. This isnt a binary transmission or whatever youre thinking. – CogitoErgoCogitoSum Jun 11 '15 at 1:04
For small values, this is only a very rough approximation. E.g., for x=1, the LHS is 0, but the RHS is 0.00015268924587. The error relative to scipy.special.gammaln is roughly proportional to x^{-4} for x in [1,3]. On the other hand, it's about four times faster. – Alex Coventry Apr 5 at 14:53

How about interpolation of the numbers in a look-up table featuring numbers drawn from a graph of the curve you are looking for? If you have the program Stella, you can literally shape by hand a curve that you desire, and it produces the numbers to be added to the table.

-
The domain of the Gamma function is the entire complex plane. You cannot use a lookup table for functions like this. – user11393 Jan 28 '13 at 9:37