# 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? Commented Jan 27, 2011 at 15:39
• complex would be better Commented Jan 27, 2011 at 15:47
• For reals it in is the GNU C library (haven't checked if it is mandated by the standard). Commented Mar 1, 2013 at 16:07
• @vonbrand: not to arbitrary precision, though. Commented Jan 8, 2014 at 10:19
• State of the art: Fredrik Johansson, "Arbitrary-precision computation of the gamma function." arXiv preprint arXiv:2109.08392 (2021). Commented Feb 20, 2022 at 21:33

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

• This note should be of interest. Commented Jul 27, 2011 at 9:49

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. Commented Mar 1, 2013 at 16:09
• This is not to arbitrary precision, though. Commented Jan 8, 2014 at 10:19

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$$

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? :) Commented Oct 9, 2011 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! Commented Oct 9, 2011 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.

• 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. Commented Jun 11, 2015 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. Commented Apr 5, 2016 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
Commented Jan 28, 2013 at 9:37

You can use Lanczos approximation to implement it. The C# implementation looks like:

using System;

class Program
{
static void Main(string[] args)
{
// Test cases for the SpecialFunctionGamma method
double[] testValues = { 0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 5.0, 6.0, 7.0 };
foreach (double value in testValues)
{
Console.WriteLine(\$"Gamma({value}) = {SpecialFunctionGamma(value)}");
}
}

//Lanczos Approximation for Gamma Function
private static double SpecialFunctionGamma(double z)
{
double[] p =
{
676.5203681218851,
-1259.1392167224028,
771.32342877765313,
-176.61502916214059,
12.507343278686905,
-0.13857109526572012,
9.9843695780195716e-6,
1.5056327351493116e-7
};
if (z < 0.5)
return Math.PI / (Math.Sin(Math.PI * z) * SpecialFunctionGamma(1 - z));
z -= 1;
double x = 0.99999999999980993;
for (int i = 0; i < p.Length; i++)
{
x += p[i] / (z + i + 1);
}

double t = z + p.Length - 0.5;
return Math.Sqrt(2 * Math.PI) * Math.Pow(t, z + 0.5) * Math.Exp(-t) * x;
}
}


Try it online!