Integral of factorial function $$
\mbox{What can we say about the integral}\quad
\int_{0}^{a} x!\,{\rm d}x\ ?.
$$
Or something like
$\displaystyle\int_{0}^{3} x!\ {\rm d}x\ ?$.
 A: This is a very amazing problem. As I said earlier, I did not find any closed form for the integral and then I only performed numerical integrations.
What is surprizing if that $y(a)=\displaystyle\int_{0}^{a} \Gamma(1+x)~ dx$ looks very much to $x(a)=\Gamma(1+a)$ (plot the two curves as a function of $a$).
Taking into account Lucian's comment, I performed a parametric plot  and I observed that, for values of $a \gt 1$, $y$ is almost linear with $x$.
A regression $\log(y)=a +b \log(x)$, for the range $2 \leq a \leq 12$, leads to $a=-0.0775761$ and $b=0.949784$ with $R^2=0.999467$.
A: The factorial function is only defined on the positive integers, so those don't make sense.  However, there is a generalization of the factorial called the Gamma function which you might want to check out.
A: Using:

Taylor Series of Gamma Function:

$$\int_0^x t! dt=  \lim_{c\to\infty} \sum_{n=0}^\infty\sum_{m=0}^\infty\frac{(-1)^{m+1}c^{n+1}\int_0^x t^n dt\text E_{-n}(-c(m+1))}{m!n!} $$
to get the following after an index shift:
$$\int_0^x t!dt= \lim_{c\to\infty} \sum_{m,n=1}^\infty\frac{(-1)^mc^nx^n\text E_{1-n}(-cm))}{\Gamma(m)n!}=\lim_{c\to\infty}\sum_{m,n=1}^\infty\frac{(-1)^{m+n}Q(n,-cm)}{\Gamma(m)n} \left(\frac xm\right)^n$$
Shown here with generalized exponential integral $E_v(z)$ and gamma regularized $Q(a,z)$. Hopefully there is a simplification.
Alternatively, we use the exponential integral Ei$(z)$
$$\int x! dx=\int\lim_{t\to\infty}\sum_{n=0}^\infty\frac{(-1)^n t^{x+n+1}}{n!(n+x+1)}dx=\boxed{C+\lim_{t\to\infty}\sum_{n=1}^\infty\frac{(-1)^{n+1}\text{Ei}((x+n)t)}{\Gamma(n)}}$$
Shown here. When truncating the sum, $n$’s upper bound must be much greater than $t$ to converge.
A: In fact, you want to compute
$$I(a)=\int_{0}^{a} \Gamma(1+x)~ dx=\int_{0}^{a}x\Gamma(x)~ dx$$
Taking into account that
$$\Gamma(x)=\int_{0}^{\infty}y^{x-1}\ e^{-y}~ dy$$
we have
$$I(a)=\int_{0}^{\infty}\ e^{-y}\left (\int_{0}^{a}x y^{x-1}~ dx\right )~ dy$$
The inner integral is easy to calculate
$$\int_{0}^{a}x y^{x-1}~ dx=\frac{y^{a}(a\ln y-1)+1}{y\ln^2 y}$$
The end result
$$I(a)=\int_{0}^{\infty}\frac{y^{a}(a\ln y-1)+1}{y\ln^2 y}e^{-y}~ dy$$
In principle, we are done.
For a particular value of $a$ we can integrate numerically.
But does there exist any analytical expression for $I(a)$ is a separate question.
I'd bet money it doesn't exist.
But I'd be happy if I was wrong.
