6

I know there are several approximations of the Gamma function that provide decent approximations of this function.

I was wondering, how can I efficiently estimate specific values of the Gamma function, like $\Gamma (\frac{1}{3})$ or $\Gamma (\frac{1}{4})$, to a high degree of accuracy (unlike Stirling's approximation and other low-accuracy methods)?

J. M. ain't a mathematician
  • 75,051
Argon
  • 25,303
  • Possible duplicate of http://math.stackexchange.com/questions/19236/algorithm-to-compute-gamma-function – lhf Jul 10 '12 at 21:34
  • Wikipedia mentions a fast algorithm for algebraic arguments by Karatsuba. – lhf Jul 10 '12 at 21:37
  • Also related – draks ... Jul 10 '12 at 22:20
  • 1
    Don't knock Stirling. The thing is that it is intended to be used for large arguments (it is an asymptotic expansion, after all), and forcing it to work for small arguments is not unlike using a shoehorn to strain tea; it can be made to work, but it's very awkward. – J. M. ain't a mathematician Jul 11 '12 at 00:13
  • 1
    It seems to me someone recommended something like this once. To compute say $\Gamma(1/4)$ accurately, first use Stirling to compute accurately some large value like $\Gamma(100+1/4)$ or something, then apply the functional equation to get down to $\Gamma(1/4)$. – GEdgar Jul 11 '12 at 00:25
  • @J.M. How big is needed to get sufficient accuracy? From my computations (bear in mind I am experiencing some truncating errors and other things that lower accuracy), estimating $\Gamma (171.5)$ gives me the estimate $\Gamma (\frac{1}{2})=1.7715877615902382$ - 3 decimal places correct. This is the largest number that my program can handle without getting errors because the estimate of $\Gamma(x)$ is too large. So far, this method is not very efficient. Am I missing something? – Argon Jul 13 '12 at 21:05
  • Argon, like I said, large arguments only. For small arguments, you have to repeatedly use the functional equation like @GEdgar says (which beings me back to my analogy; it's rather awkward to program). How are you using Stirling, BTW? Remember that to use an asymptotic series computationally, one keeps summing for as long as the terms of the series are not growing; if the $(k+1)$-th term is suddenly bigger in magnitude than the $k$-th term, you cut off at the $k$-th term. – J. M. ain't a mathematician Jul 14 '12 at 00:02
  • And of course if your computer language is limited in the size of numbers, use the logarithm. – GEdgar Jul 14 '12 at 03:12
  • @J.M. I see. I used the "truncated" version, not the series. – Argon Jul 16 '12 at 02:00

4 Answers4

7

Gourdon and Sebah pages on "Numbers, constants and computation" are usually interesting for this kind of record (start by clicking on 'constants' at the left).

They reference Shigeru Kondo and Steve Pagliarulo for having computed 10^10 digits of $\Gamma\left(\frac 14\right)$ and $\Gamma\left(\frac 13\right)$ in 2010 and 2009 (see here). Many methods for high precision are clearly exposed in the other pages of the first link (FFT, binary splitting and this kind of things should be simply explained at least in the .ps or .pdf files).

This paper of Alexander Yee concerning high precision evaluation of $\pi$ contains too references to Kondo & Pagliarulo's record (he has a page of records too but without $\Gamma$ I fear).

In the special case $\Gamma\left(\frac k{24}\right)$ a quadratic algorithm is proposed in Weisstein's Gamma Function page (see around (99)).
This method is the famous AGM method and you may find too in page 6 of Sebah and Gourdon's postscript paper on Gamma (or here) that : $$\displaystyle \Gamma\left(\frac 14\right)^2=\frac{(2\pi)^{3/2}}{AGM(\sqrt{2},1)}$$

You may find in the excellent book of the Borweins "Pi and the AGM" following derivation using the complete elliptic integral of the first kind $\rm{K}$ : page25

Hoping it helped anyway,

Raymond Manzoni
  • 43,021
  • 5
  • 86
  • 140
  • 1
    The key point is that the gamma function for small rational arguments agrees with certain values of the complete elliptic integral of the first kind, and the AGM is a very efficient way to numerically evaluate the complete elliptic integral. I might write a few more details later... – J. M. ain't a mathematician Jul 11 '12 at 00:12
  • I'll just note that in the extract from the Borweins' book, they use the modulus $k$ as the argument of the complete elliptic integrals instead of the parameter $m$ that I used in my answer. The two conventions are related, though: $m=k^2$. Thus, the $K\left(\dfrac1{\sqrt2}\right)$ of Borwein corresponds to $K\left(\dfrac12\right)$ in my answer. – J. M. ain't a mathematician Jul 13 '12 at 19:28
  • @J.M.: I noticed that too : the online definitions of $\rm{K}$ all use $\sqrt{1-k^2t^2}$ in the integral these days (including Weisstein) while A&S and Alpha/Mathematica use $\sqrt{1-kt^2}$. So that the parameters should usually be squared before being proposed to Alpha (to make things even more confusing Alpha changes the simple parameter $p$ in $p^2$ but not if written $p^1$ :-)). – Raymond Manzoni Jul 13 '12 at 19:59
  • 1
    Right; Mathematica /Wolfram Alpha and Abramowitz and Stegun use the parameter (I point to A&S as my reason for preferring the parameter convention), while Maple, MATLAB, MathWorld and DLMF use the modulus convention. Moral lesson being, figure out what convention your author is using first before applying any formula containing elliptic integrals (and elliptic functions, too), and always mention what convention you're using in your writing. – J. M. ain't a mathematician Jul 13 '12 at 23:54
  • A link to Sebah and Gourdon's "Introduction to the Gamma Function". – Raymond Manzoni May 04 '14 at 15:13
6

(This is an addendum of sorts to Raymond's answer, and an expansion of my comment there.)

For certain small rational arguments, Borwein and Zucker show that the gamma function can be evaluated using the complete elliptic integral of the first kind $K(m)$ (where $m$ is the parameter), evaluated at so-called "singular values". Since the complete elliptic integral of the first kind can be evaluated using the quadratically convergent arithmetic-geometric mean,

$$K(m)=\frac{\pi}{2\mathrm{AGM}(1,\sqrt{1-m})}$$

we also have a fast way to evaluate the gamma function for those particular rational arguments.

Letting the singular values be denoted by $k_r=\lambda^\ast(r)$ (where $\lambda^\ast(r)$ is the elliptic lambda function), e.g.

$$\begin{align*} k_1&=\frac{1}{\sqrt{2}}\\ k_2&=\sqrt{2}-1\\ k_3&=\frac{\sqrt{3}-1}{2\sqrt{2}}\\ k_4&=3-2\sqrt{2}\\ k_5&=\sqrt{\frac{1}{2}-\sqrt{\sqrt{5}-2}}\\ k_6&=\left(2-\sqrt{3}\right)\left(\sqrt{3}-\sqrt{2}\right) \end{align*}$$

we have for instance the evaluations

$$\begin{align*} \Gamma\left(\frac13\right)&=\frac{2^{7/9}\pi^{1/3}}{3^{1/12}}(K(k_3^2))^{1/3}\\ \Gamma\left(\frac14\right)&=2\pi^{1/4}\sqrt{K(k_1^2)}\\ \Gamma\left(\frac16\right)&=\frac{\sqrt{3}}{\sqrt[3]{2}\sqrt{\pi}}\left(\Gamma\left(\frac13\right)\right)^2\\ \Gamma\left(\frac18\right)\Gamma\left(\frac38\right)&=\sqrt{\sqrt{2}-1}\sqrt{\pi}2^{13/4}K(k_2^2)\\ \frac{\Gamma\left(\frac18\right)}{\Gamma\left(\frac38\right)}&=\frac{2\sqrt{\sqrt{2}+1}}{\sqrt[4]{\pi}}\sqrt{K(k_1^2)} \end{align*}$$

See the Borwein/Zucker paper for more details.

With regards to the Lanczos approximation mentioned by Rory in his answer, Paul Godfrey presents a short note on how to derive the approximation that is somewhat more straightforward than the original derivation by Lanczos.

For other methods, see this paper by Schmelzer and Trefethen, and this webpage by Peter Luschny.

J. M. ain't a mathematician
  • 75,051
  • +1: Thanks for the detailed addition J.M.! The paper of Borwein&Zucker doesn't seem freely available so that I added a scan of Borwein's book in my answer for the case $\Gamma(1/4)$. Available where the "AGM revisited" from Borwein&all and another paper a little linked to the other elliptic integral thread of these days... here – Raymond Manzoni Jul 13 '12 at 19:11
  • Yeah, it's a bit inconvenient that there's no freely accessible version of that paper, so I tried to list at least the simpler results from that paper. If you're interested in seeing it too, I suppose I could send you a copy... – J. M. ain't a mathematician Jul 13 '12 at 19:14
  • Thanks J.M. but it won't be needed. The general idea is pretty clear from the book 'Pi and the AGM' I'm reading these days (very interesting if not cheap). Cheers, – Raymond Manzoni Jul 13 '12 at 19:21
3

I'm not sure what sort of approximation you need.

Have you looked at the Lanczos approximation? It provides the technique used in the classic book Numerical Recipes. It can be used to calculate reasonable approximations for $\Gamma$ (8 decimal places or so) even on devices with constrained resources, like the 12C financial calculator from HP.

Rory Molinari
  • 48
0

I made a graph on Desmos and a c++ implementation about this earlier today if anybody else wants to take a look at it. It uses the Lanczos Approximation. Here’s the graph and here’s the C++ implementation which works with complex numbers too

Anik Patel
  • 478