For centuries, mathematicians have studied solutions to sums of powers. The basic form of the series is
\[1^k + 2^k + 3^k + ...+ n^k = \sum_{p=1}^n p^k\]
To solve for any power \(k\) summing to all integers from 1 to any number n, we need separate equations of different degrees for each power \(k\), though there are demonstrable relations between all the equations.
The simplest of the series is raising each number by the power of zero. Since the value of any number to the zeroeth power is 1, the series simply becomes a sum of 1's up to \(n\) inclusive, and so the equation to solve this series is nothing more than \(n\).
\[1^0 + 2^0 + 3^0 + ... + n^0 = n\]
The next power to be considered is \(k=1\).
\[1^1 + 2^1 + 3^1 + ... + n^1\]
This series is clearly just a sum of the integers up to \(n\), and is the basis for the oft-repeated story about Gauss as a schoolboy irritating his teacher by solving it in seconds, after the teacher hoped to saddle the class with the tedium of adding all numbers from 1 to 100, thus earning himself a few minutes of peace and quiet. But Gauss allegedley figured out how to solve it almost instantly. When I first read this story, and since I love puzzles, I immediately closed the book I was reading before it divulged the solution. It took me only seconds to visualize the answer as a spatial manipulation of the number line. Let's take a look at this, as it gives an insight to the algebraic solution to the sum of powers of one, and how that relates to the higher power series solutions.
If we want the sum of all integers from 1 to 100 (which is the same as \(1^1 + 2^1 + 3^1 + ... + 100^1\)), all we have to do is fold the number line in half, doubling it back on itself, and add the columns:
\begin{pmatrix} 1 & 2 & 3 & 4 & ...\,\,\, 50\\ 100 & 99 & 98 & 97 & ...\,\,\, 51\\ \end{pmatrix}
As we can see, each column sums to \(n+1\), or 101. We now have \(\displaystyle\frac{n}{2}\) columns (since we have folded our number line from 1 to \(n\) in half), so the answer is clearly
\[(n+1)*\frac{n}{2}\]
which can also be written as
\[\frac{n^2}{2} + \frac{n}{2}\]
which is precisely the solution for sums of powers of one.
In the 1600s, Johann Faulhaber derived equations to solve summing higher powers up to \(k=17\), based on the equation for powers of \(k=1\). Faulhaber designated the formula above for the power of \(k=1\) as the letter \(a\), so
\[a = \frac{n^2}{2} + \frac{n}{2}\]
He then went on to show that the solutions for the series where \(k\) is odd can be derived from manipulations of \(a\). For example, when \(k=3\), Faulhaber's polynomial solves the series as
\[1^3 + 2^3 + 3^3 + ... + n^3 = a^2\]
When \(k=5\), the solution becomes
\[1^5 + 2^5 + 3^5 + ... + n^5 = \frac{4a^3-a^2}{3}\]
Reducing these equations gives us
\[a^2=\biggl(\frac{n^2}{2} + \frac{n}{2}\biggl)\biggl(\frac{n^2}{2} + \frac{n}{2}\biggl) = \frac{1}{4}n^4+\frac{1}{2}n^3+\frac{1}{4}n^2\]
and
\[\frac{4a^3-a^2}{3}\;=\;\frac{4\biggl(\frac{n^2}{2} + \frac{n}{2}\biggl)\biggl(\frac{n^2}{2} + \frac{n}{2}\biggl)\biggl(\frac{n^2}{2} + \frac{n}{2}\biggl)\,-\,\biggl(\frac{n^2}{2} + \frac{n}{2}\biggl)\biggl(\frac{n^2}{2} + \frac{n}{2}\biggl)}{3} \;=\; \frac{1}{6}n^6+\frac{1}{2}n^5+\frac{5}{12}n^4-\frac{1}{12}n^2\]
The reduced equations are exactly what we will produce below with a simple iterative process requiring only two simple formulas. By repetition of the process we can generate the solutions for the series up to any power \(k\).
It turns out it is fairly easy to create the solutions for sums of powers with only two simple formulas applied in an iterative manner. The creation of the solutions for each value \(k\) will not require the Bernoulli numbers in the iterative formulas; however the Bernoulli numbers for each value \(k\) will automatically be produced through the iterative process. We can generate our solutions within a matrix that is initially populated only with zeros. It can contain as many rows as we would like. We will make it square, though we actually won't need quite as many columns as rows.
\[ \begin{matrix} k = 0 \; \\ k = 1 \; \\ k = 2 \; \\ k = 3 \; \\ k = 4 \; \\ k = 5 \; \\ \end{matrix} % \begin{matrix} \;\; \; \end{matrix} % \begin{pmatrix} 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 \\ \end{pmatrix} \]
This may look like a linear algebra matrix, but we will not treat it as such. Instead, we will use it as a construct within which to perform two simple iterative operations. The first operation will deal only with rows, and the second operaton will deal only with columns, and the process is simply repeated until no further calculations are possible within the construct.
The algorithmic end result of performing these operations within the framework will be to quite easily produce the equations for sums of powers for each power \(k\), from 0 to \(k_{max}\), with \(k_{max}\) being the number of rows minus 1 in our framework matrix. In addition, the row formula below will generate the Bernoulli number for each power \(k\).
\[B_k = 1 - (\frac{a_{c_{1}}}{b_{c_{1}}} + \frac{a_{c_{2}}}{b_{c_{2}}} +\,.. +\, \frac{a_{c_{max}}}{b_{c_{max}}})\]
where \(c_x\) denotes column number.
We first apply this formula to each element in the first row. Because all the elements are zero, the row formula will return a value of 1. This value returned by the formula is also the Bernoulli number for the \(k\) value of the row being calculated, so \(B_k=1\) when \(k=0\).
This value returned by the row formula will constitute the coefficient of \(n\). We will replace the first zero found in the row with \(n\) and its coefficient, which is in the first column. Since our coefficient is 1, and the power of n at this position is 1, the entry would simply be:
$$ \begin{pmatrix} n & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 \\ \end{pmatrix} $$
However, both of our formulas will require the coefficient in the form of a fraction. The power \(p\) to which \(n\) is raised will also be manipulated in the column formula, so we will enter our element in this form:
\[\frac{a}{b}n^{p} \]
Which gives us for our first entry as
$$ \begin{pmatrix} \frac{1}{1}n^1 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 \\ \end{pmatrix} $$The sum of the coefficients of each row will always be equal to \(1\). The value of the coefficient of the element in each row that is only raised to the first power will be the Bernoulli number associated with the the specific power \(k\) for that row. If there is no element in a row that contains \(n\) to the first power, then the coefficient that immediately follows the element \(n^p\) where \(p\geq2\), will be the Bernoulli number for that row. In the latter case the Bernoulli numbers will equal zero, and will occur for series in which \(k\) is odd. So in the first row the Bernoulli number is simply 1, as the first and only element in that row is raised only to the first power, or \(n^1\).
This row is complete because the sum of all coefficients in this row equal 1. And since \(k=0\) for the first row, and the equation for \(1^0 + 2^0 + 3^0 + ... + n^0\) is simply \(\frac{1}{1}n^1\), or \(n\), we have completed the equation for the power \(k=0\).
To continue, we now employ our column formula for the first column. The formula is iterative, and we will apply it sequentially down the column, deriving each element in the column from the element in the row \(r\) directly above. We derive the element in \(r_2\) from the element we have just entered in \(r_1\).
Since the values in the first element are \(a=1\), \(b=1\), \(p=1\), and \(k=0\), we would find the first element in the second row as
\[\frac{1*(k+1)}{1*(p+1)}n^{p+1} = \frac{1*1}{1*2}n^{1+1} = \frac{1}{2}n^{2}\]
Which gives us
$$ \begin{pmatrix} \frac{1}{1}n^1 & 0 & 0 & 0 & 0 & 0 \\ \frac{1}{2}n^2 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 \\ \end{pmatrix} $$
We now repeat the operation again basing the next element down the column in \(r_3\) upon the values held in \(r_2\). This will return the values for \(r_3,c_1\) as
\[\frac{2}{6}n^{3} \]
But we will reduce that fraction for our entry into \(r_3\) as we see below. We repeat this operation with our column formula to calculate the remaining entries, and we get:
$$ \begin{pmatrix} \frac{1}{1}n^1 & 0 & 0 & 0 & 0 & 0 \\ \frac{1}{2}n^2 & 0 & 0 & 0 & 0 & 0 \\ \frac{1}{3}n^3 & 0 & 0 & 0 & 0 & 0 \\ \frac{1}{4}n^4 & 0 & 0 & 0 & 0 & 0 \\ \frac{1}{5}n^5 & 0 & 0 & 0 & 0 & 0 \\ \frac{1}{6}n^6 & 0 & 0 & 0 & 0 & 0 \\ \end{pmatrix} $$
We now use the row formula on the second row, where \(k=1\).
\[B_k = 1 - (\frac{a_{c_{1}}}{b_{c_{1}}} + \frac{a_{c_{2}}}{b_{c_{2}}} +\,.. + \, \frac{a_{c_{max}}}{b_{c_{max}}})\, = \, 1 - \frac{1}{2} \,= \,\frac{1}{2}\]
So Bernoulli number \(B_k\) is automatically generated as \(\frac{1}{2}\), and we have our entry as \(r_2,c_2=\frac{1}{2}n\). This also completes the solution for the sum of powers where \(k=1\), which we have described above by the folding of the number line.
\(1^k + 2^k + 3^k + ... + n^k \,=\,\displaystyle\frac{n^2}{2} + \frac{n}{2}\) for \(k=1\)
We now continue the process by using the column formula again, beginning at the first row in column two with a non-zero value, and fill in all the values with our iterative process. Interestingly, the coefficient of each row in the second column always calculates to \(\frac{1}{2}\).
$$ \begin{pmatrix} \frac{1}{1}n^1 & 0 & 0 & 0 & 0 & 0 \\ \frac{1}{2}n^2 & \frac{1}{2}n\,\, & 0 & 0 & 0 & 0 \\ \frac{1}{3}n^3 & \frac{1}{2}n^2 & 0 & 0 & 0 & 0 \\ \frac{1}{4}n^4 & \frac{1}{2}n^3 & 0 & 0 & 0 & 0 \\ \frac{1}{5}n^5 & \frac{1}{2}n^4 & 0 & 0 & 0 & 0 \\ \frac{1}{6}n^6 & \frac{1}{2}n^5 & 0 & 0 & 0 & 0 \\ \end{pmatrix} $$
Now manipulating the third row with our row formula gives
\[B_k = 1 - (\frac{a_{c_{1}}}{b_{c_{1}}} + \frac{a_{c_{2}}}{b_{c_{2}}} +\,... + \, \frac{a_{c_{max}}}{b_{c_{max}}})\, = \, 1 - \frac{1}{3}+\frac{1}{2} \,= \,\frac{1}{6}\]
$$ \begin{pmatrix} \frac{1}{1}n^1 & 0 & 0 & 0 & 0 & 0 \\ \frac{1}{2}n^2 & \frac{1}{2}n\,\, & 0 & 0 & 0 & 0 \\ \frac{1}{3}n^3 & \frac{1}{2}n^2 & \frac{1}{6}n & 0 & 0 & 0 \\ \frac{1}{4}n^4 & \frac{1}{2}n^3 & 0 & 0 & 0 & 0 \\ \frac{1}{5}n^5 & \frac{1}{2}n^4 & 0 & 0 & 0 & 0 \\ \frac{1}{6}n^6 & \frac{1}{2}n^5 & 0 & 0 & 0 & 0 \\ \end{pmatrix} $$
Which correctly gives the Bernoulli number for \(k=2\) as \(\displaystyle\frac{1}{6}\), and the solution for summing at \(k=2\) as \(\displaystyle\frac{1}{3}n^3+\frac{1}{2}n^2+\frac{1}{6}n\).
So we continue the same process to finish populating our construct. Adding some rows to increase the number of solutions and better see the patterns they describe, we have now completed our calculations, and created and displayed the solutions to the sums of powers from \(k=0\) to \(k=8\), as well as the corresponding Bernoulli numbers associated with each solution.
$$ \begin{pmatrix} \frac{1}{1}n^1 & 0 & 0 & 0 & 0 & 0 \\ \frac{1}{2}n^2 & \frac{1}{2}n\,\, & 0 & 0 & 0 & 0 \\ \frac{1}{3}n^3 & \frac{1}{2}n^2 & \frac{1}{6}n & 0 & 0 & 0 \\ \frac{1}{4}n^4 & \frac{1}{2}n^3 & \frac{1}{4}n^2 & 0 & 0 & 0 \\ \frac{1}{5}n^5 & \frac{1}{2}n^4 & \frac{1}{3}n^3 & -\frac{1}{30}n & 0 & 0 \\ \frac{1}{6}n^6 & \frac{1}{2}n^5 & \frac{5}{12}n^4 & -\frac{1}{12}n^2 & 0 & 0 \\ \frac{1}{7}n^7 & \frac{1}{2}n^6 & \frac{1}{2}n^5 & -\frac{1}{6}n^3 & \frac{1}{42}n & 0 \\ \frac{1}{8}n^8 & \frac{1}{2}n^7 & \frac{7}{12}n^6 & -\frac{7}{24}n^4 & \frac{1}{22}n^2 & 0 \\ \frac{1}{9}n^9 & \frac{1}{2}n^8 & \frac{2}{3}n^7 & -\frac{14}{30}n^5 & \frac{2}{9}n^3 & -\frac{1}{30}n \\ \end{pmatrix} $$
This process can be repeated indefinitely among more rows to generate the Bernoulli numbers for higher values of \(k\) and corresponding series solutions.