Algorithm to compute Gamma function

Question

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!

Answer 1

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

Thanks for your help!

Answer 2

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

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

Answer 3

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

Answer 4

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

Answer 5

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.

Here is my reference.

Answer 6

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.

Answer 7

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!