Sums of Powers

Posted April 15, 2026

This page is all about how we can efficiently compute large sums of powers:

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

Demo

Note that large exponents will take longer to calculate, but your upper bound can be quite large - even thousands of digits!

n

(upper bound)

p

(exponent)

How it works

A straightforward approach to computing a sum like this would be something like

function sumOfPowers(n, p) {
	let total = 0;
	for (let k = 1; k <= n; k++) {
		total += Math.pow(k, p);
	}
	return total;
}

A simple algorithm for computing sums of powers

This works fine for small values of $n$ , but this algorithm requires we perform $n$ exponentiations and $n$ additions. In other words, the time to run grows as $n$ increases, which makes computing sums with a high upper bound take a long time. This is not the algorithm the demo uses.

In fact, the algorithm used in the demo looks quite different. It's not that many lines of code, but it's relying on some nifty (and slightly tricky) mathematics. If all you do is read the code, it's hard to believe it works, so rather than simply presenting it right now, I want to introduce it piece by piece and explain it as we go.

Special cases

If we want to compute a sum where each term's exponent is 1, it turns out there's a simple formula we can use (read more about why it works here):

1 + 2 + 3 + \dots + n = \frac{n (n + 1)}{2}

The nice thing about this formula is the number of operations we perform is no longer dependent on our upper bound - we'll always perform one addition, one multiplication, and one division. That means we can easily determine the sum of the numbers from 1 to some extreme limit like 1,000,000,000,000,000 much faster than with the naive algorithm.

There are similar formulas for computing sums when our exponent is 2 or 3:

1^{2} + 2^{2} + 3^{2} + \dots + n^{2} = \frac{n (n + 1) (2 n + 1)}{6}

It turns out there's a pattern that lets us construct polynomials like these for any exponent, but the pattern isn't super obvious.

Faulhaber's formula

Faulhaber's formula gives the general pattern these polynomials follow. The most commonly given version of the formula involves a special sequence called the Bernoulli numbers.

Of course, there are algorithms to compute the Bernoulli numbers, but since they're rational numbers, there are tradeoffs to consider. If we use floats to store the sequence, we need to be mindful of rounding errors; additionally, after the 258th term, the Bernoulli numbers become too large to store in a 64 bit float. This means to compute sums of large powers precisely, we'll want a rational data type, which not every language has handy.

I suppose if you need the coefficients of the polynomial for some reason, this is the formula for you, but if all you want is the sum of powers for some $n$ and $p$ , we can sidestep computing the Bernoulli numbers by taking another approach.

Stirling numbers

