Approximation of the Gamma function for small value

Question

It is well known that $\Gamma(u) \stackrel{u \to 0}{\sim} \frac{1}{u}$.

I am looking for more precise information on the behavior of $\Gamma(x)$ when $x$ is small, ie: $x\to 0$.

My question is then, are there accurate (for small value) inequalities for the Gamma function ?

Any other information on the behavior is also very welcome.

Answer 1

Using Taylor series around $u=0$, you should get $$\Gamma(u)=\frac{1}{u}-\gamma +\frac{6 \gamma ^2+\pi ^2}{12} u+O\left(u^2\right)\tag 1$$

For $u=\frac 1 {10}$, this very limited expression would give $\approx 9.52169$ while the "exact" value would be $\approx 9.51351$.

The error is lower than $0.1$% for any $0\lt x\leq \frac 1 {10}$.

Edit

It is possible to slightly improve the above approximation builging the simplest Pade approximant of $u\, \Gamma(u)$ around $u=0$. This would lead to $$\Gamma(u)=\frac 1 u \times\frac{12\gamma+(\pi^2-6\gamma^2)u}{12\gamma+(\pi^2+6\gamma^2)u}\tag 2$$ Using $(1)$ leads to overestimates while using $(2)$ leads to underestimates which makes the average much better. So, a better approximation could be $$\Gamma(u)=\frac 1 {2u}\left(1-\gamma u+\frac{1}{12} \left(6 \gamma ^2+\pi ^2\right) u^2+\frac{12\gamma+(\pi^2-6\gamma^2)u}{12\gamma+(\pi^2+6\gamma^2)u} \right)\tag 3$$ for which the error is lower than $0.01$% for any $0\lt x\leq \frac 1 {10}$.

We could also consider $$\Gamma(u)=\alpha\left(\frac{1}{u}-\gamma +\frac{6 \gamma ^2+\pi ^2}{12} u \right)+(1-\alpha)\left(\frac 1 u \times\frac{12\gamma+(\pi^2-6\gamma^2)u}{12\gamma+(\pi^2+6\gamma^2)u}\right)$$ and optimize the $\alpha$ parameter. For the consider range $\alpha\approx 0.44$ seems to be quite good leading to a maximum error lower than $0.0002$% over that range.

If we consider $0< x \leq 1$, $\alpha\approx 0.35$ leads to errors smaller than $0.6$% for the entire range.

We could also show that the Padé approximant of $u\,\Gamma(u)$ lead to relative errors lower than $1$% for the range $-0.625 \leq x \leq 1.168 $.

Answer 2

$\Gamma(z)$ has a pole at zero, with Laurent series $$ \frac{1}{z}-\gamma+ \left( {\frac {{\pi }^{2}}{12}}+{\frac {{\gamma}^{2} }{2}} \right) z+ \left( -{\frac {\zeta \left( 3 \right) }{3}}-{\frac {{\pi }^{2}\gamma}{12}}-{\frac {{\gamma}^{3}}{6}} \right) {z}^{2}+ \left( {\frac {{\pi }^{4}}{160}}+{\frac {\zeta \left( 3 \right) \gamma}{3}}+{\frac {{\pi }^{2}{\gamma}^{2}}{24}}+{\frac {{\gamma}^{4} }{24}} \right) {z}^{3}+ \left( -{\frac {\zeta \left( 5 \right) }{5}}- {\frac {{\pi }^{4}\gamma}{160}}-{\frac {\zeta \left( 3 \right) {\pi } ^{2}}{36}}-{\frac {\zeta \left( 3 \right) {\gamma}^{2}}{6}}-{\frac {{ \pi }^{2}{\gamma}^{3}}{72}}-{\frac {{\gamma}^{5}}{120}} \right) {z}^{4 }+ \left( {\frac {61\,{\pi }^{6}}{120960}}+{\frac {\zeta \left( 5 \right) \gamma}{5}}+{\frac {{\pi }^{4}{\gamma}^{2}}{320}}+{\frac { \left( \zeta \left( 3 \right) \right) ^{2}}{18}}+{\frac {\zeta \left( 3 \right) {\pi }^{2}\gamma}{36}}+{\frac {\zeta \left( 3 \right) {\gamma}^{3}}{18}}+{\frac {{\pi }^{2}{\gamma}^{4}}{288}}+{ \frac {{\gamma}^{6}}{720}} \right) {z}^{5}+O \left( {z}^{6} \right) $$ Of course $\gamma$ is Euler's constant and $\zeta$ is Riemann's zeta function.

I wonder if it looks better if those $\pi^2$ and $\pi^4$ terms are written in terms of $\zeta(2)$ and $\zeta(4)$ and so on?

