How to obtain a closed form for summation over polynomial ($\sum_{x=1}^n x^m$)? What is the method for obtaining the polynomial equal to 
\begin{equation*} \sum^{n}_{x=1}x^m \end{equation*}
for unknown $n$, and systematically for various values of $m$? I know it should be a polynomial in $n$ of degree $m+1$, but what are the coefficients?
Wolfram Alpha seems to implement this method, at least for $m$ up to around 40, but it is not clear to me how to reproduce it. Thanks.
 A: Hint: Write an equation with $x^{m+1} - (x-1)^{m+1}$ on the LHS and its expansion using the binomial theorem on the RHS. Now apply the summation operator to both sides, and the LHS telescopes. Rearrange to get $\sum x^m$ as a sum of $\sum x^k$ for each $k < m$.
This proves that $\sum x^m$ is a polynomial. You can then use polynomial interpolation to construct the polynomial.
A: Here  is an  alternate  formulation that  is  slightly different  from
Faulhaber's formula.
Suppose we are trying to evaluate
$$S_m(n) = \sum_{k=1}^n k^m$$
where $m$ is a positive integer.
Because $m$ is positive we may extend this to include zero.

As I have suggested at this MSE link, 
put $$k^m =
\frac{m!}{2\pi i}
\int_{|w|=\epsilon} \frac{1}{w^{m+1}} \exp(kw) \; dw
\\ = \frac{m!}{2\pi i}
\int_{|w|=\epsilon} \frac{1}{w^{m+1}} 
\sum_{q=0}^k {k\choose q} (\exp(w)-1)^q \; dw
\\ = \frac{m!}{2\pi i}
\int_{|w|=\epsilon} \frac{1}{w^{m+1}} 
\sum_{q=0}^k \frac{k!}{(k-q)!} \frac{(\exp(w)-1)^q}{q!} \; dw
\\ = \sum_{q=0}^k \frac{k!}{(k-q)!} \frac{m!}{2\pi i}
\int_{|w|=\epsilon} 
\frac{1}{w^{m+1}} 
\frac{(\exp(w)-1)^q}{q!} \; dw
\\ = \sum_{q=0}^k \frac{k!}{(k-q)!} {m\brace q}.$$
This yields for the sum
$$\sum_{k=0}^n \sum_{q=0}^k \frac{k!}{(k-q)!} {m\brace q}
= \sum_{q=0}^n {m\brace q} \sum_{k=q}^n \frac{k!}{(k-q)!}
\\ = \sum_{q=0}^n q! {m\brace q} 
\sum_{k=q}^n {k\choose q}.$$
Now observe that when $n>m$ the Stirling number is zero for $q>m$, so we may
may set the upper limit on the outer sum to $m$ without loosing any terms. On the other hand when $n<m$ the inner sum is zero when $q>n$ so we may again set the upper limit to $m.$
This yields
$$\sum_{q=0}^m q! {m\brace q} 
\sum_{k=q}^n {k\choose q}.$$
To evaluate the inner sum put
$${k\choose q} =
\frac{1}{2\pi i} 
\int_{|z|=\epsilon}
\frac{(1+z)^k}{z^{q+1}} \; dz.$$
This yields for the sum
$$\frac{1}{2\pi i} 
\int_{|z|=\epsilon}
\frac{1}{z^{q+1}} 
\sum_{k=q}^n (1+z)^k\; dz
= \frac{1}{2\pi i} 
\int_{|z|=\epsilon}
\frac{1}{z^{q+1}} 
\frac{(1+z)^{n+1}-(1+z)^q}{1+z-1}\; dz
\\ = \frac{1}{2\pi i} 
\int_{|z|=\epsilon}
\frac{1}{z^{q+2}} 
((1+z)^{n+1}-(1+z)^q)\; dz.$$
With $n\ge q$ this evaluates to
$${n+1\choose q+1}$$
by inspection.
Returning to the main thread we have established that
$$S_m(n) = 
\sum_{q=0}^m q! {m\brace q} {n+1\choose q+1}.$$
The Stirling number  starts at $q=m$ and thus we  have a polynomial in
$n$ of degree $m+1.$
Recall the OGF of the signed  Stirling numbers of the first kind which
is
$$\sum_{k=1}^n (-1)^{n-k} {n\brack k} z^k = n! \times {z\choose n}.$$
It follows that
$$[n^p] {n+1\choose q+1}
= [n^p] \frac{n+1}{q+1} {n\choose q}
= \frac{1}{q+1} 
\left([n^{p-1}] {n\choose q} +  [n^{p}] {n\choose q} \right)
\\ = \frac{1}{(q+1)!} 
\left((-1)^{q-p} {q\brack p} + (-1)^{q-p-1} {q\brack p-1} \right)
\\ = \frac{(-1)^{q-p}}{(q+1)!} 
\left({q\brack p} - {q\brack p-1} \right).$$
We thus have
$$[n^p] S_m(n)
= \sum_{q=0}^m q! {m\brace q}
\frac{(-1)^{q-p}}{(q+1)!} 
\left({q\brack p} - {q\brack p-1} \right)
\\ = \sum_{q=0}^m {m\brace q}
\frac{(-1)^{q-p}}{q+1} 
\left({q\brack p} - {q\brack p-1} \right).$$
This finally gives the desired polynomial in $n$ for $S_n(m):$
$$S_m(n)
= \sum_{p=1}^{m+1} n^p
\sum_{q=0}^m {m\brace q}
\frac{(-1)^{q-p}}{q+1} 
\left({q\brack p} - {q\brack p-1} \right).$$
