23
$\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
    Do you want an algorithm for its complex domain or just for real numbers?2011-01-27
  • 0
    complex would be better2011-01-27
  • 1
    For reals it in is the [GNU C library](http://www.delorie.com/gnu/docs/glibc/libc_393.html) (haven't checked if it is mandated by the standard).2013-03-01
  • 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://my.fit.edu/~gabdo/gamma.txt) should be of interest.2011-07-27
3

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

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

  • 0
    What you link to is the C++ interface to the standard C library.2013-03-01
  • 0
    This is not to arbitrary precision, though.2014-01-08
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
    Have you *tried* implementing this? :)2011-10-09
  • 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
    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](http://meta.math.stackexchange.com/questions/5020/), [here](http://meta.stackoverflow.com/a/70559/155238), [here](http://meta.math.stackexchange.com/questions/1773/) and [here](http://math.stackexchange.com/editing-help#latexhelp/notation).2015-03-06
  • 3
    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.2015-06-11
  • 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
2

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

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