**Introduction**

This post was initially motivated by an interesting recent article by Chris Wellons discussing the Blowfish cipher. The Blowfish cipher’s subkeys are initialized with values containing the first 8336 hexadecimal digits of

So, how do we compute hexadecimal– or decimal, for that matter– digits of

**Pi is easy with POSIX**

First, Chris’s implementation of the Blowfish cipher includes a script to automatically generate the code defining the subkeys. The following two lines do most of the work:

cmd='obase=16; scale=10040; a(1) * 4' pi="$(echo "$cmd" | bc -l | tr -d '\n\\' | tail -c+3)"

This computes base 16 digits of `a(1) * 4`

, or 4 times the arctangent of 1 (i.e.,

How might we do the same thing on Windows? There are plenty of approximation formulas and algorithms for computing digits of

- in any given base,
- one digit at a time “indefinitely,” i.e., without committing to a fixed precision ahead of time,
- with a relatively simple implementation,
- using only arbitrary-precision
*integer*arithmetic (such as is built into Python, or maybe C++ with a library)?

**Bailey-Borwein-Plouffe**

The Bailey-Borwein-Plouffe (BBP) formula seems ready-made for our purpose:

This formula has the nice property that it may be used to efficiently compute the

However, getting the implementation details right can be tricky. For example, this site provides source code and example data containing one million hexadecimal digits of

But suppose that we don’t want to “look ahead,” but instead want to generate *all* hexadecimal digits of

def pi_bbp(): """Conjectured BBP generator of hex digits of pi.""" a, b = 0, 1 k = 0 while True: ak, bk = (120 * k**2 + 151 * k + 47, 512 * k**4 + 1024 * k**3 + 712 * k**2 + 194 * k + 15) a, b = (16 * a * bk + ak * b, b * bk) digit, a = divmod(a, b) yield digit k = k + 1 for digit in pi_bbp(): print('{:x}'.format(digit), end='')

The idea is similar to converting a fraction to a string in a given base: multiply by the base (16 in this case), extract the next digit as the integer part, then repeat with the remaining fractional part. Here `fractions`

module doesn’t significantly improve readability, and is *much* slower than managing the numerators and denominators directly.)

Now for the interesting part: although this implementation appears to behave correctly– at least for the first 500,000 digits where I stopped testing– it isn’t clear to me that it is *always* correct. That is, I don’t see how to prove that this algorithm will continue to generate correct hexadecimal digits of

