Now, one can find the digits of e online, or one can ask Mathematica to spit them out, but that wouldn't be as fun as computing them. I had recalled so-called "spigot algorithms" for computing π digit by digit, and I thought now would be a good time to learn more about the subject and apply it to the digits of e.

A quick Google search led to Jeremy Gibbons's "An Unbounded Spigot Algorithm for the Digits of π" paper.

I found that the explanation in the paper jumped around, so I'll explain it in a different way. The idea is that you can take, e.g., a series representation of π and map each portion of that series to a function that transforms an approximation of π into a slightly better approximation. For example, given a series representation like the following (written using Horner's rule):

π = 2 + 1/3 * (2 + 2/5 * (2 + 3/7 * ( ... )))

You can then take each term like (2 + 1/3 * ...) and create a function, say, f₁ that represents it:

f₁ x = 2 + 1/3 x

and continue for the remaining terms in the infinite series:

f₂ x = 2 + 2/5 x
f₃ x = 2 + 3/7 x
f₄ x = 2 + 4/9 x
f₅ x = 2 + 5/11 x

and so on. Then, you can compute π by composing the functions:

π = f₁ . f₂ . f₃ . f₄ ...

Now, here's the cool part. With a little analysis (or even insight) we can determine that any tail of the above expression is going to represent a value in the range [3,4]. I.e., if we drop the first four terms, the resulting expression

f₅ . f₆ . f₇ ...

is going to represent a value somewhere between 3 and 4. If we drop the first million terms, the value is still going to be in [3,4]. Now why is this important? Because it lets us drop all of the terms after some point and replace them with a 3 or 4 to yield a lower or higher bound on the full computation, and as those bounds converge we can extract digits – as many as we desire.

For example, let's drop all of the terms after f₁ and replace them with 3 and 4 to yield lower and upper bounds.

f₁ 3 = 2 + 1/3 * 3 = 3
f₁ 4 = 2 + 1/3 * 4 = 3.333333

So, with just a single term, we know that the full, infinite computation of π must yield a result within [3, 3.333333]. In other words, we know that the first digit of π must be 3!

Let's add another term and see what additional precision we can obtain:

(f₁ . f₂) 3 = 3.0666666
(f₁ . f₂) 4 = 3.2

The tenths digit still differs, but we know it must be a 0, 1, or 2. Let's add a third term and see if we can resolve it:

(f₁ . f₂ . f₃) 3 = 3.10476
(f₁ . f₂ . f₃) 4 = 3.16190

Great! Now we have computed π to two places: 3.1.

That's the basic idea. We can compute successive digits by adding additional terms to the chain of functions until the digits are revealed. The rest is mostly optimization.

One optimization is that we can represent each of the functions by a 2x2 matrix such that composition becomes matrix multiplication. This lets us maintain a single matrix that represents the chain so far. We can compose additional terms by simply multiplying the chain's matrix by the new terms' matrices.

Further, we can "shift out" digits as they are revealed to keep the digit of interest in the ones position. That way, we need only test the integer part of the upper and lower bounds for equality in order to determine whether the next digit is revealed. When it becomes revealed, we subtract it out and multiply the remaining fractional part by ten to bring the next digit of interest into the ones position. Conveniently, another matrix multiplication is all we need.

Now, to turn this method to computing e, I used Google to find Sebah and Gourdon's "The constant e and its computation" and the following series representation:

e = 2 + 1/2 * (1 + 1/3 * (1 + 1/4 * (1 + 1/5) ... ))

I then observed that with the exception of the first term, each term takes values in the range [1,2] back into a subrange thereof. And, we can make the first term follow this rule, too, by computing e – 1 instead of e:

e – 1 = 1 + 1/2 * (1 + 1/3 * (1 + 1/4 * (1 + 1/5) ... ))

With this representation in hand, it was easy to use Gibbons's methods to compute e. I wrote my implementation in my favorite programming language, Haskell, which (I'm pleased to note) Gibbons also used for his π spigots.

All in all, reading Gibbons's paper and implementing the e spigot was the most fun part of solving the Google challenge task (and the follow-up task).

Thank you for reading, and a big thanks to Google for providing an entertaining problem to solve.