Files
dairydemon.net/junk/faulhaber/faulhaber.html

1080 lines
28 KiB
HTML
Raw Normal View History

2026-04-15 18:15:08 -04:00
<!doctype html>
<html lang="en-US">
<head>
<meta charset="utf-8" />
<meta name="viewport" content="width=device-width,initial-scale=1" />
<link rel="stylesheet" href="/styles/default.css" />
<link rel="stylesheet" href="/styles/junk.css" />
<link rel="stylesheet" href="/styles/faulhaber.css" />
<link rel="icon" href="/images/favicon.ico" />
<title>Sums of Powers</title>
</head>
<body>
<header>
<a id="homelink" href="/">dairydemon.net</a>
</header>
<main>
<article>
<hgroup>
<h1>Sums of Powers</h1>
<p>Posted <time datetime="2026-04-15">April 15, 2026</time></p>
</hgroup>
<p>This page is all about how we can efficiently compute large sums of
powers:</p>
<div class="math-block">
<math display="block">
<mrow>
<msup>
<mn>1</mn>
<mi>p</mi>
</msup>
<mo>+</mo>
<msup>
<mn>2</mn>
<mi>p</mi>
</msup>
<mo>+</mo>
<msup>
<mn>3</mn>
<mi>p</mi>
</msup>
<mo>+</mo>
<mo></mo>
<mo>+</mo>
<msup>
<mi>n</mi>
<mi>p</mi>
</msup>
</mrow>
</math>
</div>
<section>
<h2>Demo</h2>
<p>Note that large exponents will take longer to calculate, but your
upper bound can be quite large - even thousands of digits!</p>
<div id="demo-inputs">
<label id="limit-label" for="limit"><math><mi>n</mi></math> (upper bound)</label>
<input inputmode="numeric" pattern="\d+" id="limit" />
<label id="exponent-label" for="exponent"><math><mi>p</mi></math> (exponent)</label>
<input type="number" id="exponent" />
<input type="button" id="enter" value="Enter" />
</div>
<output id="result" for="limit exponent"></output>
</section>
<section>
<h2>How it works</h2>
<p>A straightforward approach to computing a sum like this would be
something like</p>
<figure>
<pre><code>function sumOfPowers(n, p) {
let total = 0;
for (let k = 1; k &lt;= n; k++) {
total += Math.pow(k, p);
}
return total;
}</code></pre>
<figcaption>A simple algorithm for computing sums of powers</figcaption>
</figure>
<p>This works fine for small values of <math><mi>n</mi></math>, but
this algorithm requires we perform <math><mi>n</mi></math>
exponentiations and <math><mi>n</mi></math> additions. In other
words, the time to run grows as <math><mi>n</mi></math> increases,
which makes computing sums with a high upper bound take a long time.
This is <em>not</em> the algorithm the demo uses.</p>
<p>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.</p>
<section>
<h3>Special cases</h3>
<p>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
(<a href="https://en.wikipedia.org/wiki/Triangular_number">read more about why it works here</a>):</p>
<div class="math-block">
<math display="block">
<mrow>
<mn>1</mn>
<mo>+</mo>
<mn>2</mn>
<mo>+</mo>
<mn>3</mn>
<mo>+</mo>
<mo></mo>
<mo>+</mo>
<mi>n</mi>
<mo>=</mo>
<mfrac>
<mrow>
<mi>n</mi>
<mo>(</mo>
<mi>n</mi>
<mo>+</mo>
<mn>1</mn>
<mo>)</mo>
</mrow>
<mn>2</mn>
</mfrac>
</mrow>
</math>
</div>
<p>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.</p>
<p>There are similar formulas for computing sums when our exponent
is 2 or 3:</p>
<div class="math-block">
<math display="block">
<mrow>
<msup>
<mn>1</mn>
<mn>2</mn>
</msup>
<mo>+</mo>
<msup>
<mn>2</mn>
<mn>2</mn>
</msup>
<mo>+</mo>
<msup>
<mn>3</mn>
<mn>2</mn>
</msup>
<mo>+</mo>
<mo></mo>
<mo>+</mo>
<msup>
<mi>n</mi>
<mn>2</mn>
</msup>
<mo>=</mo>
<mfrac>
<mrow>
<mi>n</mi>
<mo>(</mo>
<mi>n</mi>
<mo>+</mo>
<mn>1</mn>
<mo>)</mo>
<mo>(</mo>
<mn>2</mn>
<mi>n</mi>
<mo>+</mo>
<mn>1</mn>
<mo>)</mo>
</mrow>
<mn>6</mn>
</mfrac>
</mrow>
</math>
<math display="block">
<mrow>
<msup>
<mn>1</mn>
<mn>3</mn>
</msup>
<mo>+</mo>
<msup>
<mn>2</mn>
<mn>3</mn>
</msup>
<mo>+</mo>
<msup>
<mn>3</mn>
<mn>3</mn>
</msup>
<mo>+</mo>
<mo></mo>
<mo>+</mo>
<msup>
<mi>n</mi>
<mn>3</mn>
</msup>
<mo>=</mo>
<msup>
<mrow>
<mo>(</mo>
<mfrac>
<mrow>
<mi>n</mi>
<mo>(</mo>
<mi>n</mi>
<mo>+</mo>
<mn>1</mn>
<mo>)</mo>
</mrow>
<mn>2</mn>
</mfrac>
<mo>)</mo>
</mrow>
<mn>2</mn>
</msup>
</mrow>
</math>
</div>
<p>It turns out there's a pattern that lets us construct
polynomials like these for <em>any</em> exponent, but the pattern
isn't super obvious.</p>
</section>
<section>
<h3>Faulhaber's formula</h3>
<p>
<a href="https://en.wikipedia.org/wiki/Faulhaber%27s_formula">Faulhaber's formula</a>
gives the general pattern these polynomials follow. The most
commonly given version of the formula involves a special sequence
called the
<a href="https://en.wikipedia.org/wiki/Bernoulli_number">Bernoulli numbers</a>.
</p>
<p>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.</p>
<p>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 <math><mi>n</mi></math> and
<math><mi>p</mi></math>, we can sidestep computing the Bernoulli
numbers by taking another approach.</p>
</section>
<section>
<h3>Stirling numbers</h3>
<p>To begin explaining this alternate approach, we need to
introduce the
<a href="https://en.wikipedia.org/wiki/Stirling_numbers_of_the_second_kind">Stirling numbers of the second kind</a>
(as you might expect, there are Stirling numbers of the
<em>first</em> kind too, but they aren't super relevant here).
These numbers give the number of ways to partition a set of
<math><mi>n</mi></math> elements into <math><mi>k</mi></math>
non-empty subsets.</p>
<p>For instance, there's 6 ways to partition the letters a, b, c,
and d into 3 non-empty subsets:</p>
<ul>
<li>{{a,b},{c},{d}}</li>
<li>{{a,c},{b},{d}}</li>
<li>{{a},{b,c},{d}}</li>
<li>{{a,d},{b},{c}}</li>
<li>{{a},{b,d},{c}}</li>
<li>{{a},{b},{c,d}}</li>
</ul>
<p>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.</p>
<div class="math-block">
<math display="block">
<mrow>
<mrow>
<mo>{</mo>
<mfrac linethickness="0px">
<mi>n</mi>
<mi>k</mi>
</mfrac>
<mo>}</mo>
</mrow>
<mo>=</mo>
<mi>k</mi>
<mrow>
<mo>{</mo>
<mfrac linethickness="0px">
<mrow>
<mi>n</mi>
<mo></mo>
<mn>1</mn>
</mrow>
<mi>k</mi>
</mfrac>
<mo>}</mo>
</mrow>
<mo>+</mo>
<mrow>
<mo>{</mo>
<mfrac linethickness="0px">
<mrow>
<mi>n</mi>
<mo></mo>
<mn>1</mn>
</mrow>
<mrow>
<mi>k</mi>
<mo></mo>
<mn>1</mn>
</mrow>
</mfrac>
<mo>}</mo>
</mrow>
</mrow>
</math>
</div>
<p>To compute sums of powers for a given exponent
<math><mi>p</mi></math>, we'll need to know</p>
<div class="math-block">
<math display="block">
<mrow>
<mrow>
<mo fence="true" form="prefix">{</mo>
<mfrac linethickness="0px">
<mi>p</mi>
<mn>0</mn>
</mfrac>
<mo fence="true" form="postfix">}</mo>
</mrow>
<mo separator="true">,</mo>
<mrow>
<mo fence="true" form="prefix">{</mo>
<mfrac linethickness="0px">
<mi>p</mi>
<mn>1</mn>
</mfrac>
<mo fence="true" form="postfix">}</mo>
</mrow>
<mo separator="true">,</mo>
<mrow>
<mo fence="true" form="prefix">{</mo>
<mfrac linethickness="0px">
<mi>p</mi>
<mn>2</mn>
</mfrac>
<mo fence="true" form="postfix">}</mo>
</mrow>
<mo separator="true">,</mo>
<mo></mo>
<mo separator="true">,</mo>
<mrow>
<mo fence="true" form="prefix">{</mo>
<mfrac linethickness="0px">
<mi>p</mi>
<mi>p</mi>
</mfrac>
<mo fence="true" form="postfix">}</mo>
</mrow>
</mrow>
</math>
</div>
<p>(We'll learn why in a bit.) Here's a function to get all those
values.</p>
<figure>
<pre><code>function stirling2(p) {
if (p === 0) {
return [1n];
}
let s = [0n];
for (let i = 1; i &lt; p + 1; i++) {
s.push(1n);
for (let j = i - 1; j &gt; 0; j--) {
s[j] = BigInt(j) * s[j] + s[j-1];
}
}
return s;
}</code></pre>
<figcaption>An algorithm for computing Stirling numbers of the second kind</figcaption>
</figure>
<p>This algorithm is quadratic in the value of
<math><mi>p</mi></math>, which isn't great, but it's simple to
implement (and I don't know of any faster way).</p>
<section>
<h4>Stirling numbers calculator</h4>
<p>Be warned that these numbers can grow quite large, even for
small values of <math><mi>p</mi></math>.</p>
<div id="stirling-input">
<label id="stirling-label" for="stirling-limit"><math><mi>p</mi></math></label>
<input type="number" id="stirling-limit" />
<input type="button" id="stirling-enter" value="Enter" />
</div>
<output id="stirling-result" for="stirling-limit"></output>
</section>
</section>
<section>
<h3>Connection to exponents</h3>
<p>To understand what the Stirling numbers have to do with our
problem, it helps to first introduce the notation for
<a href="https://en.wikipedia.org/wiki/Falling_and_rising_factorials">falling factorials</a>.</p>
<div class="math-block">
<math display="block">
<mrow>
<msup>
<mi>x</mi>
<menclose notation="bottom">
<mi>n</mi>
</menclose>
</msup>
<mo>=</mo>
<mi>x</mi>
<mo>(</mo>
<mi>x</mi>
<mo></mo>
<mn>1</mn>
<mo>)</mo>
<mo>(</mo>
<mi>x</mi>
<mo></mo>
<mn>2</mn>
<mo>)</mo>
<mo></mo>
<mo>(</mo>
<mi>x</mi>
<mo></mo>
<mi>n</mi>
<mo>+</mo>
<mn>1</mn>
<mo>)</mo>
</mrow>
</math>
</div>
<p>It turns out there is a very elegant identity relating
exponents, Stirling numbers, and falling
factorials:<sup>1</sup></p>
<div class="math-block">
<math display="block">
<mrow>
<msup>
<mi>x</mi>
<mi>n</mi>
</msup>
<mo>=</mo>
<mrow>
<mo>{</mo>
<mfrac linethickness="0px">
<mi>n</mi>
<mn>0</mn>
</mfrac>
<mo>}</mo>
</mrow>
<msup>
<mi>x</mi>
<menclose notation="bottom">
<mn>0</mn>
</menclose>
</msup>
<mo>+</mo>
<mrow>
<mo>{</mo>
<mfrac linethickness="0px">
<mi>n</mi>
<mn>1</mn>
</mfrac>
<mo>}</mo>
</mrow>
<msup>
<mi>x</mi>
<menclose notation="bottom">
<mn>1</mn>
</menclose>
</msup>
<mo>+</mo>
<mrow>
<mo>{</mo>
<mfrac linethickness="0px">
<mi>n</mi>
<mn>2</mn>
</mfrac>
<mo>}</mo>
</mrow>
<msup>
<mi>x</mi>
<menclose notation="bottom">
<mn>2</mn>
</menclose>
</msup>
<mo>+</mo>
<mo></mo>
<mo>+</mo>
<mrow>
<mo>{</mo>
<mfrac linethickness="0px">
<mi>n</mi>
<mi>n</mi>
</mfrac>
<mo>}</mo>
</mrow>
<msup>
<mi>x</mi>
<menclose notation="bottom">
<mi>n</mi>
</menclose>
</msup>
</mrow>
</math>
</div>
<p>We can apply this formula to every term of our power sum and
rearrange to get</p>
<div class="math-block">
<math display="block">
<mtable displaystyle="true">
<mtr>
<mtd style="padding-left:0em;padding-right:0em;">
<mrow>
<msup>
<mn>1</mn>
<mi>p</mi>
</msup>
<mo>+</mo>
<msup>
<mn>2</mn>
<mi>p</mi>
</msup>
<mo>+</mo>
<msup>
<mn>3</mn>
<mi>p</mi>
</msup>
<mo>+</mo>
<mo></mo>
<mo>+</mo>
<msup>
<mi>n</mi>
<mi>p</mi>
</msup>
<mo>=</mo>
</mrow>
</mtd>
<mtd style="padding-left:0em;padding-right:0em;">
<mrow>
<mrow>
<mo>{</mo>
<mfrac linethickness="0px">
<mi>p</mi>
<mn>0</mn>
</mfrac>
<mo>}</mo>
</mrow>
<mo>(</mo>
<msup>
<mn>1</mn>
<menclose notation="bottom">
<mn>0</mn>
</menclose>
</msup>
<mo>+</mo>
<msup>
<mn>2</mn>
<menclose notation="bottom">
<mn>0</mn>
</menclose>
</msup>
<mo>+</mo>
<msup>
<mn>3</mn>
<menclose notation="bottom">
<mn>0</mn>
</menclose>
</msup>
<mo>+</mo>
<mo></mo>
<mo>+</mo>
<msup>
<mi>n</mi>
<menclose notation="bottom">
<mn>0</mn>
</menclose>
</msup>
<mo>)</mo>
<mo>+</mo>
</mrow>
</mtd>
</mtr>
<mtr>
<mtd style="padding-left:0em;padding-right:0em;">
<mrow></mrow>
</mtd>
<mtd style="padding-left:0em;padding-right:0em;">
<mrow>
<mrow>
<mo>{</mo>
<mfrac linethickness="0px">
<mi>p</mi>
<mn>1</mn>
</mfrac>
<mo>}</mo>
</mrow>
<mo>(</mo>
<msup>
<mn>1</mn>
<menclose notation="bottom">
<mn>1</mn>
</menclose>
</msup>
<mo>+</mo>
<msup>
<mn>2</mn>
<menclose notation="bottom">
<mn>1</mn>
</menclose>
</msup>
<mo>+</mo>
<msup>
<mn>3</mn>
<menclose notation="bottom">
<mn>1</mn>
</menclose>
</msup>
<mo>+</mo>
<mo></mo>
<mo>+</mo>
<msup>
<mi>n</mi>
<menclose notation="bottom">
<mn>1</mn>
</menclose>
</msup>
<mo>)</mo>
<mo>+</mo>
</mrow>
</mtd>
</mtr>
<mtr>
<mtd style="padding-left:0em;padding-right:0em;">
<mrow></mrow>
</mtd>
<mtd style="padding-left:0em;padding-right:0em;">
<mrow>
<mrow>
<mo>{</mo>
<mfrac linethickness="0px">
<mi>p</mi>
<mn>2</mn>
</mfrac>
<mo>}</mo>
</mrow>
<mo>(</mo>
<msup>
<mn>1</mn>
<menclose notation="bottom">
<mn>2</mn>
</menclose>
</msup>
<mo>+</mo>
<msup>
<mn>2</mn>
<menclose notation="bottom">
<mn>2</mn>
</menclose>
</msup>
<mo>+</mo>
<msup>
<mn>3</mn>
<menclose notation="bottom">
<mn>2</mn>
</menclose>
</msup>
<mo>+</mo>
<mo></mo>
<mo>+</mo>
<msup>
<mi>n</mi>
<menclose notation="bottom">
<mn>2</mn>
</menclose>
</msup>
<mo>)</mo>
<mo>+</mo>
</mrow>
</mtd>
</mtr>
<mtr>
<mtd style="padding-left:0em;padding-right:0em;">
<mrow></mrow>
</mtd>
<mtd style="padding-left:0em;padding-right:0em;">
<mrow>
<mo></mo>
<mo>+</mo>
</mrow>
</mtd>
</mtr>
<mtr>
<mtd style="padding-left:0em;padding-right:0em;">
<mrow></mrow>
</mtd>
<mtd style="padding-left:0em;padding-right:0em;">
<mrow>
<mrow>
<mo>{</mo>
<mfrac linethickness="0px">
<mi>p</mi>
<mi>p</mi>
</mfrac>
<mo>}</mo>
</mrow>
<mo>(</mo>
<msup>
<mn>1</mn>
<menclose notation="bottom">
<mi>p</mi>
</menclose>
</msup>
<mo>+</mo>
<msup>
<mn>2</mn>
<menclose notation="bottom">
<mi>p</mi>
</menclose>
</msup>
<mo>+</mo>
<msup>
<mn>3</mn>
<menclose notation="bottom">
<mi>p</mi>
</menclose>
</msup>
<mo>+</mo>
<mo></mo>
<mo>+</mo>
<msup>
<mi>n</mi>
<menclose notation="bottom">
<mi>p</mi>
</menclose>
</msup>
<mo>)</mo>
</mrow>
</mtd>
</mtr>
</mtable>
</math>
</div>
<p>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
<math><mi>p</mi><mo>&gt;</mo><mn>0</mn></math>,</p>
<div class="math-block">
<math display="block">
<mrow>
<msup>
<mn>1</mn>
<menclose notation="bottom">
<mi>p</mi>
</menclose>
</msup>
<mo>+</mo>
<msup>
<mn>2</mn>
<menclose notation="bottom">
<mi>p</mi>
</menclose>
</msup>
<mo>+</mo>
<msup>
<mn>3</mn>
<menclose notation="bottom">
<mi>p</mi>
</menclose>
</msup>
<mo>+</mo>
<mo></mo>
<mo>+</mo>
<msup>
<mi>n</mi>
<menclose notation="bottom">
<mi>p</mi>
</menclose>
</msup>
<mo>=</mo>
<mfrac>
<mrow>
<mo>(</mo>
<mi>n</mi>
<mo>+</mo>
<mn>1</mn>
<msup>
<mo>)</mo>
<menclose notation="bottom">
<mrow>
<mi>p</mi>
<mo>+</mo>
<mn>1</mn>
</mrow>
</menclose>
</msup>
</mrow>
<mrow>
<mi>p</mi>
<mo>+</mo>
<mn>1</mn>
</mrow>
</mfrac>
</mrow>
</math>
</div>
<p>(For <math><mi>p</mi><mo>=</mo><mn>0</mn></math>, the sum is
just <math><mi>n</mi></math>.)<sup>2</sup> This lets us
simplify.<sup>3</sup></p>
<div class="math-block">
<math display="block">
<mrow>
<msup>
<mn>1</mn>
<mi>p</mi>
</msup>
<mo>+</mo>
<msup>
<mn>2</mn>
<mi>p</mi>
</msup>
<mo>+</mo>
<msup>
<mn>3</mn>
<mi>p</mi>
</msup>
<mo>+</mo>
<mo></mo>
<mo>+</mo>
<msup>
<mi>n</mi>
<mi>p</mi>
</msup>
<mo>=</mo>
<mrow>
<mo>{</mo>
<mfrac linethickness="0px">
<mi>p</mi>
<mn>0</mn>
</mfrac>
<mo>}</mo>
</mrow>
<mi>n</mi>
<mo>+</mo>
<mrow>
<mo>{</mo>
<mfrac linethickness="0px">
<mi>p</mi>
<mn>1</mn>
</mfrac>
<mo>}</mo>
</mrow>
<mfrac>
<mrow>
<mo form="prefix" stretchy="false" lspace="0em" rspace="0em">(</mo>
<mi>n</mi>
<mo>+</mo>
<mn>1</mn>
<msup>
<mo>)</mo>
<menclose notation="bottom">
<mn>2</mn>
</menclose>
</msup>
</mrow>
<mn>2</mn>
</mfrac>
<mo>+</mo>
<mrow>
<mo>{</mo>
<mfrac linethickness="0px">
<mi>p</mi>
<mn>2</mn>
</mfrac>
<mo>}</mo>
</mrow>
<mfrac>
<mrow>
<mo form="prefix" stretchy="false" lspace="0em" rspace="0em">(</mo>
<mi>n</mi>
<mo>+</mo>
<mn>1</mn>
<msup>
<mo>)</mo>
<menclose notation="bottom">
<mn>3</mn>
</menclose>
</msup>
</mrow>
<mn>3</mn>
</mfrac>
<mo>+</mo>
<mrow>
<mo>{</mo>
<mfrac linethickness="0px">
<mi>p</mi>
<mn>3</mn>
</mfrac>
<mo>}</mo>
</mrow>
<mfrac>
<mrow>
<mo form="prefix" stretchy="false" lspace="0em" rspace="0em">(</mo>
<mi>n</mi>
<mo>+</mo>
<mn>1</mn>
<msup>
<mo>)</mo>
<menclose notation="bottom">
<mn>4</mn>
</menclose>
</msup>
</mrow>
<mn>4</mn>
</mfrac>
<mo>+</mo>
<mo></mo>
<mo>+</mo>
<mrow>
<mo>{</mo>
<mfrac linethickness="0px">
<mi>p</mi>
<mi>p</mi>
</mfrac>
<mo>}</mo>
</mrow>
<mfrac>
<mrow>
<mo form="prefix" stretchy="false" lspace="0em" rspace="0em">(</mo>
<mi>n</mi>
<mo>+</mo>
<mn>1</mn>
<msup>
<mo>)</mo>
<menclose notation="bottom">
<mrow>
<mi>p</mi>
<mo>+</mo>
<mn>1</mn>
</mrow>
</menclose>
</msup>
</mrow>
<mrow>
<mi>p</mi>
<mo>+</mo>
<mn>1</mn>
</mrow>
</mfrac>
</mrow>
</math>
</div>
<p>The first term is 0 unless
<math><mi>p</mi><mo>=</mo><mn>0</mn></math>, which is a trivial
case, so we really only need to consider the remaining terms.</p>
</section>
<section>
<h3>Implementation</h3>
<p>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.</p>
<figure>
<pre><code>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 &lt; limit; k++) {
q *= n + 1n - k;
total += s[k] * q / (k + 1n);
}
return total;
}</code></pre>
<figcaption>A function that takes an upper bound and Stirling
numbers to compute a sum of powers</figcaption>
</figure>
</section>
<p>Note that the loop is dependent on <math><mi>p</mi></math> and not
<math><mi>n</mi></math>, which allows us to compute sums with
extremely large upper bounds very quickly.</p>
<section>
<h3>Final thoughts</h3>
<p>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.</p>
<p>If you would rather use the more commonly known Faulhaber
formula, you'll need to compute the Bernoulli numbers.
<a href="https://arxiv.org/abs/1108.0286v3">Brent and Harvey</a>
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.</p>
<p>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.</p>
<p>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.</p>
<p>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.</p>
</section>
</section>
<section>
<h2>Footnotes</h2>
<ol id="footnote-list">
<li><sup>1</sup> For proofs of this identity, see chapter 1.9 of
<cite>Enumerative Combinatorics</cite> by Stanley and chapter 6.1
of <cite>Concrete Mathematics</cite> by Graham, Knuth, and
Patashnik.</li>
<li><sup>2</sup> See chapter 2.6 of <cite>Concrete
Mathematics</cite> for an explanation of this identity.</li>
<li><sup>3</sup> A version of this formula is also mentioned in
chapter 6.5 of <cite>Concrete Mathematics</cite>.</li>
</ol>
</section>
</article>
</main>
<footer>
<p id="copyright">© 2026 filifa. This page is licensed under <a href="https://creativecommons.org/licenses/by-sa/4.0/">CC BY-SA 4.0</a> and the <a href="https://spdx.org/licenses/0BSD.html">BSD Zero Clause license</a>.</p>
</footer>
<script src="faulhaber.js"></script>
</body>
</html>