How to approximate Gamma Function when inputting decimals?

Question

I was writing some code on approximating the gamma function using Stirling's method, but as far as I know it is only applicable for integer arguments. Does anyone know of another that can take decimals as arguments? Thanks

Answer 1

It is simpler and clearer if you use Stirling approximation as $$\log (\Gamma (x))=x (\log (x)-1)+\frac{1}{2} \left(\log \left(\frac{1}{x}\right)+\log (2 \pi )\right)+\frac{1}{12 x}-\frac{1}{360 x^3}+O\left(\frac{1}{x^5}\right)$$ and then $$\Gamma (x)=e^{\log (\Gamma (x)) }$$

You can use any real number.

Answer 2

One can derive many variants of the Stirling approximation. One variant that I have found useful for practical computation is $$\Gamma(x) \approx \sqrt{2\pi}\exp\left(\left(x-\frac{1}{2}\right)\ln x-x\right)\left(1+\frac{1}{12x}+\frac{1}{288x^2}-\ldots\right)$$

The Digital Library of Mathematical Functions lists this in §5.11(i) Poincaré-Type Expansions, provides additional terms, and includes pointers to publications by Wrench and Spira that give even more terms.

For numerical computation it is important to remember that $\exp$ magnifies errors in its argument. Achieving fully accurate results for any given target precision therefore requires computing the argument to $\exp$ and in particular $\ln x$ to more than the target precision. Also, the approximation above will not deliver very accurate results for arguments near zero. An approximation $\frac{1}{xp(x)}$ could be used for that, where $p$ is a polynomial: a series approximation of $\frac{1}{\Gamma(x)}$ around zero. For negative arguments the well-known reflection formula can be used.

When creating an implementation for a specific floating-point format, one typically replaces series expansions with the related minimax approximations. Standard tools like Maple and Mathematica have built-in facilities to generate these, or one can use the open-source Sollya tool, in particular its fpminimax command.

I have recently learned that there is a formalization of the above approach due to Nico Temme. The first appearance of this I can find is in N. M. Temme, "A set of algorithms for the incomplete gamma functions," Probability in the Engineering and Informational Sciences, Vol. 8, No. 2, April 1994, pp. 291-307, where $\Gamma(x) = \sqrt{2\pi}e^{-x}x^{x-\frac{1}{2}}\Gamma^{\star}(x), x\gt0$, and where $\Gamma^{\star}$ is called the "tempered gamma function". Next it appeared in Amparo Gil, Javier Segura, and Nico M. Temme, Numerical Methods for Special Functions, SIAM 2007, p. 243 as $\Gamma^{\star}(a)=\sqrt{\frac{a}{2\pi}}e^{a}a^{-a}\Gamma(a), a > 0$. More recently the name of $\Gamma^{\star}$ was changed to "regulated gamma function" in Amparo Gil, Javier Segura, and Nico M. Temme, "GammaCHI: a package for the inversion and computation of the gamma and chi-square cumulative distribution functions (central and noncentral)," Computer Physics Communications, Vol. 191, June 2015, pp. 132-139. In the authors' GNSTLIB library this function is included as gammastar.

Answer 3

The asymptotic expansions of the gamma function work for real (and complex) variable as well. For each $N\geq 1$, let us write $$ \Gamma (z) = \sqrt {2\pi } z^{z - 1/2} \mathrm{e}^{ - z} \left( {\sum\limits_{n = 0}^{N - 1} {( - 1)^n \frac{{\gamma _n }}{{z^n }}} + R_N (z)} \right) $$ and $$ \log \Gamma (z) = \tfrac{1}{2}\log (2\pi ) + \left( {z - \tfrac{1}{2}} \right)\log z - z + \sum\limits_{n = 1}^{N - 1} {\frac{{B_{2n} }}{{2n(2n - 1)z^{2n - 1} }}} + \widetilde R_N (z). $$ Here $\gamma_n$ denote the so-called Stirling coefficients and $B_n$ denote the Bernoulli numbers. The former may be computed via the recurrence relations $$ \gamma _{2n - 1} = - \frac{1}{{2n - 1}}\sum\limits_{k = 1}^n {\frac{{B_{2k} }}{{2k}}\gamma _{2n - 2k} } ,\quad \gamma _{2n} = - \frac{1}{{2n}}\sum\limits_{k = 1}^n {\frac{{B_{2k} }}{{2k}}\gamma _{2n - 2k + 1} } $$ for $n\geq 1$ with $\gamma_0=1$ (see this paper). The remainder terms satisfy the bounds $$ \left |R_N(z) \right | \le \left(\frac{\left |\gamma_N \right |}{|z|^N} + \frac{\left |\gamma_{N+1} \right |}{|z|^{N+1}}\right)\times \begin{cases} 1 & \text{ if }\; |\arg z| \leq \frac{\pi}{4}, \\ |\csc(2\arg z)| & \text{ if }\; \frac{\pi}{4} < |\arg z| < \frac{\pi}{2}, \end{cases} $$ and $$ |\widetilde R_N (z)| \le \frac{\left| B_{2N} \right|}{2N(2N - 1)\left| z \right|^{2N - 1} } \\ \times \begin{cases} 1 & \text{ if } \; \left|\arg z\right| \leq \tfrac{\pi}{4}, \\ \min \!\big(\left| \csc (2\arg z) \right|,\tfrac{1}{2}\chi (2N - 1) + 1\big) & \text{ if } \; \tfrac{\pi}{4} < \left|\arg z\right| \leq \tfrac{\pi}{2}, \\ \cfrac{\sqrt {2\pi (2N - 1)} }{2\left| \sin (\arg z) \right|^{2N - 1} } + \tfrac{1}{2}\chi (2N - 1) + 1 & \text{ if } \; \tfrac{\pi}{2} < \left|\arg z\right| < \pi. \end{cases} $$ Here $\chi(p)=\sqrt \pi \Gamma \big( \frac{p}{2} + 1 \big)/\Gamma \big( \frac{p}{2} + \frac{1}{2} \big)<\sqrt{\frac{\pi(p+1)}{2}}$ with any $p>0$. See here and here. We also have $$ \left| {R_N (z)} \right| \le \left( {\sqrt N + \tfrac{1}{2}} \right)\frac{{(1 + \zeta (N))\Gamma (N)}}{{(2\pi )^{N + 1} \left| z \right|^N }} $$ for $N\geq 2$ and $|\arg z|\leq \frac{\pi}{2}$, with $\zeta$ being the Riemann zeta function. Again, see here.

Other approximations include the Lánczos approximation and the Spouge approximation.