How is Faulhaber's formula derived? I have been wanting to understand how to find the sum of this series. 

$$1^p + 2^p + 3^p +{\dots} + n^p$$

I am familiar with Gauss' diagonalised adding trick for the sum of the first $n$ natural numbers.

I can prove the formulas for

$$\begin{align} 
 \sum_{1}^{n} k^2 &amp;= \frac{n(n+1)(2n+1)}{6}\\
 \sum_{1}^{n} k^3 &amp;= \frac{n^2(n+1)^2}{4}\\
\sum_{1}^{n} k^4 &amp;= \frac{n(n+1)(2n+1)(3n^2+3n-1)}{30} \\
\end{align}$$

With mathematical induction. But, beyond that even proofs with mathematical induction are difficult.

I'm interested in learning the theory and the proof behind Faulhaber's formula. What is the knowledge required to understand this proof ?
The following interesting derivation is from Aigner's "A Course in Enumeration" (Springer, 2007).

Remember that if we have the following exponential generating functions:

$\begin{align}
   \widehat{A}(z)
     &amp;= \sum_{n \ge 0} a_n \frac{z^n}{n!} \\
   \widehat{B}(z)
     &amp;= \sum_{n \ge 0} b_n \frac{z^n}{n!}
\end{align}$

then also:

$\begin{align}
   \widehat{A}(z) \cdot \widehat{B}(z)
     &amp;= \sum_{n \ge 0}
          \left(
            \sum_{0 \le k \le n} \frac{a_k}{k!} \frac{b_{n - k}}{(n - k)!}
          \right) z^n \\
     &amp;= \sum_{n \ge 0}
          \left(
            \sum_{0 \le k \le n} \frac{n!}{k! (n - k)!} a_k b_{n - k}
          \right) \frac{z^n}{n!} \\
     &amp;= \sum_{n \ge 0}
          \left(
            \sum_{0 \le k \le n} \binom{n}{k} a_k b_{n - k}
          \right) \frac{z^n}{n!}
\end{align}$

Let's define:

$\begin{align}
S_m(n)
      = \sum_{1 \le k \le n - 1} k^m
\end{align}$

We can define the exponential generating function:

$\begin{align}
\widehat{S}_n(z)
      &amp;= \sum_{m \ge 0} S_m(n) \frac{z^m}{m!} \\
      &amp;= \sum_{1 \le k \le n - 1} \sum_{m \ge 0} \frac{k^m z^m}{m!} \\
      &amp;= \sum_{1 \le k \le n - 1} \mathrm{e}^{k z} \\
      &amp;= \frac{\mathrm{e}^{n z} - 1}{\mathrm{e}^z - 1} \\
\end{align}$

This is almost the exponential generating function of the powers of $n$:

$\begin{align}
  \widehat{P}_n(z)
      &amp;= \sum_{m \ge 0} n^m \frac{z^m}{m!} \\
      &amp;= \mathrm{e}^{n z}
\end{align}$

We can write:

$\begin{align}
   (\widehat{P}_n(z) - 1) \widehat{B}(z)
     = z \widehat{S}_n(z) \tag{1}
\end{align}$

where we have the exponential generating function of the Bernoulli numbers:

$\begin{align}
   \widehat{B}(z)
      &amp;= \frac{z}{\mathrm{e}^z - 1} \\
      &amp;= \sum_{n \ge 0} B_n \frac{z^n}{n!}
\end{align}$

Comparing coefficients of $z^{m + 1}$ in (1):

$\begin{align}
 \sum_{m \ge 1} z^m
   \sum_{0 \le k \le m}
      \binom{m}{k} \frac{(n z)^{m - k}}{(m - k)!} B_k
      = \sum_{m \ge 0} S_m(n) \frac{z^{m + 1}}{m!}
\end{align}$

we get after simpĺifying:

$\begin{align}
   S_m(n)
      = \frac{1}{m + 1} \,
           \sum_{0 \le k \le m} \binom{m + 1}{k} B_k n^{m + 1 - k}
\end{align}$

Note: The formula given is often associated with Faulhaber, but Faulhaber's formulas where quite different (and computationally more efficient). This formula is due to Bernoulli.
Here I give an induction formula: from which all sums can be derived:

Let $$ S_n(p)=\sum_{k=1}^{n} k^p\qquad n, p\in\mathbb N ~~~~~\text{called Cavalieri sum of oder p}$$

then, 
We know the following Binomial formula 

$$ (k+1)^p = k^p+ \sum_{i=0}^{p-1}\binom{p}{i} k^i$$
where $\binom{p}{i}= \frac{p!}{i!(p-i)!}$.
Which implies that, 

$$\sum_{k=1}^{n} (k+1)^p =\sum_{k=1}^{n} k^p+\sum_{i=0}^{p-1}\binom{p}{i} \sum_{k=1}^{n} k^i = S_n(p) +\sum_{i=0}^{p-1}\binom{p}{i} S_n(i) $$

But $$\sum_{k=1}^{n} (k+1)^p = \sum_{k=2}^{n+1} k^p = S_{n+1}(p) -1 = S_n(p) +(n+1)^p -1$$

Hence finally we get the formula :

$$\color{red}{(n+1)^p -1  =\sum_{i=0}^{p-1}\binom{p}{i} S_n(i)} $$


  From this it is possible to compute the sum for any $p\ge 1 $ in $ \mathbb N $.