Spigot Algorithms

Introduction

Recently I was talking with my girlfriend about how her unit calculator tool, Rink, calculates $\pi$. This reminded me of some algorithms that are not very well known but are very interesting, so I decided to search about them more but I couldn't find much information apart from a section in a Wikipedia article and a few research papers. So I'm writing this post to condense all the information I could find or derive myself, explaining everything as thoroughly as it is reasonable in the hope it will be interesting.

There are many widely known methods for approximating $\pi$ and other important irrational numbers with infinite series, continued fractions or other procedures involving limits. The following alternating series was independently discovered by Madhava of Sangamagrama in the 14-th century, James Gregory in 1671 and Gottfried Leibniz in 1673 and is known as the Leibniz formula for $\pi$

How good is this formula for approximating $\pi$? Sadly, not very good, as it needs five thousand terms just to get four decimal places and literally billions of terms for ten decimal places. In constrast, Ramanujan's formula for $\pi$ converges extremely fast:

Indeed, it is so fast that it only needs one term to calculate $\pi$ up to six decimal places and two terms for fifteen. This formula is an example of a Ramanujan-Sato series which are the basis of the Chudnovsky algorithm, an extraordinarily fast algorithm for calculating $\pi$, used for many record breaking computations of up to 202 trillion digits in 2024 using y-cruncher, a very well optimized program that computes $\pi$ and other constants.

Although these series converge really fast, they all suffer from two problems:

• They require floating point numbers. Floating point operations suffer from precision issues and it can be hard to predict how much precision we need to get a specific number of digits;
• You need to compute all the earlier digits before computing some later digit. It is not possible to extract the billionth digit of $\pi$ using a traditional algorithm without computing and storing the 999,999,999 other digits.

The second problem may sound very esoteric, you may have never considered that it is even possible to extract any arbitrary digit of an irrational number without calculating previous digits, but I will soon show how this is, very surprisingly, actually possible. As for the first problem, there is a way to compute $\pi$ digit by digit without any floating point operations, but it won't be faster than an algorithm like Chudnovsky's due to slower convergence.

We call algorithms that approximate a number digit by digit spigot algorithms, and if the algorithm also allows the digits to be calculated without first calculating the earlier digits it is a digit extraction algorithm. Before presenting a spigot algorithm for $\pi$ I will first use Euler's number as an example of how we can construct such algorithms since it is easier than the algorithm for $\pi$.

The Idea

Euler's number is defined as

Euler proved that

and we can rewrite this as

We know that $2$ is the integer part of $e$ so of course the remaining expression is its fractional part, which means we can multiply it by $10$ to shift the next digit to the integer part. For this example it is enough to use only six terms:

$a/b = q + r/b$ for any integer $a,b$ where $q$ is the quotient and $r$ the remainder so we can apply this to all the fractions from right to left like this:

and, indeed, $7$ is the next digit of $e$. You can repeat this procedured as many times as you want and it will keep producing more digits as long as you use enough terms. In general, if we have a number defined by a recursive sum of the form

we can apply this method to produce $a$ digit by digit, $a_0$ is always its integer part and the rest of the recursive sum is its fractional part. So, by induction, if we apply the procedure above we will get the next digit $d_1$ so the integer part of $10a$ is $a_0d_1$ and the rest of the expression is the fractional part, whose first digit is of course the second digit of $a$. So we can keep applying this procedure to generate $a$ digit by digit. Hopefully you can see this pattern by looking at the example. And as it turns out, computers are very good at patterns!

However, this is not necessarily true for the most general case when the recursive sum is not of the above form. In a sum with different fractions it is possible that the sum after $a_0$ is not always the fractional part. For example, a similar formula for $\pi$ is

and of course $2$ is not the integer part of $\pi$, so the rest of the sum still contributes to the integer part. In fact, after we get some number of digits the remaining sum will again contribute to the last digit we produced. This is the complication I mentioned in the Introduction section and we will see how to deal with it later.

In the next section I will formalize the idea presented above and prove that it works as expected for sums of the form (*).

Generalizing Number Bases

We are all used to numbers written in base 10, such as $1729$, $277777788888899$ and $5.28571428571$, and I believe many of you are also fluent in base 2 and base 16. In general, for an integer $b > 1$ then we have a positional numeral system of base $b$ where a string of digits $d_n \dots d_2d_1d_0.d_{-1}d_{-2}d_{-3} \dots$ represents the number

and when $d_i < b$ for all $i$ we say that the representation is regular. This representation does not need to be unique. To see this, let's write $0.999...$ as

which follows from the definition given above. We also have the inequality

