How to prove the formula for the Reciprocal Multifactorial constant?

Question

https://mathworld.wolfram.com/ReciprocalMultifactorialConstant.html

The sum of reciprocal multifactorials can be given in closed form by the beautiful formula \begin{align} m(k) &= \sum_{n=0}^\infty\frac1{n\underset{k}{\underbrace{!\dots!}}} \\ &= 1+\frac{e^{1/k}}k\sum_{r=1}^k k^{r/k}\gamma\Big(\frac rk,\frac 1k\Big) \\ &= \frac{e^{1/k}}k\left[k+\sum_{r=1}^{k-1}k^{r/k}\gamma\Big(\frac rk,\frac 1k\Big)\right] \end{align} where $\gamma(a,x)$ is a lower incomplete gamma function (E. W. Weisstein, Aug. 6, 2008). This gives the special cases \begin{align} m(2) &= \sqrt e \left[1+\sqrt{\frac\pi2}\operatorname{erf}\Big(\frac{1}{\sqrt2}\Big)\right] \\ m(3) &= \frac13 e^{1/3}\left[3+3^{1/3}\gamma\Big(\frac13,\frac13\Big)+3^{2/3}\gamma\Big(\frac23,\frac13\Big)\right]. \end{align}

I saw the formula on the above Mathworld link.

Since the reciprocal of $n$ (from $n=0$ to $\infty$) with $k$ number of factorials converge. This page gives a general formula to find the value of the constant for any $k$.

I've tried searching online for the "Reciprocal Multifactorial Constant" and how to derive its formula but I cannot find it anywhere.

Does anyone know of a proof for this formula? Or at least how one can approach proving it?

Answer 1

Let's "decipher" the notation of "multifactorial". Assume $k>0$. Any $n>0$ can be written in the form $n=kq+r$ with $q\geqslant 0$ and $1\leqslant r\leqslant k$; then $n\underbrace{!\ldots!}_{k}:=\prod_{j=0}^q(kj+r)$, and $0\underbrace{!\ldots!}_{k}:=1$. But $$\prod_{j=0}^q(kj+r)=k^{q+1}\prod_{j=0}^q\left(j+\frac rk\right)=k^{q+1}\frac{\Gamma(q+1+r/k)}{\Gamma(r/k)}=\frac{k^{q+1}q!}{\mathrm{B}(r/k,q+1)}$$ using the gamma and beta functions. Using the integral representation of the latter, we get

\begin{align*} m(k)&=1+\sum_{r=1}^k\sum_{q=0}^\infty\frac{1}{(kq+r)\underbrace{!\ldots!}_{k}} \\&=1+\sum_{r=1}^k\sum_{q=0}^\infty\frac{1}{k^{q+1}q!}\mathrm{B}\left(\frac rk,q+1\right) \\&=1+\sum_{r=1}^k\sum_{q=0}^\infty\frac{1}{k^{q+1}q!}\int_0^1 t^{r/k-1}(1-t)^q\,dt \\&=1+\frac1k\sum_{r=1}^k\int_0^1 t^{r/k-1}\sum_{q=0}^\infty\frac1{q!}\left(\frac{1-t}{k}\right)^q dt \\&=1+\frac1k\sum_{r=1}^k\int_0^1 t^{r/k-1}e^{(1-t)/k}\,dt\quad\color{gray}{[t=kx]} \\&=1+\frac{e^{1/k}}k\sum_{r=1}^k k^{r/k}\int_0^{1/k}x^{r/k-1}e^{-x}\,dx, \end{align*}

which is the first of the two closed forms. The second one is then easy to obtain: at $r=k$ $$k^{r/k}\int_0^{1/k}x^{r/k-1}e^{-x}\,dx=k\int_0^{1/k}e^{-x}\,dx=k(1-e^{-1/k}).$$