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
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]
, or 
, which do not necessarily produce the same result. That is, double-precision addition is not associative:
might be computed in, say, 80-bit precision, before clamping the final sum to 64-bit double precision. This has bitten me in 
Possibly Wrong 