In school, you probably had to memorize the formula \[1 + 2 + \dots + n = \frac{n(n+1)}{2}.\] If your math teacher went above and beyond, you may also have seen the formulas \[1^2 + 2^2 + \dots + n^2 = \frac{n(n+1)(2n+1)}{6}\] and \[1^3 + 2^3 + \dots + n^3 = \left(\frac{n(n+1)}{2}\right)^2.\] From these few formulas, we might conjecture that such a formula exists for each natural number $k$. More specifically, for the general sum \[S_k(n) := 1^k + 2^k + \dots + n^k,\] the pattern seems to be that $S_k(n)$ is always a polynomial in $n$ of degree $k + 1$. This follows fairly easily by induction, as we shall now see.

Theorem 1. For all integers $k\geq 0$, there exists a polynomial $p_k$ of degree $k + 1$ such that $p_k(n) = \sum_{j=1}^n j^k$ for all integers $n\geq 0$. Moreover, the leading coefficient of $p_k$ is $\frac{1}{k+1}$. (By convention, the empty sum $\sum_{j=1}^0 j^k$ is defined to be $0$.)

Proof. In the base case $k = 0$, the polynomial $p_0(x) = x$ satisfies the conclusion of the theorem. Indeed, we have $\sum_{j=1}^n j^0 = n$ for all $n$, and the leading coefficient of $p_0$ is $\frac{1}{0 + 1} = 1$.

Now pick an arbitrary value $k\geq 1$, and suppose the theorem holds for all smaller values of $k$. Then for all $n\geq 0$,

\begin{align*} S_{k+1}(n) &= \sum_{j=1}^n j^{k+1} \\ &= \sum_{j=0}^{n-1} (j+1)^{k+1} \\ &= 1 - (n + 1)^{k+1} + \sum_{j=1}^n (j + 1)^{k+1} \\ &= 1 - \sum_{\ell=0}^{k+1} \binom{k+1}{\ell}n^{\ell} + \sum_{j=1}^n \sum_{\ell=0}^{k+1} \binom{k+1}{\ell}j^{\ell}\quad\text{(by the binomial theorem)} \\ &= 1 - \sum_{\ell=0}^{k+1} \binom{k+1}{\ell}n^{\ell} + \sum_{\ell=0}^{k+1}\sum_{j=1}^n \binom{k+1}{\ell}j^{\ell}\quad\text{(interchanging the independent double sums)} \\ &= 1 - \sum_{\ell=0}^{k+1} \binom{k+1}{\ell}n^{\ell} + \sum_{\ell=0}^{k+1} \binom{k+1}{\ell}\sum_{j=1}^n j^{\ell} \\ &= 1 - \sum_{\ell=0}^{k+1} \binom{k+1}{\ell}n^{\ell} + \sum_{\ell=0}^{k+1} \binom{k+1}{\ell}S_{\ell}(n) \\ &= 1 - \sum_{\ell=0}^{k+1} \binom{k+1}{\ell}n^{\ell} + \sum_{\ell=0}^k \binom{k+1}{\ell}S_{\ell}(n) + S_{k+1}(n). \end{align*}

Therefore, \[0 = 1 - \sum_{\ell=0}^{k+1} \binom{k+1}{\ell}n^{\ell} + \sum_{\ell=0}^k \binom{k+1}{\ell}S_{\ell}(n).\] By extracting the $S_k(n)$ term from the last sum and isolating it, we see that \[(k + 1)S_k(n) = \sum_{\ell=0}^{k+1} \binom{k+1}{\ell}n^{\ell} - \sum_{\ell=0}^{k-1} \binom{k+1}{\ell}S_{\ell}(n) - 1.\] Now we invoke the inductive hypothesis. For each $0\leq \ell\leq k - 1$, write the polynomial $p_{\ell}$ generically as \begin{equation}\label{p_k coefficient definition} p_{\ell}(x) = c_0^{(\ell)}x^{\ell + 1} + c_1^{(\ell)}x^{\ell} + \cdots + c_{\ell + 1}^{(\ell)} = \sum_{j=0}^{\ell + 1} c_j^{(\ell)}x^{\ell + 1 - j} \end{equation} where the $c_j^{(\ell)}$ are the coefficients. (Note that the “$(\ell)$” superscript is not an exponent, but rather a second index.) We know that $S_{\ell}(n) = p_{\ell}(n)$ for all $\ell < k$. Hence, the constant term is $c_{\ell + 1}^{(\ell)} = p_{\ell}(0) = S_{\ell}(0) = 0$. Therefore,

