# Pi with the BBP formula (Python)

The BBP formula allows us to extract isolated hexadecimal digits of π. The formula is

$\pi = 4S(1) - 2S(4) - S(5) - S(6)\,$

where

$S(j) = \sum_{k=0}^\infty \frac{1}{16^k} \left( \frac{1}{8k+j} \right).$

## Basic implementation

Calculating the n-th digit of π is equivalent to calculating the first digit in the fractional part of 16nπ. Normally, we would have to perform the entire summation with at least n digits of precision and afterwards extract the fractional part — the efficiency of the BBP algorithm is due to the fact that the entire computation can be performed on the fractional part, and we can even use ordinary floating-point arithmetic.

If we use {} to denote fractional parts, the summation can be rewritten as follows:

$\left\{16^n S(j)\right\} = \left\{\sum_{k=0}^{n} \frac{16^{n-k} \mod 8k+j}{8k+j} + \sum_{k = n + 1}^{\infty} \frac{16^{n-k}}{8k+j}\right\}.$

The key step is the modulo operation in the left-hand sum, which is justified since it simply throws away the useless integer part of the term. Integer exponentiation modulo a small integer can be performed in time O(log p) where p is the exponent, using Python's built in pow() function, so it is easy to see that only O(n log n) time will be required for the entire calculation. The right-hand sum converges rapidly and is trivial to calculate.

The following function implements the S summation as written above:

<<pibbp.py>>=
def S(j, n):
# Left sum
s = 0.0
k = 0
while k <= n:
r = 8*k+j
s = (s + pow(16, n-k, r) / r) % 1.0
k += 1
# Right sum
t = 0.0
k = n + 1
while 1:
newt = t + pow(16, n-k) / (8*k+j)
# Iterate until t no longer changes
if t == newt:
break
else:
t = newt
k += 1
return s + t


Note the reduction to a fractional part (using % 1.0), which must be performed each time a new term is added to the first sum to minimize rounding error.

The main function evaluates the sum and ouputs the value as a string of hexadecimal digits. Converting to hexadecimal can be done by multiplying with a power of 16, rounding to an integer, and converting the integer to hexadecimal. The summation of S is performed with double precision using a 53-bit mantissa, corresponding to roughly 13 hex digits, so taking 14 hex digits will extract all information from the number (along with some accumulated rounding error!).

<<pibbp.py>>=
def pi(n):
n -= 1
x = (4*S(1, n) - 2*S(4, n) - S(5, n) - S(6, n)) % 1.0
return "%014x" % int(x * 16**14)


The decrement of n serves to index the digit immediately after the radix point with n = 1, rather than with n = 0. This is an arbitrary choice; both conventions are in use elsewhere. The advantage of the convention chosen here is that the '3' to the left of the radix point can be indexed with a nonnegative number, and that the weight for the n-th digit simply becomes 16n.

## Test output

The correct digits used for reference are taken from [1] and from Borwein & Bailey.

Starting with n = 0:

Calculated: 3243f6a8885a36
Correct:    3243f6a8885a30


Starting with n = 1:

Calculated: 243f6a8885a300
Correct:    243f6a8885a308


Starting with n = 10000:

Calculated: 68ac8fcfb80310
Correct:    68ac8fcfb8016c


Starting with n = 106:

Calculated: 26c65e52cb3188
Correct:    26c65e52cb4593


Starting with n = 107:

Calculated: 17af5863f0f5b0
Correct:    17af5863efed8d


The rounding error is significant but not disastrous; it should be possible to extract a correct digit or two around the billionth. Timings on a 1400 MHz Athlon suggest that such a computation would require about two days of CPU time with the present implementation.

## Using integer arithmetic

Instead of using floating-point arithmetic, one can perform fixed-point arithmetic on integers. This way a custom precision level can be used, enabling arbitrarily large computations.

<<pibbp2.py>>=
D = 14        # number of digits of working precision
M = 16 ** D
SHIFT = (4 * D)

def S(j, n):
# Left sum
s = 0
k = 0
while k <= n:
r = 8*k+j
s = (s + (pow(16,n-k,r)<<SHIFT)//r) & MASK
k += 1
# Right sum
t = 0
k = n + 1
while 1:
xp = int(16**(n-k) * M)
newt = t + xp // (8*k+j)
# Iterate until t no longer changes
if t == newt:
break
else:
t = newt
k += 1
return s + t

def pi(n):
n -= 1
x = (4*S(1, n) - 2*S(4, n) - S(5, n) - S(6, n)) & MASK
return "%014x" % x


Output starting with n = 108:

Calculated: ecb840e2188552
Correct:    ecb840e21926ec


## Sequential generation

The following result can be derived from the BBP formula: if one sets x0 = 0 and for n = 1, 2, 3, ... computes

$x_n = \left\{ 16x_{n-1} + \frac{120n^2 - 89n + 16}{512n^4-1024n^3+712n^2-206n+21}\right\},$

the n-th hexadecimal digit of π − 3 is (according to experimental evidence) given by

$\lfloor 16 x_n \rfloor.$

We may use this to write a generator that yields hexadecimal digits of π one by one. If we want to generate lots of digits, this method should be faster than repeatedly calculating isolated digits the "ordinary" way with the BBP formula. Rational arithmetic is required, so we use the mpq (rational) and mpz (integer) types from GMPY.

<<sequential_pi.py>>=
from gmpy import mpq, mpz

def mod1(x):
return x-mpz(x)

def pi():
x = 0
n = 1
while 1:
p = mpq((120*n-89)*n+16, (((512*n-1024)*n+712)*n-206)*n+21)
x = mod1(16*x + p)
n += 1
yield int(16*x)


This pi function can for instance be used as follows:

<<sequential_pi.py>>=
def allpi():
from sys import stdout
p = pi()
n = 1
while 1:
stdout.write("%x" % p.next())
if n % 1000 == 0:
print "\n\n%i\n\n" % n
n += 1


## References

• Jonathan Borwein & David Bailey, Mathematics by Experiment - Plausible Reasoning in the 21st Century, A K Peters 2003, 118-125

hijacker
hijacker
hijacker
hijacker