Closed form of $\sum_{k=1}^{n}k!$

Question

I was trying to get a closed form of the sum $$\sum_{k=1}^{n}k!$$ Mathematica gives the answer $$(-1)^n\Gamma(n+2)(!(-n-2))-!(-1)-1$$ Here the ! sign before a number is the subfactorial. Here is my try: \begin{align} \sum_{k=1}^{n}k!&=\sum_{k=1}^{n}\Gamma(k+1)\\ &=\sum_{k=1}^{n}\int_{0}^{\infty}x^{k}e^{-x}dx\\ &=\int_{0}^{\infty}e^{-x}\sum_{k=1}^{n}x^kdx\\ &=\int_{0}^{\infty}\frac{e^{-x}(x-x^{n+1})}{1-x}dx \end{align} But I don't know what to do further. So, how can that closed form given by mathematica be derived?
Also, that closed form expression looks complex for some numbers to me. Is it even right?

Answer 1

I believe that this is possibly how Mathematica got its answer.

---

Extend the definition of the sum to negative upper indices: $$ \sum_{k=0}^nf(k)=-\sum_{k=n+1}^{-1}f(k)\tag1 $$ Extend the definition of factorial to negative numbers: $$ n!=\frac{(-1)^{-n-1}}{(-n-1)!}(-1)!\tag2 $$ Applying these extensions to the definition for $!n$, we get $$ \begin{align} !n &=\sum_{k=0}^n(-1)^k\frac{n!}{k!}\tag{3a}\\ &=-\sum_{k=n+1}^{-1}(-1)^k\frac{n!}{k!}\tag{3b}\\ &=-\sum_{k=n+1}^{-1}(-1)^n\frac{(-k-1)!}{(-n-1)!}\tag{3c} \end{align} $$ Explanation:
$\text{(3a)}$: normal definition of $!n$
$\text{(3b)}$: apply $(1)$
$\text{(3c)}$: apply $(2)$

Substituting $n\mapsto-n$ and $k\mapsto-k$, we get $$ !(-n)=\frac{(-1)^{n-1}}{(n-1)!}\sum_{k=1}^{n-1}(k-1)!\tag4 $$ Finally, substituing $n\mapsto n+2$ and $k\mapsto k+1$, we get $$ !(-n-2)=\frac{(-1)^{n+1}}{(n+1)!}\sum_{k=0}^nk!\tag5 $$ which then leads immediately to $$ \sum_{k=0}^nk!=(-1)^{n+1}!(-n-2)\overbrace{\ (n+1)!\ }^{\Gamma(n+2)}\tag6 $$ Formula $(6)$ agrees with what Mathematica gives, using $!(-1)=0$.

I note that my answer and Mathematica's gives a different exponent of $-1$ than in the question. Also, my answer includes $k=0$, so it is $1$ greater.

---

How Mathematica Actually Computes Subfactorial

Let $m\in\mathbb{Z}^+$ and $n+m\gt-1$.

Integrating by parts, we get $$ \frac1{n!}\int_{-1}^\infty x^ne^{-x}\,\mathrm{d}x=\frac{(-1)^n}{n!}e+\frac1{(n-1)!}\int_{-1}^\infty x^{n-1}e^{-x}\,\mathrm{d}x\tag7 $$ Thus, by induction we get $$ \frac1{n!}\int_{-1}^\infty x^ne^{-x}\,\mathrm{d}x=\sum_{k=0}^n\frac{(-1)^k}{k!}e\tag8 $$ and therefore, $$ \begin{align} !n &=\sum_{k=0}^n(-1)^k\frac{n!}{k!}\tag{9a}\\ &=\frac1e\int_{-1}^\infty x^ne^{-x}\,\mathrm{d}x\tag{9b}\\ &=\frac1e\,\Gamma(n+1,-1)\tag{9c}\\ &=\frac1e\,\Gamma(n+1)+\frac{e^{\pi in}}e\int_0^1x^ne^x\,\mathrm{d}x\tag{9d} \end{align} $$ Explanation:
$\text{(9a)}$: equation $(5)$ from this answer
$\text{(9b)}$: apply $(8)$
$\text{(9c)}$: definition of the Incomplete Gamma Function
$\text{(9d)}$: write $(-1)^n=e^{\pi in}$

I think $(-1)^n=e^{-\pi in}$ is just as reasonable a replacement, and better is to average the two, getting $(-1)^n=\cos(\pi n)$. This also has the benefit of keeping $!n\in\mathbb{R}$ for $n\in\mathbb{R}$. However, we will continue with this, as this seems to be what Mathematica does.