\begin{align*} (k + 1)S_k(n) &= \sum_{\ell=0}^{k+1} \binom{k+1}{\ell}n^{\ell} - \sum_{\ell=0}^{k-1} \binom{k+1}{\ell}\sum_{j=0}^{\ell}c_j^{(\ell)}n^{\ell + 1 - j} - 1 \\ &= n^{k+1} + \sum_{\ell=0}^k \binom{k+1}{\ell}n^{\ell} - \sum_{\ell=0}^{k-1} \binom{k+1}{\ell}\sum_{j=0}^{\ell}c_j^{(\ell)}n^{\ell + 1 - j} - 1. \end{align*}

At this point, the proof is completed by letting

\begin{equation}\label{p_k equation} p_k(x) = \frac{1}{k+1}\left[x^{k+1} + \sum_{\ell=0}^k \binom{k+1}{\ell}x^{\ell} - \sum_{\ell=0}^{k-1} \binom{k+1}{\ell}\sum_{j=0}^{\ell}c_j^{(\ell)}x^{\ell + 1 - j} - 1\right]. \end{equation}

Indeed, $p_k$ is a polynomial because it is the sum of polynomials. Also, observe that none of the summations contribute any terms of degree greater than $k$, so the highest-degree term is $x^{k+1}$ with coefficient $\frac{1}{k+1}$. QED.

Our next goal is to study the coefficients $c_j^{(k)}$ of $p_k$ (following the notation introduced in equation \eqref{p_k coefficient definition}). We already know that $c_0^{(k)} = \frac{1}{k+1}$ and that $c_{k+1}^{(k)} = S_k(0) = 0$. Now let $1\leq m\leq k$. We will compute $c_m^{(k)}$ by collecting all the $x^{k + 1 - m}$ terms in the equation \eqref{p_k equation} and summing their coefficients.

The summation $\sum_{\ell=0}^k \binom{k+1}{\ell}x^{\ell}$ contributes one $x^{k + 1 - m}$ term with coefficient $\binom{k+1}{k + 1 - m} = \binom{k + 1}{m}$.

The double sum $\sum_{\ell=0}^{k-1} \binom{k+1}{\ell}\sum_{j=0}^{\ell}c_j^{(\ell)}x^{\ell + 1 - j}$ contributes one $x^{k + 1 - m}$ term with coefficient $\binom{k+1}{\ell}c_{\ell-k+m}^{(\ell)}$ for every iteration of the outer sum $\sum_{\ell=0}^{k-1}$ where $\ell\geq k-m$. To see why, consider which values of $j\in [0, \ell]$ satisfy $\ell + 1 - j = k + 1 - m$. Solving for $j$ gives $j = \ell - k + m$. Since $j\geq 0$, we require $\ell\geq k - m$, and since $j\leq\ell$, we require $m\leq k$ which is already true by assumption. Hence, for each value of $\ell\geq k - m$, we get one $x^{k + 1 - m}$ term when $j = \ell - k + m$.

Therefore, the sum of the coefficients of all $x^{k + 1 - m}$ terms is \begin{equation}\label{coefficient recursive equation} c_m^{(k)} = \frac{1}{k+1}\left[\binom{k + 1}{m} - \sum_{\ell=k-m}^{k-1} \binom{k+1}{\ell}c_{\ell-k+m}^{(\ell)}\right]. \end{equation} At this point, we are ready to compute $c_m^{(k)}$ for specific values of $1\leq m\leq k$. Let’s start with $m = 1$: we have

\begin{align*} c_1^{(k)} &= \frac{1}{k+1}\left[\binom{k+1}{1} - \sum_{\ell=k-1}^{k-1} \binom{k+1}{\ell}c_{\ell-k+1}^{(\ell)}\right] \\ &= \frac{1}{k+1}\left[(k + 1) - \binom{k+1}{k-1}c_0^{(k-1)}\right] \\ &= \frac{1}{k+1}\left((k + 1) - \frac{(k+1)k}{2}\cdot \frac{1}{k}\right) \\ &= 1 - \frac{1}{2} \\ &= \frac{1}{2} \end{align*}

for all $k\geq 1$. So now we know that the second coefficient of $p_k$ is always $\frac{1}{2}$ (for example, $p_1(x) = \frac{1}{2}x^2 + {\color{red}\frac{1}{2}} x$ and $p_2(x) = \frac{1}{3}x^3 + {\color{red}\frac{1}{2}}x^2 + \frac{1}{6}x)$. Let’s compute the $m = 2$ and $m = 3$ cases:

\begin{align*} c_2^{(k)} &= \frac{1}{k+1}\left[\binom{k+1}{2} - \binom{k+1}{k-2}c_0^{(k-2)} - \binom{k+1}{k-1}c_1^{(k-1)}\right] \\ &= \frac{1}{k+1}\left[\frac{(k+1)k}{2} - \frac{(k+1)(k)(k-1)}{6}\cdot\frac{1}{k-1} - \frac{(k+1)k}{2}\cdot\frac{1}{2}\right] \\ &= k\left(\frac{1}{2} - \frac{1}{6} - \frac{1}{4}\right) \\ &= \frac{1}{12}k \end{align*}

