Summing the first <span class="math">\(n\)</span> terms (arbitrary power)

We know that:

\begin{equation*} \sum_{k=1}^{n} k=\frac{n(n+1)}{2} \end{equation*}

We want to find the expression for

\begin{equation*} \sum_{k=1}^{n} k^{m} \end{equation*}

Recursive Derivation

A general prescription for doing this is looking at one higher power.

\begin{equation*} k^{m+1}-(k-1)^{m+1} \end{equation*}

We can expand this to get an expression that has \(k^{m}\) as its highest power.

Now let’s start putting in values beginning with \(k=n\) going downwards. I’ll use \(m=2\) as an example:

\begin{align*} n^{3}-(n-1)^{3} & = 3n^{2}-3n+1 \\ (n-1)^{3}-(n-2)^{3} & = 3(n-1)^{2}-3(n-1)+1 \\ ... \\ 2^{3}-1^{3} & = 3(2)^{2}-3(1)+1 \\ 1^{3}-0^{3} & = 3(1)^{3}-3(0)+1 \end{align*}

Now add up all these equations. Note that on the LHS, almost all the terms cancel. We end up with:

\begin{equation*} n^{3}=3\sum k^{2}-3\sum k+\sum 1 \end{equation*}

Rearrange for \(\sum k^{2}\) and you have your result.

Clearly, this method works only if you know the formula for all the sums of lower powers.

Below is an algorithm to compute the expression all the way up to \(\sum k^{10}\)

n = var('n')
sums = [n, n*(n+1)/2]

def get_sum(power, sums):
    x = var('x')
    A = x^(power+1) - (x-1)^(power+1)
    coeffs = A.coefficients(sparse=False)
    terms = (n^(power+1) - sum(coeff * s for (coeff, s) in zip(coeffs[:-1], sums)))/coeffs[-1]
    return terms

for index in range(2, 11):
    sums.append(get_sum(index, sums))
    print(f"n^{index}: {sums[-1].factor()}")

n^2: 1/6*(2*n + 1)*(n + 1)*n
n^3: 1/4*(n + 1)^2*n^2
n^4: 1/30*(3*n^2 + 3*n - 1)*(2*n + 1)*(n + 1)*n
n^5: 1/12*(2*n^2 + 2*n - 1)*(n + 1)^2*n^2
n^6: 1/42*(3*n^4 + 6*n^3 - 3*n + 1)*(2*n + 1)*(n + 1)*n
n^7: 1/24*(3*n^4 + 6*n^3 - n^2 - 4*n + 2)*(n + 1)^2*n^2
n^8: 1/90*(5*n^6 + 15*n^5 + 5*n^4 - 15*n^3 - n^2 + 9*n - 3)*(2*n + 1)*(n + 1)*n
n^9: 1/20*(2*n^4 + 4*n^3 - n^2 - 3*n + 3)*(n^2 + n - 1)*(n + 1)^2*n^2
n^10: 1/66*(3*n^6 + 9*n^5 + 2*n^4 - 11*n^3 + 3*n^2 + 10*n - 5)*(n^2 + n - 1)*(2*n + 1)*(n + 1)*n

Special case: \(m=3\)

\begin{equation*} \sum_{k=1}^{n} k^{3}=\frac{n^{2}\left(n+1\right)^{2}}{4}=\left(\sum_{k=1}^{n} k\right)^{2} \end{equation*}

I’m sure I once read an intuitive argument on why the sum of cubes is the square of the sum of the first \(n\) numbers, but I can’t seem to recall it.

Derivation Using A System of Equations

The previous derivation requires you to know the formulae for the sum of all lower powers. What if you want to get the expression for \(\sum k^{10}\) without knowing \(\sum k^{9}\)?

One approach is to guess that it is a polynomial of order \(m+1\) (one order higher). Specifically:

\begin{equation*} \sum_{k=1}^{n}k^{m}=a_{0}+a_{1}n+a_{2}n^{2}\cdots+a_{m+1}n^{m+1} \end{equation*}

and the goal is to determine \(a_{i}\).

I don’t have a rigorous argument on why it would be of one order higher, or why it would even be a polynomial. The only thing I can suggest is realizing that:

\begin{equation*} \int k^{m}\ dk=\frac{1}{m+1}k^{m+1} \end{equation*}

and using that as a guide. Do note that once you have an expression, you can use induction to verify it.

How do we determine \(a_{i}\)? Write out the equations for the first \(m+1\) terms. As an example, for \(m=2\):

\begin{align*} 0^{2} & = 0 & = a_{0} \\ 1^{2} & = 1 & = a_{0}+a_{1}+a_{2}+a_{3} \\ 1^{2}+2^{2} & = 5 & = a_{0}+a_{1}(2)+a_{2}(2)^{2}+a_{3}(2)^{3} \\ 1^{2}+2^{2}+3^{2} & = 14 & = a_{0}+a_{1}(3)+a_{2}(3)^{2}+a_{3}(3)^{3} \end{align*}

We have 4 equations in 4 unknowns. Solve this set of equations. to determine the coefficients.

Here’s some Sage code:

def make_A(power):
    num_terms = power + 2
    rows = []
    for index in range(num_terms):
        row = [index**k for k in range(num_terms)]
        rows.append(row)
    return Matrix(rows)