(Initial thoughts: Since it’s relatively easy to show that each term

Interestingly, Bailey *et. al.* conjecture (see Reference 1 below) a similar algorithm that they have verified out to 10 million hexadecimal digits. The algorithm involves a strangely similar-but-slightly-different approach and formula:

def pi_bbmw(): """Conjectured BBMW generator of hex digits of pi.""" a, b = 0, 1 k = 0 yield 3 while True: k = k + 1 ak, bk = (120 * k**2 - 89 * k + 16, 512 * k**4 - 1024 * k**3 + 712 * k**2 - 206 * k + 21) a, b = (16 * a * bk + ak * b, b * bk) a = a % b yield 16 * a // b

Unfortunately, this algorithm is slower, requiring one more expensive arbitrary-precision division operation per digit than the BBP version.

**Proven algorithms**

Although the above two algorithms are certainly short and sweet, (1) they only work for generating hexadecimal digits (vs. decimal, for example), and (2) we don’t actually know if they are correct. Fortunately, there are other options.

Gibbons (Reference 2) describes an algorithm that is not only proven correct, but works for generating digits of

def pi_gibbons(base=10): """Gibbons spigot generator of digits of pi in given base.""" q, r, t, k, n, l = 1, 0, 1, 1, 3, 3 while True: if 4 * q + r - t < n * t: yield n q, r, t, k, n, l = (base * q, base * (r - n * t), t, k, (base * (3 * q + r)) // t - base * n, l) else: q, r, t, k, n, l = (q * k, (2 * q + r) * l, t * l, k + 1, (q * (7 * k + 2) + r * l) // (t * l), l + 2)

The bad news is that this is by far the slowest algorithm that I investigated, nearly an order of magnitude slower than BBP on my laptop.

The good news is that there is at least one other algorithm, that is not only competitive with BBP in terms of throughput, but is also general enough to easily compute– in any base– not just the digits of

The idea is to express the desired value as a generalized continued fraction:

where in particular

Then digits may be extracted similarly to the BBP algorithm above: iteratively refine the *convergent* (i.e., approximation) of the continued fraction until the integer part doesn’t change; extract this integer part as the next digit, then multiply the remaining fractional part by the base and continue. In Python:

def continued_fraction(a, b, base=10): """Generate digits of continued fraction a(0)+b(1)/(a(1)+b(2)/(...).""" (p0, q0), (p1, q1) = (a(0), 1), (a(1) * a(0) + b(1), a(1)) k = 1 while True: (d0, r0), (d1, r1) = divmod(p0, q0), divmod(p1, q1) if d0 == d1: yield d1 p0, p1 = base * r0, base * r1 else: k = k + 1 x, y = a(k), b(k) (p0, q0), (p1, q1) = (p1, q1), (x * p1 + y * p0, x * q1 + y * q0)

This approach is handy because not only

for digit in continued_fraction(lambda k: 0 if k == 0 else 2 * k - 1, lambda k: 4 if k == 1 else (k - 1)**2, 10): print(digit, end='')

Or to generate digits of the golden ratio

for digit in continued_fraction(lambda k: 1, lambda k: 1, 10): print(digit, end='')

**Consuming blocks of a generator**

Finally, once I got around to actually using the above algorithm to try to reproduce Chris’s original code generation script, I accidentally injected a bug that took some thought to figure out. Recall that the Blowfish cipher has a couple of (sets of) subkeys, each populated with a segment of the sequence of hexadecimal digits of *subsequent* block of digits, do something else, etc.

A simple way to do this in Python is with the built-in zip function, that takes multiple iterables as arguments, and returns a single generator that outputs tuples of elements from each of the inputs… and “truncates” to the length of the shortest input. In this case, to extract a fixed number of digits of `zip`

the “infinite” digit generator together with a range of the desired length.

To more clearly see what happens, let’s simplify the context a bit and just try to print the first 10 *decimal* digits of

digits = continued_fraction(lambda k: 0 if k == 0 else 2 * k - 1, lambda k: 4 if k == 1 else (k - 1)**2, 10) for digit, k in zip(digits, range(5)): print(digit, end='') print() for digit, k in zip(digits, range(5)): print(digit, end='')

This doesn’t work: the resulting output blocks are (31415) and (26535)… but

The problem is that `zip`

evaluates each input iterator in turn, stopping only when one of them is exhausted. In this case, during the 6th iteration of the first loop, we “eat” the 9 from the digit generator *before* we realize that the `range`

iterator is exhausted. When we continue to the second block of 5 digits, we can’t “put back” the 9.

This is easy to fix: just reverse the order of the `zip`

arguments, so the `range`

is exhausted first, *before* eating the extra element of the “real” sequence we’re extracting from.

for k, digit in zip(range(5), digits): print(digit, end='') print() for k, digit in zip(range(5), digits): print(digit, end='')

This works as desired, with output blocks (31415) and (92653).

**References:**

- Bailey, D., Borwein, J., Mattingly, A., Wightwick, G., The Computation of Previously Inaccessible Digits of
and Catalan’s Constant, *Notices of the American Mathematical Society*,**60**(7) 2013, p. 844-854 [PDF] - Gibbons, J., Unbounded Spigot Algorithms for the Digits of Pi,
*American Mathematical Monthly*,**113**(4) April 2006, p. 318-328 [PDF]