Integral involving Gamma function: density of Kendall-Ressel family of distributions I came across the following function when reading the famous paper of Letac and Mora "Natural exponential families with cubic variance functions", i.e.,
$$f(x) = \frac{x^x e^{-x}}{\Gamma(x+2)}$$ for $x \ge 0$. Proposition 5.5 there showed via a transform of a Levy process that $f$ is a density on $(0,\infty)$ without the need to compute the integral
$$\int_{0}^{\infty} f(x)dx$$
In fact, $f$ is called Kendall-Ressel density.
Tonight, I used Mathematica to compute the above integral but Mathematica did not give the answer $1$. Instead, it was stuck and returned the orginal integral. So, I guess a direct computation is not trivial!?
I am curious on how to directly show that $$\int_{0}^{\infty} f(x)dx = 1$$ Any suggestions? Thank you.
Update on June 21, 2018: to follow up on Nemo's solution: is there a general formula or recursive formula for $$I(b,\alpha) = \int_{C} \frac{t^{-1}dt}{(t + \ln(-t)+b)^\alpha},$$ where $\alpha$ is a positive constant? Say, $\alpha$ is a natural number.
 A: According to Hankel's formula
$$
\frac{1}{\Gamma(z)}=\frac{i}{2\pi}\int_C(-t)^{-z}e^{-t}dt,
$$
where $C$ is Hankel contour. So
$$
\frac{x^x}{\Gamma(x+1)}=\frac{i}{2\pi}\int_C(-t)^{-x-1}e^{-xt}dt,\quad x>0.
$$
Consider the integral
$$
I(b)=\int_0^\infty \frac{x^xe^{-bx}}{\Gamma(x+1)}dx=\int_0^\infty e^{-bx}dx\frac{i}{2\pi}\int_C(-t)^{-x-1}e^{-xt}dt.
$$
Changing the order of integration and calculating the integral over $x$ we get
$$
I(b)=-\frac{i}{2\pi}\int_C\frac{dt}{t(t+\ln(-t)+b)}.
$$
This integral can be easily calculated using residue theory
$$
I(b)=-\frac{1}{1+W_{-1}(-e^{-b})},
$$
where $W_{-1}(z)$ is the Lambert W-function, satisfying $W_{-1}(z)e^{W_{-1}(z)}=z$, which has the derivative
$$
W_{-1}'(z)=\frac{e^{-W_{-1}(z)}}{1+W_{-1}(z)}.\tag{*}
$$
Now we write the initial integral as follows
\begin{align}
\int_0^\infty \frac{x^x e^{-x}}{\Gamma(x+2)}dx&=e\int_1^\infty e^{-b}I(b)db\\
&=-e\int_1^\infty e^{-b}\frac{1}{1+W_{-1}(-e^{-b})}db\\
&=-e\int_0^{1/e} \frac{ds}{1+W_{-1}(-s)}\\
&=e\int_0^{-1/e} \frac{ds}{1+W_{-1}(s)}\\
&=e\int_{-\infty}^{-1} e^{W}d{W}=e\cdot 1/e=1.
\end{align}
Here we used the equation $(*)$ to write $\frac{1}{1+W_{-1}(s)}=e^{W_{-1}(s)}W_{-1}'(s)$ and also the facts $W_{-1}(-1/e)=-1~$, $W_{-1}(0)=-\infty$.