def make_b(power):
    num_terms = power + 2
    b = []
    for terminal_k in range(num_terms):
        b.append(sum(k**(power) for k in range(1, terminal_k+1)))
    return vector(b)

def get_sum_powers(power):
    """Sum terms of k^power"""
    A = make_A(power)
    b = make_b(power)
    x = A \ b

    n = var('n')
    polynomial = 0
    for index, coeff in enumerate(x):
        polynomial += coeff * n^index
    return polynomial.factor()

for power in range(1, 11):
    print(f"k^{power}: {get_sum_powers(power)}")

k^1: 1/2*(n + 1)*n
k^2: 1/6*(2*n + 1)*(n + 1)*n
k^3: 1/4*(n + 1)^2*n^2
k^4: 1/30*(3*n^2 + 3*n - 1)*(2*n + 1)*(n + 1)*n
k^5: 1/12*(2*n^2 + 2*n - 1)*(n + 1)^2*n^2
k^6: 1/42*(3*n^4 + 6*n^3 - 3*n + 1)*(2*n + 1)*(n + 1)*n
k^7: 1/24*(3*n^4 + 6*n^3 - n^2 - 4*n + 2)*(n + 1)^2*n^2
k^8: 1/90*(5*n^6 + 15*n^5 + 5*n^4 - 15*n^3 - n^2 + 9*n - 3)*(2*n + 1)*(n + 1)*n
k^9: 1/20*(2*n^4 + 4*n^3 - n^2 - 3*n + 3)*(n^2 + n - 1)*(n + 1)^2*n^2
k^10: 1/66*(3*n^6 + 9*n^5 + 2*n^4 - 11*n^3 + 3*n^2 + 10*n - 5)*(n^2 + n - 1)*(2*n + 1)*(n + 1)*n

This should match my previous code.

Derivation Equating Coefficients

Let

\begin{equation*} f(n)=\sum_{k=1}^{n}k^{m} \end{equation*}

Then note that:

\begin{equation*} f(n)-f(n-1)=n^{m} \end{equation*}

Once again, without the best justifications, let’s assume \(f\) is a polynomial of order \(m+1\).

Let’s look at \(m=2\) again. Then \(f(n)=a_{0}+a_{1}n+a_{2}n^{2}+a_{3}n^{3}\). If we write out \(f(n)-f(n-1)\), we get:

\begin{equation*} f(n)-f(n-1)=3a_{3}n^{2} + (2a_{2}-3a_{3})n + a_{1} - a_{2} + a_{3}=n^{2} \end{equation*}

From here we can deduce that:

\begin{equation*} 3a_{3}=1 \end{equation*}

\begin{equation*} 2a_{2}-3a_{3}=0 \end{equation*}

\begin{equation*} a_{1} - a_{2} + a_{3}=0 \end{equation*}

And solve for these coefficients.

Below is Sage code that does just that:

def make_poly_func(power):
    num_terms = power + 2
    n = var('n')
    # For convenience, I'm ignoring a0. We know it will always be 0
    coeff_var_string = " ".join([f"a{index}" for index in range(1, num_terms)])
    coeff_vars = var(coeff_var_string)
    def f(k):
        powers = [k^m for m in range(1, num_terms)]
        return sum(a*km for a, km in zip(coeff_vars, powers))
    return f, coeff_vars

def powers_sum_diff(power):
    f, coeff_vars = make_poly_func(power)
    n = var('n')
    diff = f(n) - f(n-1)
    coeffs = (diff - n^power).coefficients(n, sparse=False)
    # Note: Solution will be something like:
    # [[a1 == (1/6), a2 == (1/2), a3 == (1/3)]]
    solution = solve(coeffs, *coeff_vars)[0]
    # Get just the values.
    solutions = [e.right_hand_side() for e in solution]
    polynomial = 0
    for index, coeff in enumerate(solutions, 1):
        polynomial += coeff * n^index
    return polynomial.factor()

for power in range(1, 11):
    print(f"n^{power}: {powers_sum_diff(power)}")

n^1: 1/2*(n + 1)*n
n^2: 1/6*(2*n + 1)*(n + 1)*n
n^3: 1/4*(n + 1)^2*n^2
n^4: 1/30*(3*n^2 + 3*n - 1)*(2*n + 1)*(n + 1)*n
n^5: 1/12*(2*n^2 + 2*n - 1)*(n + 1)^2*n^2
n^6: 1/42*(3*n^4 + 6*n^3 - 3*n + 1)*(2*n + 1)*(n + 1)*n
n^7: 1/24*(3*n^4 + 6*n^3 - n^2 - 4*n + 2)*(n + 1)^2*n^2
n^8: 1/90*(5*n^6 + 15*n^5 + 5*n^4 - 15*n^3 - n^2 + 9*n - 3)*(2*n + 1)*(n + 1)*n
n^9: 1/20*(2*n^4 + 4*n^3 - n^2 - 3*n + 3)*(n^2 + n - 1)*(n + 1)^2*n^2
n^10: 1/66*(3*n^6 + 9*n^5 + 2*n^4 - 11*n^3 + 3*n^2 + 10*n - 5)*(n^2 + n - 1)*(2*n + 1)*(n + 1)*n

Remember that we arbitrarily guessed that the polynomial is of order one higher than the power? We can drop the assumption in the code, and increase this line to a higher power:

f, coeff_vars = make_poly_func(power+10)

You will see that the results don’t change. However, if you use a smaller power, the solve function returns an empty list.