25
$\begingroup$

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!

  • 0
    @vonbrand: not to arbitrary precision, though.2014-01-08

6 Answers 6

14

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

Thanks for your help!

  • 1
    [This note](http://$m$y.fit.e$d$u/~gabdo/gamma.txt) should be of interest.2011-07-27
5

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

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

  • 0
    This is not to arbitrary precision, though.2014-01-08
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$

2

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

  • 0
    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!2011-10-09
2

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.

  • 0
    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.2016-04-05
1

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.

  • 1
    The domain of the Gamma function is the entire complex plane. You cannot use a lookup table for functions like this.2013-01-28