and

\begin{align*} c_3^{(k)} &= \frac{1}{k+1}\left[\binom{k+1}{3} - \binom{k+1}{k-3}c_0^{(k-3)} - \binom{k+1}{k-2}c_1^{(k-2)} - \binom{k+1}{k-1}c_2^{(k-1)}\right] \\ &= \frac{1}{k+1}\left[\frac{(k+1)(k)(k-1)}{6} - \frac{(k+1)(k)(k-1)(k-2)}{24}\cdot\frac{1}{k-2} - \frac{(k+1)(k)(k-1)}{6}\cdot\frac{1}{2} - \frac{(k+1)k}{2}\cdot \frac{1}{12}(k-1)\right] \\ &= k(k-1)\left(\frac{1}{6} - \frac{1}{24} - \frac{1}{12} - \frac{1}{24}\right) \\ &= 0. \end{align*}

A pattern is emerging: it looks like $c_j^{(k)}$ is always equal to $k(k-1)\cdots (k-j+2)$ multiplied by some constant. Indeed, we have $\binom{k+1}{j} = (k+1)[{\color{red}k(k-1)\cdots (k-j+2)}]$, and every term in the sum $\sum_{\ell=k-m}^{k-1} \binom{k+1}{\ell}c_{\ell-k+m}^{(\ell)}$ always seems to simplify to $(k+1)[{\color{red}k(k-1)\cdots (k-j+2)}]$ times some constant. Now it is time for another proof by induction.

Theorem 2. There exists a sequence of constants $a_1, a_2, a_3, \dots$ such that for all integers $j, k$ where $1\leq j\leq k$, \[c_j^{(k)} = \frac{k!}{(k-j+1)!}a_j.\]

Remark. Note the order of the quantifiers: the sequence $a_1, a_2, a_3, \dots$ is independent of any variables. They are pure constants. From our computation of $c_1^{(k)}$, $c_2^{(k)}$, and $c_3^{(k)}$, we can see that $a_1 = \frac{1}{2}$, $a_2 = \frac{1}{12}$, and $a_3 = 0$.

Before embarking on the proof, it will be helpful to figure out exactly what the constants $a_j$ are. Fix $k\geq 1$. If the theorem is true, then by taking $j = k$, we have \[c_k^{(k)} = \frac{k!}{(k - k + 1)!}a_k = k!a_k,\] so $a_k = \frac{1}{k!}c_k^{(k)}$. This is what $a_k$ must be if there is any hope of the theorem being true, so now we know how to define the sequence $a_1, a_2, a_3, \dots$ in the proof.

Proof. For all $j\geq 0$, define $a_j = \frac{1}{j!}c_j^{(j)}$. Note that we have defined $a_0 = \frac{1}{0!}c_0^{(0)} = 1$ as well; while $a_0$ is not part of the theorem’s conclusion, we will find this definition useful. We now prove the theorem by induction on $k\geq 1$.

Let $k = 1$ and suppose $1\leq j\leq k$. Then $j = 1$. Hence, $a_j = a_1 = \frac{1}{1!}c_1^{(1)} = c_1^{(1)}$, so \[c_j^{(k)} = c_1^{(1)} = \frac{1!}{(1 - 1 + 1)!}c_1^{(1)} = \frac{k!}{(k - j + 1)!}a_j.\] Therefore, the theorem is true for $k = 1$.

Pick an arbitrary value $k\geq 2$ and suppose the theorem holds for all smaller values of $k$. Fix $1\leq j\leq k$. Using equation \eqref{coefficient recursive equation} and the inductive hypothesis, we have

\begin{align*} c_j^{(k)} &= \frac{1}{k+1}\left[\binom{k + 1}{j} - \sum_{\ell=k-j}^{k-1} \binom{k+1}{\ell}c_{\ell-k+j}^{(\ell)}\right] \\ &= \frac{1}{k+1}\left[\binom{k + 1}{j} - \sum_{\ell=k-j}^{k-1} \binom{k+1}{\ell}\left(\frac{\ell!}{(\ell-(\ell-k+j)+1)!}a_{\ell-k+j}\right)\right] \\ &= \frac{1}{k+1}\left[\binom{k+1}{j} - \sum_{\ell=k-j}^{k-1} \frac{(k+1)!}{(k + 1 - \ell)!}\left(\frac{1}{(k-j+1)!}a_{\ell-k+j}\right)\right] \\ &= \frac{1}{k+1}\left[\binom{k+1}{j} - \frac{(k+1)!}{(k-j+1)!}\sum_{\ell=k-j}^{k-1} \frac{a_{\ell-k+j}}{(k+1-\ell)!}\right] \\ &= \frac{k!}{j!(k+1-j)!} - \frac{k!}{(k-j+1)!}\sum_{\ell=k-j}^{k-1} \frac{a_{\ell-k+j}}{(k+1-\ell)!} \\ &= \frac{k!}{(k-j+1)!}\left(\frac{1}{j!} - \sum_{\ell=k-j}^{k-1} \frac{a_{\ell-k+j}}{(k+1-\ell)!}\right). \end{align*}

