There is no closed formula for the given sum
$$ \sum_{1\le k\le n} k^2\cdot k!\ . $$
(At least i cannot see any at the first glance.) As in the comment i committed some seconds in advance, the implemented algorithm works for
$$ \sum_{1\le k\le n} P(k)\cdot k!\ . $$
iff the polynomial $P$ that appears is in the linear span of the polynomials:

- $P_1=(k+1)-1$ ,
- $P_2=(k+2)(k+1)-1$ ,
- $P_3=(k+3)(k+2)(k+1)-1$ ,
- $P_4=(k+4)(k+3)(k+2)(k+1)-1$ and so on.

(They all satisfy the necessary condition $P(-1)=-1$.)

But $P(k)=k^2$ is not in this span. Instead, $k^2+1=P_2-3P_1$ is.

So the following works for some simple polynomials in the span:

```
sage: for P in ( k, k^2+1, k^3-1 ):
....: print ( "P = %-7s and sum( P(k).k! for k=1..n ) is %s"
....: % ( P, sum( P*factorial(k), k, 1, n ) ) )
....:
P = k and sum( P(k).k! for k=1..n ) is factorial(n + 1) - 1
P = k^2 + 1 and sum( P(k).k! for k=1..n ) is n*factorial(n + 1)
P = k^3 - 1 and sum( P(k).k! for k=1..n ) is (n^2 - 2)*factorial(n + 1) + 2
```

**Bonus:** (Out of scope. Please ignore if annoying.)

Let us do the job completely in $\LaTeX$ for polynomials of the shape $k^{power}+$constant:

Code:

```
D = 9 # stop at this dimension of vector space to be used
R.<x> = QQ[]
L = [ R(1), ] + [ prod( x+j for j in [1..k] ) - 1 for k in [1..D-1] ]
L = [ p.coefficients( sparse=0 ) for p in L ]
L = [ p+[ QQ(0) for _ in range(D-len(p)) ] for p in L ]
V = (QQ^D) . span_of_basis( L )
for kk in [1..D-1]:
v = vector( QQ, [ j==kk for j in [0..D-1] ] )
print "%s -> %s" % ( v, V.coordinate_vector(v) )
# and now for the problem
import re
var( 'k,n' );
for kk in [1..D-1]:
v = vector( QQ, [ j==kk for j in [0..D-1] ] )
c = V.coordinate_vector(v)[0]
P = k^kk - c
print ( r" - $ P = %s$ and $\sum_1^n P\cdot k! = %s$"
% ( P, re.sub( r'\\,', '', latex( sum( P*factorial(k), k, 1, n ) ) ) ) )
```

Relevant results are now copy+pasted here:

- $ P = k$ and $\sum_1^n P\cdot k! = \left(n + 1\right)! - 1$
- $ P = k^2 + 1$ and $\sum_1^n P\cdot k! = n \left(n + 1\right)!$
- $ P = k^3 - 1$ and $\sum_1^n P\cdot k! = {\left(n^{2} - 2\right)} \left(n + 1\right)! + 2$
- $ P = k^4 - 2$ and $\sum_1^n P\cdot k! = {\left(n^{3} - 3 n + 3\right)} \left(n + 1\right)! - 3$
- $ P = k^5 + 9$ and $\sum_1^n P\cdot k! = {\left(n^{4} - 4 n^{2} + 6 n + 4\right)} \left(n + 1\right)! - 4$
- $ P = k^6 - 9$ and $\sum_1^n P\cdot k! = {\left(n^{5} - 5 n^{3} + 10 n^{2} + 5 n - 30\right)} \left(n + 1\right)! + 30$
- $ P = k^7 - 50$ and $\sum_1^n P\cdot k! = {\left(n^{6} - 6 n^{4} + 15 n^{3} + 4 n^{2} - 66 n + 55\right)} \left(n + 1\right)! - 55$
- $ P = k^8 + 267$ and $\sum_1^n P\cdot k! = {\left(n^{7} - 7 n^{5} + 21 n^{4} - 119 n^{2} + 175 n + 126\right)} \left(n + 1\right)! - 126$

All above results are computed and displayed via `sage`

.

Note that the following works:

Now somebody has to implement the trick also for all sums of the shape $$\sum_{1\le k\le n} P(k)\cdot k!$$ with an arbitrary polynomial $P$. (Well,

`simplify_factorial`

did not help...)