To begin explaining this alternate approach, we need to introduce the Stirling numbers of the second kind (as you might expect, there are Stirling numbers of the first kind too, but they aren't super relevant here). These numbers give the number of ways to partition a set of $n$ elements into $k$ non-empty subsets.

For instance, there's 6 ways to partition the letters a, b, c, and d into 3 non-empty subsets:

{{a,b},{c},{d}}
{{a,c},{b},{d}}
{{a},{b,c},{d}}
{{a,d},{b},{c}}
{{a},{b,d},{c}}
{{a},{b},{c,d}}

The Stirling numbers of the second kind (or just "Stirling numbers" from here on) happen to have a nice recursive formula that we can apply to compute them.

{\binom{n}{k}} = k {\binom{n - 1}{k}} + {\binom{n - 1}{k - 1}}

To compute sums of powers for a given exponent $p$ , we'll need to know

{\binom{p}{0}}, {\binom{p}{1}}, {\binom{p}{2}}, \dots, {\binom{p}{p}}

(We'll learn why in a bit.) Here's a function to get all those values.

function stirling2(p) {
	if (p === 0) {
		return [1n];
	}

	let s = [0n];
	for (let i = 1; i < p + 1; i++) {
		s.push(1n);
		for (let j = i - 1; j > 0; j--) {
			s[j] = BigInt(j) * s[j] + s[j-1];
		}
	}

	return s;
}

An algorithm for computing Stirling numbers of the second kind

This algorithm is quadratic in the value of $p$ , which isn't great, but it's simple to implement (and I don't know of any faster way).

Stirling numbers calculator

Be warned that these numbers can grow quite large, even for small values of $p$ .

p

Connection to exponents

To understand what the Stirling numbers have to do with our problem, it helps to first introduce the notation for falling factorials.

x^{n} = x (x - 1) (x - 2) \dots (x - n + 1)

It turns out there is a very elegant identity relating exponents, Stirling numbers, and falling factorials:¹

x^{n} = {\binom{n}{0}} x^{0} + {\binom{n}{1}} x^{1} + {\binom{n}{2}} x^{2} + \dots + {\binom{n}{n}} x^{n}

We can apply this formula to every term of our power sum and rearrange to get

\begin{matrix} 1^{p} + 2^{p} + 3^{p} + \dots + n^{p} = & {\binom{p}{0}} (1^{0} + 2^{0} + 3^{0} + \dots + n^{0}) + \\ {\binom{p}{1}} (1^{1} + 2^{1} + 3^{1} + \dots + n^{1}) + \\ {\binom{p}{2}} (1^{2} + 2^{2} + 3^{2} + \dots + n^{2}) + \\ \dots + \\ {\binom{p}{p}} (1^{p} + 2^{p} + 3^{p} + \dots + n^{p}) \end{matrix}

This might seem more complicated at first, but despite the similar notation to normal exponentiation, it's much easier to compute sums of falling factorials. For $p > 0$ ,

1^{p} + 2^{p} + 3^{p} + \dots + n^{p} = \frac{(n + 1)^{p + 1}}{p + 1}

(For $p = 0$ , the sum is just $n$ .)² This lets us simplify.³

1^{p} + 2^{p} + 3^{p} + \dots + n^{p} = {\binom{p}{0}} n + {\binom{p}{1}} \frac{(n + 1)^{2}}{2} + {\binom{p}{2}} \frac{(n + 1)^{3}}{3} + {\binom{p}{3}} \frac{(n + 1)^{4}}{4} + \dots + {\binom{p}{p}} \frac{(n + 1)^{p + 1}}{p + 1}

The first term is 0 unless $p = 0$ , which is a trivial case, so we really only need to consider the remaining terms.

Implementation

This formula might still seem a little complicated, but you may be surprised by how simple it is to implement. We already have a method to compute the Stirling numbers we need. With those values available, we can write a simple loop to compute each term and add to our total.

function faulhaber(n, s) {
	const limit = BigInt(s.length);
	if (limit === 1n) {
		return n;
	}

	let q = n + 1n;
	let total = 0n;
	for (let k = 1n; k < limit; k++) {
		q *= n + 1n - k;
		total += s[k] * q / (k + 1n);
	}

	return total;
}

A function that takes an upper bound and Stirling numbers to compute a sum of powers

Note that the loop is dependent on $p$ and not $n$ , which allows us to compute sums with extremely large upper bounds very quickly.

Final thoughts

This approach spends most of its time computing the Stirling numbers. Additionally, it seems like the performance is largely hampered by how large the Stirling numbers get - computing them with a modulus is considerably quicker. The demo is a little clever, in that it only recomputes the Stirling numbers if the value of the exponent changes.

If you would rather use the more commonly known Faulhaber formula, you'll need to compute the Bernoulli numbers. Brent and Harvey give a simple algorithm that is as efficient as how we're calculating the Stirling numbers. The needed binomial coefficients can also be calculated easily.

Furthermore, Brent and Harvey give another algorithm for the Bernoulli numbers that is even more optimized (although it is considerably more complicated). I believe using that algorithm in conjunction with the more traditional formula would be more efficient than the Stirling number approach.

However, in other ways, I like this formula better. For one thing, whichever algorithm you use, Bernoulli numbers are a bit harder to calculate, and require more lines of code. You'll have to compute the binomial coefficients, as well.

One nice thing about the Stirling approach is the terms are integers throughout. I think that makes this algorithm a little better for modular calculations.

Footnotes

¹ For proofs of this identity, see chapter 1.9 of Enumerative Combinatorics by Stanley and chapter 6.1 of Concrete Mathematics by Graham, Knuth, and Patashnik.
² See chapter 2.6 of Concrete Mathematics for an explanation of this identity.
³ A version of this formula is also mentioned in chapter 6.5 of Concrete Mathematics.