The proof will be complete if we can prove the recursive relation \begin{equation}\label{a_j recursive equation} a_j = \frac{1}{j!} - \sum_{\ell=k-j}^{k-1} \frac{a_{\ell-k+j}}{(k+1-\ell)!} \end{equation} for all $1\leq j\leq k$. Note that this is why we had to define a value for $a_0$ at the beginning of the proof; the first term of the summation (where $\ell = k - j$) involves the quantity $a_0$. We split the proof of \eqref{a_j recursive equation} to an independent Lemma.

Lemma 3. Define $a_j = \frac{1}{j!}c_j^{(j)}$ for all $j\geq 0$. Then equation \eqref{a_j recursive equation} holds for all integers $j, k$ where $1\leq j\leq k$.

Proof. This is a simple computation. First, observe that

\begin{align*} \frac{1}{j!} - \sum_{\ell=k-j}^{k-1} \frac{a_{\ell-k+j}}{(k+1-\ell)!} &= \frac{1}{j!} - \sum_{\ell=k-j}^{k-1} \frac{c_{\ell-k+j}^{(\ell-k+j)}}{(\ell-k+j)!(k+1-\ell)!} \\ &= \frac{1}{j!} - \sum_{\ell=0}^{j-1} \frac{c_{\ell}^{(\ell)}}{\ell!(j+1-\ell)!}\quad\text{(by shifting the index $\ell\mapsto \ell-k+j$)} \\ &= \frac{1}{j!} - \sum_{\ell=0}^{j-1} \frac{1}{(j+1)!}\binom{j+1}{\ell}c_{\ell}^{(\ell)} \\ &= \frac{1}{j!}\left(1 - \frac{1}{j+1}\sum_{\ell=0}^{j-1} \binom{j+1}{\ell}c_{\ell}^{(\ell)}\right). \end{align*}

We now invoke equation \eqref{coefficient recursive equation} to deduce that \[c_j^{(j)} = \frac{1}{j+1}\left(\binom{j+1}{j} - \sum_{\ell=0}^{j-1} \binom{j+1}{\ell} c_{\ell}^{(\ell)}\right) = 1 - \frac{1}{j+1}\sum_{\ell=0}^{j-1} \binom{j+1}{\ell}c_{\ell}^{(\ell)}.\] Therefore, \[\frac{1}{j!} - \sum_{\ell=k-j}^{k-1} \frac{a_{\ell-k+j}}{(k+1-\ell)!} = \frac{1}{j!}c_j^{(j)} = a_j\] as required. This completes the proof of the Lemma and consequently Theorem 2. QED.

To conclude, we now have that

\begin{align*} p_k(x) &= \frac{1}{k+1}x^{k+1} + \sum_{j=1}^k c_j^{(k)}x^{k-j+1}\quad\text{(recall that $c_{k+1}^{(k)} = 0$)} \\ &= \frac{1}{k+1}x^{k+1} + \sum_{j=1}^k \frac{k!}{(k-j+1)!}a_jx^{k-j+1} \\ &= \frac{1}{k+1}x^{k+1} + \sum_{j=1}^k \frac{k!}{j!(k-j+1)!}c_j^{(j)}x^{k-j+1} \\ &= \frac{1}{k+1}\left(x^{k+1} + \sum_{j=1}^k \binom{k+1}{j}c_j^{(j)}x^{k-j+1}\right) \\ &= \frac{1}{k+1}\sum_{j=0}^k \binom{k+1}{j}c_j^{(j)}x^{k-j+1}\quad\text{(since $c_0^{(0)} = 1$)}. \end{align*}

for all $k\geq 0$.

The final expression we have obtained for $p_k(x)$ is called Faulhaber’s formula. Note that the linked Wikipedia article uses the notation $B_j^+$ instead of $c_j^{(j)}$; these numbers are commonly known as Bernoulli numbers.

Faulhaber’s formula leads to a fairly easy dynamic-programming algorithm to compute the polynomial $p_k$ (represented as a list of coefficients $[c_0^{(k)}, c_1^{(k)}, \dots, c_k^{(k)}]$) in $O(k^2)$ time and $O(k)$ space.