# 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

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

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. is a poor mathematician 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

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 '16 at 14:53

Take the divergent series

$$\hat J_\nu(z):=\sqrt{\dfrac{2}{\pi z}}\left\{p(z)\cos\left(z-\dfrac{2\nu+1}{4}\pi\right)-q(z)\sin\left(z-\dfrac{2\nu+1}{4}\pi\right)\right\}$$

with

$$p(z):=1-\dfrac{(4\nu^2-1^2)(4\nu^2-3^2)}{2!(8z)^2}+\dfrac{(4\nu^2-1^2)(4\nu^2-3^2)(4\nu^2-5^2)(4\nu^2-7^2)}{4!(8z)^4}-+\dots$$

and

$$q(z):=\dfrac{(4\nu^2-1^2)}{1!(8z)^1}-\dfrac{(4\nu^2-1^2)(4\nu^2-3^2)(4\nu^2-5^2)}{3!(8z)^3}+-\dots$$

For the Bessel function $J_\nu(z)$ we get the estimation

$$\left\vert J_\nu(z)-\hat J_\nu^{(k)}(z)\right\vert\le O\left(z^{-k-3/2}\right)$$

for $\vert z\vert\ge1$ and $Re\{z\}\ge0$ with the partial sum $\hat J_\nu^{(k)}(z)$ whose last summand has the term

$$\dfrac{(4\nu^2-1^2)\cdots(4\nu^2-(2k-1)^2)}{k!(8z)^k}.$$

The Bessel function $J_\nu(z)$ can be calculated to

$$J_\nu(z)=\sum_{j=0}^\infty\dfrac{(-1)^j}{j!\Gamma(\nu+j+1)}\left(\dfrac{z}{2}\right)^{\nu+2j}.$$

With $\Gamma(x+1)=x\Gamma(x)$ we obtain

$$J_\nu(z)=\dfrac{1}{\Gamma(\nu+1)}\sum_{j=0}^\infty\dfrac{(-1)^j}{j!(\nu+1)\cdots(\nu+j)}\left(\dfrac{z}{2}\right)^{\nu+2j}.$$

Take

$$K_\nu(z)=\sum_{j=0}^\infty\dfrac{(-1)^j}{j!(\nu+1)\cdots(\nu+j)}\left(\dfrac{z}{2}\right)^{\nu+2j}.$$

Then we obtain

$$\lim_{n\to\infty}\dfrac{K_\nu\left(n\pi+\dfrac{2\nu+1}{4}\pi\right)}{\hat J_\nu^{(k)}\left(n\pi+\dfrac{2\nu+1}{4}\pi\right)}=\Gamma(\nu+1)$$

and the gamma function drops out of the game.

As an example we get for $\nu=0.5$, $k=10$ and $z=4\pi+\dfrac{\pi}{2}$

$$\dfrac{K_{0.5}\left(z\right)}{\hat J_{0.5}^{(10)}\left(z\right)}=0.88622692546003057$$

and

$$\dfrac{K_{0.5}\left(z+0.1\right)}{\hat J_{0.5}^{(10)}\left(z+0.1\right)}=0.88622692544073867.$$

But

$$\Gamma\left(\dfrac{3}{2}\right)=\dfrac{1}{2}\sqrt{\pi}= 0.88622692545275794.$$

Because the series $\hat J_{0.5}(z)$ terminates for $\nu=\dfrac{1}{2}$ we even have $\hat J_{0.5}(z)=J_{0.5}(z)$ and the difference is only due to the precision of our calculation that was made with the type $\mathbf{double}$ (or System.Double) in C#.

For $\nu=0.25$, $k=10$ and $z=4\pi+\dfrac{3\pi}{8}$ we get

$$\dfrac{K_{0.25}\left(z\right)}{\hat J_{0.25}^{(10)}\left(z\right)}=0.90640247708308364\approx\Gamma(1.25)=\dfrac{1}{4}\Gamma(0.25)$$

and for $\nu=0.75$, $k=10$ and $z=4\pi+\dfrac{5\pi}{8}$ we get

$$\dfrac{K_{0.75}\left(z\right)}{\hat J_{0.75}^{(10)}\left(z\right)}=0.91906252683450718\approx\Gamma(1.75)=\dfrac{3}{4}\Gamma(0.75).$$

With

$$\dfrac{\pi}{\sin\pi x}=\Gamma(x)\Gamma(1-x)$$

we obtain

$$(4\cdot 0.90640247708308364)\cdot(4\cdot 0.91906252683450718\,/\,3)-\dfrac{\pi}{\sin\pi/4}=0.000000000065\dots$$

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