First, a recurrence for the integral in $\text{(9d)}$: $$ \begin{align} \frac1{n!}\int_0^1x^ne^x\,\mathrm{d}x &=e\frac1{(n+1)!}-\frac1{(n+1)!}\int_0^1x^{n+1}e^x\,\mathrm{d}x\tag{10a}\\ &=-e\sum_{k=1}^m\frac{(-1)^k}{(n+k)!}+\frac{(-1)^m}{(n+m)!}\int_0^1x^{n+m}e^x\,\mathrm{d}x\tag{10b} \end{align} $$ Explanation:
$\text{(10a)}$: integration by parts
$\text{(10b)}$: induction on $m$ using $\text{(10a)}$ $$ \begin{align} \frac1e\int_0^1x^ne^x\,\mathrm{d}x &=-\sum_{k=1}^m\frac1{k!\binom{-n-1}{k}}+\frac1{m!\binom{-n-1}{m}}\frac1e\int_0^1x^{n+m}e^x\,\mathrm{d}x\tag{11a}\\ &=-\sum_{k=1}^{m-1}\frac1{k!\binom{-n-1}{k}}-\frac1{m!\binom{-n-1}{m}}\left(1-\frac1e\int_0^1x^{n+m}e^x\,\mathrm{d}x\right)\tag{11b} \end{align} $$ Explanation:
$\text{(11a)}$: multiply $(10)$ by $n!/e$ and use $\frac{(-1)^k(n+k)!}{n!}=k!\binom{-n-1}{k}$
$\text{(11b)}$: move the $k=m$ term from the sum to the right
$$ \begin{align} !n &=\overbrace{\frac{(-1)^m\Gamma(n+m+1)}{e\,m!\binom{-n-1}{m}}}^{\frac1e\Gamma(n+1)}\\ &-\sum_{k=1}^{m-1}\frac{e^{\pi in}}{k!\binom{-n-1}{k}}-\frac{e^{\pi in}}{m!\binom{-n-1}{m}}\left(1-\frac1e\int_0^1x^{n+m}e^x\,\mathrm{d}x\right)\tag{12a}\\[6pt] &=\frac{e^{\pi in}\left(e-\int_0^1x^{n+m}e^x\,\mathrm{d}x\right)-(-1)^m\Gamma(n+m+1)}{e\,(m-1)!\binom{-n-1}{m-1}(n+m)}\\ &-\sum_{k=1}^{m-1}\frac{e^{\pi in}}{k!\binom{-n-1}{k}}\tag{12b} \end{align} $$ Explanation:
$\text{(12a)}$: $\text{(9d)}$, $\text{(11b)}$, and $\Gamma(z+1)=z\,\Gamma(z)$
$\text{(12b)}$: commute terms and use $m!\binom{-n-1}{m}=-(m-1)!\binom{-n-1}{m-1}(n+m)$

Let $n\to-m^+$, then $\text{(12b)}$ and L'Hôpital yield $$ \begin{align} !(-m) &=(-1)^m\left[\frac{\gamma-\int_0^1\log(x)\,e^x\,\mathrm{d}x}{e\,(m-1)!}-\sum_{k=1}^{m-1}\frac1{k!\binom{m-1}{k}}\right]+\frac{(-1)^m\pi i}{e\,(m-1)!}\tag{13a}\\ &=\frac{(-1)^m}{(m-1)!}\left[\frac{\gamma+\pi i-\int_0^1\log(x)\,e^x\,\mathrm{d}x}e-\sum_{k=1}^{m-1}(k-1)!\right]\tag{13b} \end{align} $$ where $\gamma$ is the Euler-Mascheroni Constant, which is computed in this answer, and $$ \begin{align} \int_0^1\log(x)\,e^x\,\mathrm{d}x &=\gamma-\operatorname{Ei}(1)\tag{14a}\\ &=-\sum_{k=1}^\infty\frac1{k\,k!}\tag{14b}\\ \end{align} $$ where $\operatorname{Ei}(x)$ is the Exponential Integral.

The values given in $\text{(13b)}$ match those returned by Mathematica for Subfactorial.

Thus, substituting $m\mapsto m+2$ and $k\mapsto k+1$ in $\text{(13b)}$, we get $$ \begin{align} \sum_{k=0}^mk! &=(-1)^{m+1}!(-m-2)(m+1)!+\frac{\gamma+\pi i-\int_0^1\log(x)\,e^x\,\mathrm{d}x}e\tag{15a}\\ &=(-1)^{m+1}!(-m-2)(m+1)!\,-\,!(-1)\tag{15b} \end{align} $$