A power sum is a summmation \(\sum_{k=1}^nk^p\) for some nonnegative integer \(p\). Each one is equal to a polynomial of degree \(p+1\) and asymptotic to \(n^{p+1}/(p+1)\). The following are a few well known power sum formulas.
\[\begin{align} &\sum_{k=1}^n1=n&\sum_{k=1}^nk=\frac{n(n+1)}{2}\\ &\sum_{k=1}^nk^2=\frac{n(n+1)(2n+1)}{6}&\sum_{k=1}^nk^3=\frac{n^2(n+1)^2}{4}\\ \end{align}\]The following are a few more which can be proved.
\[\begin{align} &\sum_{k=1}^nk^4={1\over30}n(n+1)(2n+1)(3n^2+3n-1) ={1\over30}(6n^5+15n^4+10n^3-n)\\ &\sum_{k=1}^nk^5={1\over12}n^2(n+1)^2(2n^2+2n-1) ={1\over12}(2n^6+6n^5+5n^4-n^2)\\ &\sum_{k=1}^nk^6={1\over42}n(n+1)(2n+1)(3n^4+6n^3-3n+1) ={1\over42}(6n^7+21n^6+21n^5-7n^3+n)\\ &\sum_{k=1}^nk^7={1\over24}n^2(n+1)^2(3n^4+6n^3-n^2-4n+2) ={1\over24}(3n^8+12n^7+14n^6-7n^4+2n^2)\\ \end{align}\]There are several methods for generating these formulas. The first method is one I have not found anywhere on the internet, but it is something that others certainly have discovered.
One way is to rewrite the summation of \(p\)th powers as a double summation of \(p-1\) powers. To see how we can do this, consider a table like below.
| \(1^{p-1}\) | ||||
| \(2^{p-1}\) | \(2^{p-1}\) | |||
| \(3^{p-1}\) | \(3^{p-1}\) | \(3^{p-1}\) | ||
| \(\vdots\) | \(\ddots\) | |||
| \(n^{p-1}\) | \(\ldots\) | \(n^{p-1}\) |
This represents rewriting the summation as follows, where each row represents the inner summation.
\[\sum_{i=1}^ni^p=\sum_{i=1}^ni\cdot i^{p-1}=\sum_{i=1}^n\sum_{j=1}^ii^{p-1}\]But now let's change the summation order so the inner summation goes over columns instead of rows. Then we will have the following.
\[\begin{align}&\sum_{i=1}^n\sum_{j=i}^nj^{p-1} =\sum_{i=1}^n\left(\sum_{j=1}^nj^{p-1}-\sum_{j=1}^ij^{p-1}+i^{p-1}\right) =n\sum_{j=1}^nj^{p-1}+\sum_{i=1}^ni^{p-1}-\sum_{i=1}^n\sum_{j=1}^ij^{p-1}\\ &=(n+1)\sum_{i=1}^ni^{p-1}-\sum_{i=1}^n\sum_{j=1}^ij^{p-1}\\ \end{align}\]Now, we can use the formula for summmation of \(p-1\) powers in both summations here. We have to do it a second time for the double summation. What will happen is the leading term is \(i^p/p\) and every other term will be powers of \(i\) from \(1,2,\ldots,p-1\) so we can use formulas for lower powers. We can move that leading term to the other side so we end up with \(\frac{p+1}{p}\sum_{i=1}^ni^p\) and what remains is entirely algebra that can be evaluated using formulas of lower power sums. This is probably easier to understand with an example so here is a derivation of the sum for \(p=2\).
\[\begin{align} &\sum_{i=1}^ni^2=(n+1)\sum_{i=1}^ni-\sum_{i=1}^n\sum_{j=1}^ij =(n+1)\frac{n(n+1)}{2}-\sum_{i=1}^n\frac{i(i+1)}{2}\\ &=\frac{n(n+1)^2}{2}-\frac{1}{2}\sum_{i=1}^ni^2-\frac{1}{2}\sum_{i=1}^ni =\frac{n(n+1)^2}{2}-\frac{1}{2}\frac{n(n+1)}{2}-\frac{1}{2}\sum_{i=1}^ni^2\\ \end{align}\]Now move that summation of \(i^2\) to the other side and we have this.
\[\begin{align} &\frac{3}{2}\sum_{i=1}^ni^2=\frac{n(n+1)^2}{2}-\frac{n(n+1)}{4} =n(n+1)\left(\frac{n+1}{2}-\frac{1}{4}\right)\\ &=n(n+1)\left(\frac{2n+2}{4}-\frac{1}{4}\right)=\frac{n(n+1)(2n+1)}{4}\\ \end{align}\]Then we just multiply each side by \(2/3\) and the result is the well known formula for sum of squares.
\[\sum_{i=1}^ni^2=\frac{n(n+1)(2n+1)}{6}\]Here is another example for \(p=3\).
\[\begin{align}& \sum_{i=1}^ni^3=(n+1)\sum_{i=1}^ni^2-\sum_{i=1}^n\sum_{j=1}^ij^2 =(n+1)\frac{n(n+1)(2n+1)}{6}-\sum_{i=1}^n\frac{i(i+1)(2i+1)}{6}\\& =\frac{n(n+1)^2(2n+1)}{6}-\frac{1}{6}\sum_{i=1}^n(2i^3+3i^2+i) \Rightarrow\\& \frac{4}{3}\sum_{i=1}^ni^3=\frac{n(n+1)^2(2n+1)}{6} -\frac{1}{2}\sum_{i=1}^ni^2-\frac{1}{6}\sum_{i=1}^ni\\& =\frac{n(n+1)^2(2n+1)}{6}-\frac{n(n+1)(2n+1)}{12}-\frac{n(n+1)}{12}\\& =\frac{n(n+1)}{12}\left(2(n+1)(2n+1)-(2n+1)-1\right)\\& =\frac{n(n+1)}{12}\left(4n^2+6n+2-2n-2\right)=\frac{n(n+1)}{12}(4n(n+1)) =\frac{n^2(n+1)^2}{3} \Rightarrow\\& \sum_{i=1}^ni^3=\frac{n^2(n+1)^2}{4} =\left(\frac{n(n+1)}{2}\right)^2=\left(\sum_{i=1}^ni\right)^2 \end{align}\]That last part showing that the sum of cubes is the square of the sum of consecutive integers is also known as Nicomachus's theorem. There is a nice visual proof of it on Wikipedia.
Since the formula for the sum of \(p\)th powers is a \(p+1\) degree polynomial, we could also find the sum \(\sum_{i=1}^ni^p\) for \(n=1,2,\ldots,p+2\) and fit a \(p+1\) degree polynomial to these points. For example, for \(p=2\), we would be looking for a cubic polynomial \(P(n)=an^3+bn^2+cn+d\) that satisfies the following.
\[P(1)=1^2=1,P(2)=1^2+2^2=5,P(3)=1^2+2^2+3^2=14,P(4)=1^2+2^2+3^2+4^2=30\]This would give us the following linear system.
\[\begin{bmatrix}1&1&1&1\\8&4&2&1\\27&9&3&1\\64&16&4&1\\\end{bmatrix} \begin{bmatrix}a\\b\\c\\d\\\end{bmatrix} =\begin{bmatrix}1\\5\\14\\30\end{bmatrix} \]The solution is \((a,b,c,d)=(1/3,1/2,1/6,0)\) which tells us \(\sum_{i=1}^ni^2=\frac{1}{3}n^3+\frac{1}{2}n^2+\frac{1}{6}n\) which is equivalent to the formula already shown. We could also find the linear system another way, which is described on this page. For example, with \(p=3\), we would want the following.
\[\sum_{i=1}^ni^3+(n+1)^3=\sum_{i=1}^{n+1}i^3\]For some polynomial \(P(n)=an^4+bn^3+cn^2+dn\). There is no constant term because we would want \(P(0)=0\). If we substitute this in, we have this.
\[\begin{align}& P(n)+(n+1)^3=P(n+1) \Rightarrow \\& an^4+bn^3+cn^2+dn+(n+1)^3=a(n+1)^4+b(n+1)^3+c(n+1)^2+d(n+1)\\& n^4(a)+n^3(b+1)+n^2(c+3)+n(d+3)\\& \qquad=n^4(a)+n^3(4a+b)+n^2(6a+3b+c)+n(4a+3b+2c+d)+(a+b+c+d)\\& n^3(4a)+n^2(6a+3b)+n(4a+3b+2c)+(a+b+c+d)=n^3+3n^2+3n+1 \end{align}\]By matching the polynomial coefficients, this gives us a linear system in 4 variables.
\[\begin{bmatrix}4&0&0&0\\6&3&0&0\\4&3&2&0\\1&1&1&1\\\end{bmatrix} \begin{bmatrix}a\\b\\c\\d\\\end{bmatrix} =\begin{bmatrix}1\\3\\3\\1\end{bmatrix} \]We have a lower triangle matrix so it can be solved with forward substitution and the solution is \((a,b,c,d)=(1/4,1/2,1/4,0)\). This tells us that \(\sum_{i=1}^ni^3=\frac{1}{4}n^4+\frac{1}{2}n^3+\frac{1}{4}n^2\) which is equivalent to the previously shown formula.
This matrix has the pattern of Pascal's triangle in it. This comes from the way we formed the linear equations. Using this method, we can quickly generate a lower triangle matrix and solve for the polynomial coefficients with forward substitution. This is faster than the simple polynomial interpolation which requires more work to solve a full system with a Vandermonde matrix. In terms of computational complexity, this would be \(O(p^2)\) operations. The first polynomial interpolation and the rewriting as a double summation both require \(O(p^3)\) operations.
To derive this formally, we start with a polynomial and equation. The goal is to find coefficients of the \(p+1\) degree polynomial equivalent to \(\sum_{i=1}^ni^p\).
\[P(n)=\sum_{i=0}^{p+1}c_in^i=c_{p+1}n^{p+1}+c_pn^p+\ldots+c_1n+c_0\] \[(n+1)^p=P(n+1)-P(n)\]Since we will have \(c_0=0\), we could start the summation for the polynomial terms at \(i=1\). Any choice of \(c_0\) will satisfy this equation since they cancel in the subtraction, but \(c_0=0\) is the one that makes sense.
\[\begin{align}& (n+1)^p=\sum_{i=1}^{p+1}c_i(n+1)^i-\sum_{i=1}^{p+1}c_in^i\\& \sum_{i=0}^p{p\choose i}n^i=\sum_{i=1}^{p+1}c_i\sum_{j=0}^i{i\choose j}n^j -\sum_{i=1}^{p+1}c_in^i \end{align}\]On the left are the binomial coefficients for each power of \(n\). To match stuff on the right for forming a linear system, we need to group them by powers of \(n\). Below is a table for the double summation where each \(i\) corresponds to a row and each \(j\) corresponds to a column.
| \(c_1{1\choose0}\) | \(c_1{1\choose1}n\) | ||||
|---|---|---|---|---|---|
| \(c_2{2\choose0}\) | \(c_2{2\choose1}n\) | \(c_2{2\choose2}n^2\) | |||
| \(c_3{3\choose0}\) | \(c_3{3\choose1}n\) | \(c_3{3\choose2}n^2\) | \(c_3{3\choose3}n^3\) | ||
| \(\vdots\) | \(\vdots\) | \(\ddots\) | \(\ddots\) | ||
| \(c_{p+1}{p+1\choose0}\) | \(c_{p+1}{p+1\choose1}n\) | \(\ldots\) | \(\ldots\) | \(c_{p+1}{p+1\choose p}n^p\) | \(c_{p+1}{p+1\choose p+1}n^{p+1}\) |
Notice how the summation for \(P(n)\) cancels the highlighted diagonal of terms shown in the table. So we can rewrite the right side by changing the summation order so the outer summation index corresponds to a column, which groups the powers of \(n\) together.
\[\sum_{i=0}^p{p\choose i}n^i=\sum_{i=0}^pn^i\sum_{j=i+1}^{p+1}c_j{j\choose i}\]Each \(i=0,1,\ldots,p\) corresponds to a linear equation for the coefficient of \(n^i\). We can form a linear system \(Ax=b\) in \(p+1\) variables. Here, \(x\) is the vector of the variables \(c_1,c_2,\ldots,c_{p+1}\), \(b\) is the vector of \({p\choose0},{p\choose1},\ldots,{p\choose p}\), and \(A\) is the \((p+1)\times(p+1)\) matrix with the following entries (\(1\leq i,j\leq p+1\)).
\[A_{ij}=\begin{cases} {j\choose i-1}&j\geq i\\ 0&j<i\\ \end{cases}\]This gives us an upper triangle system shown below which can be solved with backward substitution. In the example, we ordered the variables differently so the system was lower triangular.
\[ \begin{bmatrix} {1\choose0}&{2\choose0}&{3\choose0}&\ldots&{p+1\choose0}\\ 0&{2\choose1}&{3\choose1}&\ldots&{p+1\choose1}\\ 0&0&{3\choose2}&\ldots&{p+1\choose2}\\ \vdots&\vdots&\ddots&\ddots&\vdots\\ 0&0&\ldots&0&{p+1\choose p}\\ \end{bmatrix} \begin{bmatrix}c_1\\c_2\\\vdots\\c_{p+1}\\\end{bmatrix} = \begin{bmatrix}{p\choose0}\\{p\choose1}\\\vdots\\{p\choose p}\\\end{bmatrix} \]This system is guaranteed to have a solution since it is triangular with a main diagonal of nonzero entries. With the solution, the following formula is determined.
\[\sum_{i=1}^ni^p=\sum_{i=1}^{p+1}c_in^i =c_1n+c_2n^2+\ldots+c_pn^p+c_{p+1}n^{p+1}\]Using integrals, we can create a system that can be solved for all the formulas for powers \(0,1,\ldots,p\). It still requires \(O(p^3)\) operations to find all the formulas like the previous methods would with triangular systems. Start with this integral, using the substitutions \(u_i=x+i\) and expanding with the binomial theorem. This method is shown on this page and below is the work shown to derive it in more detail.
\[\begin{align} &\int_0^nx^rdx=\int_0^1x^rdx+\int_1^2x^rdx+\ldots+\int_{n-1}^nx^rdx= \sum_{i=0}^{n-1}\int_i^{i+1}x^rdx\\& =\sum_{i=0}^{n-1}\int_0^1u_i^rdu_i=\int_0^1(x+0)^rdx+\int_0^1(x+1)^rdx +\ldots+\int_{n-1}^n(x+n-1)^rdx\\& =\int_0^1\left[\sum_{i=0}^{n-1}(x+i)^r\right]dx =\int_0^1\left[\sum_{i=0}^{n-1}\sum_{j=0}^r{r\choose j}x^ji^{r-j}\right]dx \end{align}\]Now, swap the summation order (fine since it sums over the rectangle \([0,\ldots,n-1]\times[0,\ldots,r]\)), evaluate the integration, and reverse the order of the terms in the final summation, letting \(S_r=\sum_{i=0}^{n-1}i^r\).
\[\begin{align}& \int_0^1\left[\sum_{j=0}^r\left[{r\choose j}x^j\sum_{i=0}^{n-1}i^{r-j} \right]\right]dx =\int_0^1\left[\sum_{j=0}^r{r\choose j}x^jS_{r-j}\right]dx =\left[\sum_{j=0}^r{r\choose j}{x^{j+1}\over j+1}S_{r-j}\right]_0^1\\& =\sum_{j=0}^r{r\choose j}{S_{r-j}\over j+1} =\sum_{i=0}^r{r\choose i}{S_i\over r+1-i} \end{align}\]Finally, this whole thing is equal to the integral we started with, so we have the following equality.
\[\sum_{i=0}^r{r\choose i}{S_i\over r+1-i}=\int_0^nx^rdx={n^{r+1}\over r+1}\]For each choice of \(r=0,1,\ldots,p\), we can form a linear equation in the variables \(S_0,S_1,\ldots,S_p\). This gives us a linear system \(Ax=b\) where \(x=(S_0,S_1,\ldots,S_p)^T\), \(b=({n^1\over1},{n^2\over2},\ldots,{n^{p+1}\over p+1})^T\) and \(A\) is defined by
\[A_{ij}=\begin{cases} {{i-1\choose j-1}\over i+1-j} & j\leq i \\ 0 & j>i \\ \end{cases}\]The linear system looks like this, and is a triangular system. Although it is triangular, each operation for subtracting a multiple of \(S_i\) takes \(O(p)\) operations since each \(S_i\) will be represented as a polynomial so the overall work is still \(O(p^3)\) operations.
\[\begin{bmatrix} {{0\choose0}\over1}&0&\ldots&\ldots&0\\ {{1\choose0}\over2}&{{1\choose1}\over1}&0&&\vdots\\ {{2\choose0}\over3}&{{2\choose1}\over2}&{{2\choose2}\over1}&\ddots&\vdots\\ \vdots&&&\ddots&0\\ {{p\choose0}\over p+1}&{{p\choose1}\over p}&{{p\choose2}\over p-1}&\ldots &{{p\choose p}\over1}\\ \end{bmatrix} \begin{bmatrix}S_0\\S_1\\\vdots\\S_p\\\end{bmatrix} =\begin{bmatrix}{n^1\over1}\\{n^2\over2}\\\vdots\\ {n^{p+1}\over p+1}\\\end{bmatrix} \]Note that solving this gives the formulas for the summations from \(i=0\) to \(i=n-1\). To adjust for summing to \(i=n\), simply add the appropriate power of \(n\) to the formulas obtained.
This method gives an invertible matrix which, upon inversion, directly gives the coefficients for the polynomial. Since it requires inverting a matrix, it is still \(O(p^3)\) time. Consider the following sum when \(r\) and \(n\) are positive integers.
\[\begin{align}&(1^r-0^r)+(2^r-1^r)+(3^r-2^r)+\ldots+(n^r-(n-1)^r)=n^r-0^r\\& \sum_{i=0}^{n-1}((i+1)^r-i^r)=n^r =\sum_{i=0}^{n-1}\sum_{j=0}^{r-1}{r\choose j}i^j \end{align}\]The difference of \(r\)th powers has been replaced by the appropriate part of the sum of the binomial series, which is every term other than \(i^r\). Next, the summation order can be swapped since the inner summation bounds do not depend on the outer summation index. This allows us to isolate power summations which we will substitute with \(S_r=\sum_{i=0}^{n-1}i^r\).
\[\begin{align}& \sum_{j=0}^{r-1}\sum_{i=0}^{n-1}{r\choose j}i^j =\sum_{j=0}^{r-1}{r\choose j}\sum_{i=0}^{n-1}i^j =\sum_{j=0}^{r-1}{r\choose j}S_j \end{align}\]Now, for each \(r=1,2,\ldots,p+1\), we can create a linear equation in terms of each \(S_r\) as variables. By collecting them together, we have a \((p+1)\times(p+1)\) matrix system \(Ax=b\) with \(x=(S_0,S_1,\ldots,S_p)^T\), \(b=(n^1,n^2,\ldots,n^{p+1})^T\), and \(A\) defined below.
\[A_{ij}=\begin{cases}{i\choose j-1}&j\leq i\\0&j>i\\\end{cases}\]The system looks like what is shown below.
\[\begin{bmatrix} {1\choose0}&0&\ldots&\ldots&0\\ {2\choose0}&{2\choose1}&0&&\vdots\\ {3\choose0}&{3\choose1}&{3\choose2}&\ddots&\vdots\\ \vdots&&&\ddots&0\\ {p+1\choose0}&{p+1\choose1}&{p+1\choose2}&\ldots&{p+1\choose p}\\ \end{bmatrix} \begin{bmatrix}S_0\\S_1\\\vdots\\S_p\\\end{bmatrix} =\begin{bmatrix}n^1\\n^2\\\vdots\\n^{p+1}\\\end{bmatrix} \]Just like the matrix system resulting from the integral, this is for sums up to \(n-1\) so they can be adjusted the same way. That could be avoided by starting with a slightly different summation and working out a similar linear system, this time substituting \(S'_r=\sum_{i=1}^ni^r\).
\[\begin{align}&\sum_{i=1}^n(i^r-(i-1)^r)=n^r =-\sum_{i=1}^n((i-1)^r-i^r) =-\sum_{i=1}^n\sum_{j=0}^{r-1}{r\choose j}i^{j}(-1)^{r-j}\\& =-\sum_{j=0}^{r-1}{r\choose j}(-1)^{r-j}\sum_{i=1}^ni^j =-\sum_{j=0}^{r-1}{r\choose j}(-1)^{r-j}S'_j \end{align}\]This defines a similar matrix system, but with the signs of the numbers in Pascal's triangle alternating. The \(x\) and \(b\) vectors are the same as above and the matrix \(A\) is defined below.
\[A_{ij}=\begin{cases}(-1)^{i+j}{i\choose j-1}&j\leq i\\ 0&j>i\end{cases}\]