Since $\lim_{n \to \infty} 1/10^n = 0$, by the squeeze theorem

Therefore $0.\overline{9} = 1$. And the same is true for any other base, $0.\overline{1} = 1$ in binary and $0.\overline{F} = 1$ in hexadecimal.

Rewriting (4) as

gives us a way to generalized numeral systems. Instead of having a base based on a fixed $b$, we can use any numbers for each $b$ and $b^{-1}$ in the above expression, giving us

where $\dots, b_{-3}, b_{-2}, b_{-1}, b_1, b_2, b_3, \dots, b_n$ are all possibly different numbers, maybe even fractions. We write this base as $(b_n, \dots, b_3, b_2, b_1 ; b_{-1}, b_{-2}, b_{-3}, \dots)$, and we're going to omit the integer part when it is in decimal

So we can use the base $\boldsymbol{b} = (; 2, 3, 4, \dots)$, that is, the radices are all the natural numbers greater than $1$, to write (2) as $(a_0 ; a_1, a_2, a_3\cdots)_{\boldsymbol{b}}$. In particular, if we set $a_0$ to $2$ and $a_{-1}, a_{-2}, a_{-3}, \dots$ to $1$ we can write Euler's number in base $\boldsymbol{b}$ as $(2 ; 1, 1, 1, \dots)_{\boldsymbol{b}}$. We have created a numeral system where an irrational number can be periodic! We could even write Euler's number as $2.\overline{1}$ if the context is clear.

George Cantor was the first to formalize and generalize mixed radix numeral systems as we just defined, publishing his work in Über einfache Zahlensysteme by the Zeitschrift für Math journal. Remember that these numeral systems are nothing but a way to write a number, and as we saw earlier numbers don't even need to have a unique representation, so this is a perfectly valid way of representing a number.

The same way the representation of a number in some fixed natural radix is regular when all the digits are smaller than the base, in mixed-radix a representation is regular when each $d_i$ is smaller than the numerator of $b_i$. So if we now look at (1) we can see that multiplication by an integer in mixed-radix works the same way as in fixed-base: we perform the multiplication in the right-most digit $d_i$ modulo the numerator of $b_i$ and carry the quotient of $d_i / b_i$, add the quotient to the next digit and repeat for the next digits.

For example, in base 10 to do $0.87 * 2$ we multiply $7 * 2 \equiv 4 \bmod 10$, so the first digit is $4$, and carry the quotient of $7 * 2 / 10$, which is $1$, and we keep doing this until we get the integer part $1$ and the new fractional part $74$. So the manipulations we did in (1) are equivalent to multiplying in base $\boldsymbol{b}$! This becomes even clearer when we write $0.87$ in base 10 as

I'm not going to develop mixed-radix systems any further because we already have all we need to continue.

The Algorithm for $e$

As we disucussed, the only problem with our main idea is that the rest of the recursive sum might not actually be the fractional part of the number in question, it might be the case that the sum after $a_0$ in the first expression of (2) is greater than $1$. What this means in terms of base $\boldsymbol{b}$ is that $(0;a_1, a_2, \dots)$ might not be less than $1$ in decimal even though it has integer part $0$ in $\boldsymbol{b}$.

So we want to know when $(0; a_1, a_2, \dots) \geq 1$. This is a special case of the following lemma.

Lemma 1. If $i \geq 1$ then $(0; 0, 0, \dots, 0, a_i, a_{i + 1}, \dots)_{\boldsymbol{b}} \leq \frac{1}{i}$

Proof. This is equivalent to showing that

so we let $a_i, a_{i + 1}, \dots$ to be their maximum values and we get

Since

the series is absolutely convergent by the ratio test and we are allowed to rearrange it:

The fact that $1/i! \leq 1/i$ completes the proof.

Corollary 1.1. $(0;a_1a_2a_3\dots)_{\boldsymbol{b}} = 1$ if and only if $a_i = b_i - 1$ for all $i$, and in this case the integer part of $(a_0;a_1a_2a_3\dots)_{\boldsymbol{b}}$ is $a_0 + 1$, otherwise the integer part is $a_0$. The fractional part is never greater than $1$.

The above corollary implies that if we choose some amount of precision to compute $e$ like we did in the introduction section, which is equivalent to setting $(a_i, a_{i + 1}, \dots)$ to $0$ for some $i \geq 1$, and apply the procedure we will never get an integer greater than $9$. The only question left is how many digits in $\boldsymbol{b}$ we actually need. It turns out that we need $n$ digits is enough, and the proof of this fact is in the proof of the correctness of the following algorithm.

