Approximation of the Gamma function for small value 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.
 A: 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 $.
A: $\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?
A: 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}$.
A: 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!