$$ \frac{1}{z}-\gamma+ \left( {\frac {\zeta \left( 2 \right) }{2}}+{\frac { {\gamma}^{2}}{2}} \right) z+ \left( -{\frac {\zeta \left( 3 \right) }{3}}-{\frac {\zeta \left( 2 \right) \gamma}{2}}-{\frac {{\gamma}^{3} }{6}} \right) {z}^{2}+ \left( {\frac {9\,\zeta \left( 4 \right) }{16}} +{\frac {\zeta \left( 3 \right) \gamma}{3}}+{\frac {\zeta \left( 2 \right) {\gamma}^{2}}{4}}+{\frac {{\gamma}^{4}}{24}} \right) {z}^{3}+ \left( -{\frac {\zeta \left( 5 \right) }{5}}-{\frac {9\,\zeta \left( 4 \right) \gamma}{16}}-{\frac {\zeta \left( 3 \right) \zeta \left( 2 \right) }{6}}-{\frac {\zeta \left( 3 \right) {\gamma}^{2}}{ 6}}-{\frac {\zeta \left( 2 \right) {\gamma}^{3}}{12}}-{\frac {{\gamma} ^{5}}{120}} \right) {z}^{4}+ \left( {\frac {61\,\zeta \left( 6 \right) }{128}}+{\frac {\zeta \left( 5 \right) \gamma}{5}}+{\frac {9 \,\zeta \left( 4 \right) {\gamma}^{2}}{32}}+{\frac { \left( \zeta \left( 3 \right) \right) ^{2}}{18}}+{\frac {\zeta \left( 3 \right) \zeta \left( 2 \right) \gamma}{6}}+{\frac {\zeta \left( 3 \right) { \gamma}^{3}}{18}}+{\frac {\zeta \left( 2 \right) {\gamma}^{4}}{48}}+{ \frac {{\gamma}^{6}}{720}} \right) {z}^{5}+O \left( {z}^{6} \right) $$

Maybe part of the $\zeta(4)$ should be a $\zeta(2)^2$ and similarly for the $\zeta(6)$ to make the pattern recognizable?

Answer 3

I'm a little late to the party but I derived something for this exact purpose.

Published here: http://albert.life/Math/FactorialApproximation/

Wolfram Alpha compatible: x!≈(0.44265*x+e^-eulergamma)^(0.986261*x) for -1/e<x<0

For applying above to Gamma, shift one, via Γ(x+1)==x!

Answer 4

Laurent Series

Here is one way to show how the coefficients in GEdgar's answer were derived.

In this answer, it is shown that for $f(x)=\log(\Gamma(x))$, $f(1)=0$, $f'(1)=-\gamma$, and for $k\ge2$, $$ f^{(k)}(1)=(-1)^k(k-1)!\,\zeta(k)\tag1 $$ From which it can be recursively derived, using the Leibniz Rule, that $$ \begin{align} \Gamma'(1)&=f'(1)\Gamma(1)&&=-\gamma\tag{2a}\\[12pt] \Gamma''(1)&=f''(1)\Gamma(1)+f'(1)\Gamma'(1)&&=\frac{\pi^2}6+\gamma^2\tag{2b}\\ \Gamma'''(1)&=f'''(1)\Gamma(1)+2f''(1)\Gamma'(1)\\ &+f'(1)\Gamma''(1)&&=-2\zeta(3)-\frac{\pi^2}2\gamma-\gamma^3\tag{2c}\\ \Gamma^{(4)}(1)&=f^{(4)}\Gamma(1)+3f'''(1)\Gamma'(1)\\ &+3f''(1)\Gamma''(1)+f'(1)\Gamma'''(1)&&=8\gamma\zeta(3)+\frac{3\pi^4}{20}+\gamma^2\pi^2+\gamma^4\tag{2d}\\ &\ \,\vdots\\ \Gamma^{(n)}(1)&=-\gamma\,\Gamma^{(n-1)}(1)+\sum_{k=2}^n(-1)^k\frac{(n-1)!}{(n-k)!}\,\zeta(k)\,\Gamma^{(n-k)}(1)\hspace{-4cm}\tag{2e} \end{align} $$ Thus, Taylor's Theorem says, $$ \begin{align} \Gamma(1+x) &=1-\gamma x+\left(\frac{\pi^2}6+\gamma^2\right)\frac{x^2}2-\left(2\zeta(3)+\frac{\pi^2}2\gamma+\gamma^3\right)\frac{x^3}6\\ &+\left(8\gamma\zeta(3)+\frac{3\pi^4}{20}+\gamma^2\pi^2+\gamma^4\right)\frac{x^4}{24}+\dots\tag3 \end{align} $$ and therefore, since $\Gamma(x)=\frac1x\Gamma(1+x)$, $$ \begin{align} \Gamma(x) &=\frac1x-\gamma+\left(\frac{\pi^2}6+\gamma^2\right)\frac{x}2-\left(2\zeta(3)+\frac{\pi^2}2\gamma+\gamma^3\right)\frac{x^2}6\\ &+\left(8\gamma\zeta(3)+\frac{3\pi^4}{20}+\gamma^2\pi^2+\gamma^4\right)\frac{x^3}{24}+\dots\tag4 \end{align} $$ Here is a Mathematica implementation of $\text{(2e)}$ to compute $\Gamma^{(n)}(1)$:

G[0] = 1; G[n_] := G[n] = 
  Sum[(-1)^k (n-1)!/(n-k)! Zeta[k] G[n-k], {k, 2, n}] - EulerGamma G[n-1]

Simplify[G[4]] should give the value in $\text{(2d)}$.

---

A Decent Approximation

I have often used the slight modification of Stirling's formula $$ \Gamma(x)\approx\sqrt{2\pi\left(x+\tfrac16\right)}\,\frac{x^{x-1}}{e^x}\tag5 $$ which gives very good results, even for $x$ close to $0$. Although the Laurent expansion, $(4)$, gives a better approximation for $0\lt x\lt\frac13$, approximation $(5)$ stays within $2\frac13\%$ of the correct answer all the way to $0$.

Although we computed to the third order in $(4)$, the graph only contains terms to the second order.

For large $x$, the relative error of the approximation in $(5)$ is $\sim-\frac1{144x^2}$, whereas Stirling has a relative error of $\sim-\frac1{12x}$.