Algorithm A. (Generate the digits of Euler's number in base 10). Given the number of desired digits $\texttt{N}$ this algorithm outputs $\texttt{N}$ digits of $e$.

• Step 1: Set $\texttt{q}$ $\leftarrow$ 2.
• Step 2: Set $\texttt{l}$ $\leftarrow \texttt{N} + 1$.
• Step 3: Set $\texttt{A[t]}$ to 1 for $0 \leq t \leq \texttt{l}$.
• Step 4: Set $\texttt{i}$ to $1$.
• Step 5: Output $\texttt{q}$.
• Step 6: Set $\texttt{A[\(t\)]}$ $\leftarrow 10\texttt{A[t]}$ for $0 \leq t \leq \texttt{l}$
• Step 7: If $\texttt{i} = \texttt{N}$ the algorithm terminates.
• Step 8: Set $\texttt{j} \leftarrow \texttt{l}$.
• Step 9: If $\texttt{j} < 0$ set $\texttt{i} \leftarrow \texttt{i} + 1$ and go to step 5.
• Step 10: Divide $\texttt{A[j]}$ by $\texttt{j} + 2$ and let $\texttt{q}$ be the quotient and $\texttt{r}$ the remainder.
• Step 11: Set $\texttt{A[j]} \leftarrow \texttt{r}$.
• Step 12: Set $\texttt{A[j - 1]} \leftarrow \texttt{A[j - 1]} + \texttt{q}$
• Step 13: Set $\texttt{j} \leftarrow \texttt{j} + 1$.
• Step 14: Go to step 7.

Theorem 1. With input $n$ Algorithm A produces $n$ digits of $e$ in $\Theta(n^2)$ time and $\Theta(n)$ space.

Proof. The algorithm maintains an array of the digits of the fractional part in base $\boldsymbol{b}$ and does exactly what we did in (2). Lemma 1 shows that this $\texttt{q}$ in Step 4 will never be greater than $9$. So all we have to do now is show that $n + 3$ digits in $\boldsymbol{b}$ is enough precision to get $n$ digits in base 10.

First we will prove that for large enough $n$, $n$ digits of precision is enough. So we need to show that $e - e_n < 1/10^n$ where $e_n = (2;1, 1, 1, \dots, 1)$:

Stirling's approximation says that

so

By taking the natural logarithm of

we get

and taking the exponential we get

So $e - e_n < 1/10^n$ for every $n \geq \lceil 10e \rceil = 28$. This means that $n$ digits in $\boldsymbol{b}$ is enough for $n$ digits in decimal when $n \geq 28$, and we can manually verify that $n + 3$ is enough for $n < 28$. We have to keep in mind that we are counting the integer part, which is already in decimal so we just print it directly in the algorithm. This is why we use an array of length $n + 2$ to store the digits in base $\boldsymbol{b}$.

The fact that the algorithm uses $\Theta(n)$ space is trivial given that we always use an array of length $n + 2$. Step 2 takes $\Theta(n)$ time and steps from Step 5 to Step 16 take $\Theta(n^2)$ time so the total running time is $\Theta(n^2)$.

Here is an implementation of this algorithm in C and Haskell:

#include <stdio.h>
#include <stdlib.h>

void generate_digits_of_e(int n)
{
        int q = 2;
        int len = n + 2;
        int *a = malloc(len * sizeof(int));
        if (!a) {
                fprintf(stderr, "Out of memory.");
                exit(1);
        }
        for (int i = 0; i < len; ++i)
                a[i] = 10;
        printf("%d.", q);
        for (int i = 0; i < n - 1; ++i) {
                for (int j = len - 1; j >= 0; --j) {
                        const int d = a[j];
                        a[j] = d % (j + 2);
                        q = d / (j + 2);
                        a[j - 1] += q;
                }
                printf("%d", q);
                fflush(stdout);
                for (int j = 0; j < len; ++j)
                        a[j] *= 10;
        }
        free(a);
}

digits n = 2 : digits' (n - 1) (replicate (n + 2) 1)
digits' n a@(h:t) = let (qs, rs) = unzip $ scanl f (f (0, 0) (h * 10, initB)) t'
                    in last qs : if n == 1 then [] else digits' (n - 1) rs
  where initB = length a + 1
        t' = zip (map (*10) t) [initB - 1, initB - 2..]
        f (q, _) (n, b) = divMod (n + q) b

This algorithm was first published by A. H. J. Sale in his 1967 paper “The Calculation of $e$ to Many Significant Digits”.

The Algorithm for $\pi$

In 1995 Stanley Rabinowitz and Stan Wagon published a similar algorithm for $\pi$. We defined a base $\boldsymbol{b} = (2, 3, 4, \dots)$ for representing $e$ and constructed the algorithm for $e$ by converting this base to decimal, likewise we will use the base $\boldsymbol{c} = (3, 5/2, 7/3, \dots)$ based on (3) to represent $\pi$ as $(2, 2, 2\dots)_{\boldsymbol{c}}$. Base $\boldsymbol{c}$ is, unfortunately, not as nice as $\boldsymbol{b}$. For example, this base with maximal digits and no integer part (meaning every digit $d_i$ in the fractional part is one less than the denominator of $b_i^{-1}$ ) is the number

Where $x!!$ is the double factorial, see this Wikipedia article for what it means (it is not just factorial applied twice). Using the identity

we get

Here I used the fact that

which can be proven using integration by parts, and

This shows that numbers of the form $(0 ; a_1, a_2, a_3, \dots)_{\boldsymbol{c}}$ are in the range $[0, 2]$ and, in contrast with base $\boldsymbol{b}$, the integer part of $(a_0 ; a_1, a_2, a_3, \dots)$ may be $a_0 + 1$ with limited precision and we may get an integer in the range $[10, 19]$ after the multiply-reduce-carry stage of the algorithm.

This is not hard to fix. The algorithm is almost the same, but instead of printing the digit right away we hold it for as long as the next digits are 9's. If we then get an integer greater than 9 we add 1 to the pre-digit we first got and print it, then we print 0's instead of 9's, otherwise we print the pre-digit, all the 9's and set the new pre-digit to the integer modulo $10$. There is nothing complicated going on here, just multiplication with carry. Of course, we only print as many digits as needed, so we may not print all the 9's or all the 0's.

So this algorithm has a very serious flaw compared to the algorithm for Euler's number: at the beginning we need to specify a number of digits to be printed, but we may reach a sequence of 9's that goes beyond this limit, in which case the algorithm would terminate without fixing the pre-digit, possibly printing one minus the correct digit, and a bunch of 9's. We can't fix this by simply allocating more memory for the array and continuing as if we started with more precision since this would mess up the array's state. It is also not possible (or is it?) to predict in advance where and how long the sequences of 9's are, coming up with such a method could even lead to a proof that $\pi$ is or is not normal, an open problem in mathematics. So if you get a sequence of 9's, you need to call the algorithm again with more digits to be sure the digits are correct.

That being said, the authors of the original paper say that we need $\lfloor 10n/3 \rfloor$ base $\boldsymbol{c}$ digits but, as you can check with $n = 1$, this is wrong and we actually need $\lfloor 10n/3 \rfloor + 1$ digits. Here is the code in C:

#include <stdio.h>
#include <stdlib.h>

void generate_digits_of_pi(int n)
{
        const int len = (10 * n) / 3 + 1;
        int *const a = malloc(len * sizeof(len));
        if (!a) {
                fprintf(stderr, "Out of memory.");
                exit(1);
        }
        for (int i = 0; i < len; ++i)
                a[i] = 2;
  
        int nines = 0;
        int predigit = -1;
        while (n >= 0) {
                for (int j = 0; j < len; ++j)
                        a[j] *= 10;
                int q;
                for (int j = len - 1; j >= 1; --j) {
                        const int d = a[j];
                        a[j] = d % (2 * j + 1);
                        q = d / (2 * j + 1);
                        a[j - 1] += q * j;
                }
                // Remember that the integer part of base c uses base 10.
                const int d = a[0];
                a[0] = d % 10;
                q = d / 10;
                if (q == 9) {
                        ++nines;
                } else {
                        // This prints predigit + 1 if the quotient is greater
                        // than 9 and predigit otherwise.
                        if (predigit >= 0)
                                printf("%d", predigit + q / 10);
                        if (n < nines)
                                nines = n;
                        // Decrease n by the number of digits we are about to
                        // print plus the digit we printed.
                        n -= nines + 1;
                        while (--nines >= 0)
                                printf("%d", q == 10 ? 0 : 9);
                        nines = 0;
                        predigit = q == 10 ? 0 : q;
                }
                fflush(stdout);
        }
       free(a);
}

I'm not going to go over the correctness of this algorithm since the idea is basically the same as the algorithm for $e$, if you are interested check the paper by Stanley Rabinowitz and Stan Wagon. Jeremy Gibbons published an algorithm that doesn't have the problems I mentioned, in the paper “Unbounded Spigot Algorithms for the Digits of Pi”. The purpose of the last sections was to present the idea of spigot algorithms and how they can be constructed, if you really are interested in this kind of thing I recommend reading this paper, as I'm about to move on to something a lot more interesting.

Digit Extraction Algorithms

In 1995 Simon Plouffe discovered the remarkable formula

This is the BBP formula, named after David Bailey, Peter Borwein and Plouffe. This formula is special because we can use it to construct a spigot algorithm that is a lot more interesting than the previous algorithm: not only can we generate $\pi$ digit by digit, but we can also calculate the nth hexadecimal digit of $\pi$ without calculating the preceding digits!

To see how this is possible first let

and notice that the digit of $\pi$ in position $(n - 2)$ from 1 (counting the integer part) is the first digit in the fractional part of

where $\{x\} = x - \lfloor x \rfloor$ means the fractional part of $\{x\}$. We are just shifting digits and removing the integer part. From this definition of fractional part we can see that if $n$ is an integer then

which gives us

We can also show that $\{\{a\} - \{b\}\} = \{a - b\}$ the same way. Now we can write

Let's isolate $\{16^n S_i\}$:

The second sum didn't have an integer part in the first place so we don't need to take the fractional part. To see this, just consider the geometric series

So only the first sum has terms with integer part. Since we are taking the fractional part we can remove the integer part of every term by reducing the numerator modulo the denominator (the remainder of $a/b$ divided by $b$ is equal to the fractional part of $a/b$ ). So we now have

Putting everything together:

Of course, the digit we want is $\lfloor 16 \{16^n\pi\} \rfloor$, and now we have a great way to calculate this. The right sum of $\Sigma_i$ is very small and for large $n$ it quickly gets smaller than our floating-point precision of choice can represent. We can also use binary modulo exponentiation to compute $16^{n - k} \bmod (8k + i)$ efficiently, this is the full source code:

#include <math.h>
#include <stdint.h>
#include <stdio.h>
    
/*
 * Returns 16 raised to n, modulo m.
 */
static double expm2(long n, double m)
{
        if (m == 1)
                return 0;
        double result = 1.;
        long double base = fmod(16, m);
        while (n > 0) {
                if (n % 2 == 1)
                        result = fmod(result * base, m);
                n /= 2;
                base = fmod(base * base, m);
        }
        return result;
}
    
/*
 * We need the modulo to be always non-zero since although the digits in the
 * fractional part are all correct the ones in the integer part are not, which
 * means we might end up with a negative number in the last stage of the
 * algorithm where we use modulo 1 to take the fractional part.
 */
static double mod(double x, long m)
{
        x = fmod(x, (double)m);
        return x < 0 ? x + m : x;
}
   
static double sigma(long i, long n)
{
        double s = 0;
    
#pragma omp parallel for reduction(+:s)
        for (long k = 0; k <= n; ++k) {
                const double m = 8 * k + i;
                const double a = expm2(n - k, m) / m;
                s += mod(a, 1);
        }
    
        for (long k = n + 1; k < 20; ++k)
            s += pow(16, n - k) / (8 * k + i);
   
        return mod(s, 1);
}
    
int get_digit(long n) {
        n -= 2;
        double acc = 0;
        acc += 4 * sigma(1, n);
        acc -= 2 * sigma(4, n);
        acc -= sigma(5, n);
        acc -= sigma(6, n);
        return (int)floor(16 * mod(acc, 1));
}

Since this algorithm is embarrassingly parallel I took the liberty of adding an OpenMP pragma to speed things up. The amount of precision required depends on how far the digit is, at some point using 64-bit floating point numbers will not be enough and you would have to use a library like MPFR with some numerical analysis to estimate how much precision you need, but a 64-bit float was enough for me to calculate the 10,000,000th digit (it is $A$!) in 2 seconds. Also, by computing the digit we want we inevitably compute the next few digits so with enough floating-point precision and terms in the sums we can get not only the nth digit but also $k$ digits from that position by multiplying by $16^k$ instead of just $16$. It is important to understand that the reason why this algorithm works is because of the division by $1/16^k$ in the BBP formula, which allows us to multiply by $16^n$ and split each sum into two, this is why this formula is special.

In 1997 Fabrice Bellard discovered the formula

which allows to calculate the hexadecimal or binary digits of $\pi$ 43% faster than the BBP formula. Using 20 computers in parallel Bellard calculated 152 binary digits of $\pi$ starting from the trillionth position. Colin Percival calculated the quadrillionth binary digit $0$ in his PiHex project using a distributed network of 1,734 computers in 56 different countries.

There are many known BBP-type formulas for other constants:

In the $\arctan$ one $b$ is an integer. In general, this kind of formula has the form

where $p$ and $q$ are polynomials with integer coefficients and $b$ is an integer greater than one (this is also the base for which we can construct a digit extraction algorithm).