# Pi with the BBP formula (Python)

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

where

## Contents |

## [edit] Basic implementation

Calculating the *n*-th digit of π is equivalent to calculating the first digit in the fractional part of 16^{n}π. 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:

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>>=defS(j, n): # Left sum s = 0.0 k = 0whilek <= n: r = 8*k+j s = (s + pow(16, n-k, r) / r) % 1.0 k += 1 # Right sum t = 0.0 k = n + 1while1: newt = t + pow(16, n-k) / (8*k+j) # Iterate until t no longer changesift == newt:breakelse: t = newt k += 1returns + 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>>=defpi(n): n -= 1 x = (4*S(1, n) - 2*S(4, n) - S(5, n) - S(6, n)) % 1.0return"%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 16^{−n}.

## [edit] Test output

In hexadecimal, π begins 3.243f6a8885a308d31319...

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* = 10^{6}:

Calculated:26c65e52cb3188 Correct: 26c65e52cb4593

Starting with *n* = 10^{7}:

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.

## [edit] 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) MASK = M - 1defS(j, n): # Left sum s = 0 k = 0whilek <= n: r = 8*k+j s = (s + (pow(16,n-k,r)<<SHIFT)//r) & MASK k += 1 # Right sum t = 0 k = n + 1while1: xp = int(16**(n-k) * M) newt = t + xp // (8*k+j) # Iterate until t no longer changesift == newt:breakelse: t = newt k += 1returns + tdefpi(n): n -= 1 x = (4*S(1, n) - 2*S(4, n) - S(5, n) - S(6, n)) & MASKreturn"%014x" % x

Output starting with *n* = 10^{8}:

Calculated:ecb840e2188552 Correct: ecb840e21926ec

## [edit] Sequential generation

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

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

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>>=fromgmpyimportmpq, mpzdefmod1(x):returnx-mpz(x)defpi(): x = 0 n = 1while1: 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>>=defallpi():fromsysimportstdout p = pi() n = 1while1: stdout.write("%x" % p.next())ifn % 1000 == 0:

## [edit] References

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

Download code |