Pi with the BBP formula (Python)

From LiteratePrograms

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

where

Basic implementation

Calculating the n-th digit of π is equivalent to calculating the first digit in the fractional part of 16ⁿπ. 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>>=
def S(j, n):
    # Left sum
    s = 0.0
    k = 0
    while k <= n:
        r = 8*k+j
        s = (s + float(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 16⁻ⁿ.

Test output

In hexadecimal, π begins 3.243f6a8885a308d31319...

The correct digits used for reference are taken from 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⁶:

Calculated: 26c65e52cb3188
Correct:    26c65e52cb4593

Starting with n = 10⁷:

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)
MASK = M - 1
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 = 10⁸:

Calculated: ecb840e2188552
Correct:    ecb840e21926ec

Sequential generation

The following result can be derived from the BBP formula: if one sets x₀ = 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>>=
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():
    for n, p in enumerate(pi()):
        print "%x" % p
        if n % 1000 == 0:
            print "\n\n%i\n\n" % n

References

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

Download code

Retrieved from "http://literateprograms.org/Pi_with_the_BBP_formula_%28Python%29"

Categories: Pi with the BBP formula | Programming language:Python