sage_bbp.ipynb

Kernel: SageMath 7.2.rc2

SageMath version of BBP

Experiments with BBP algorithm in Sage.

First version

def ratio(d, k, j):
    return power_mod(16, d - k, 8*k + j) / (8.0*k + j)

def frac_part(x):
    return x - floor(x)

def sigma(d, j):
    r = 0.0
    for k in xrange(0, d + 1):
        r = frac_part(r + ratio(d, k, j))

    k = d + 1
    u = 1.0 / (16^(k-d) * (8*k + j))

    while u > RR.epsilon():
        r += u
        k += 1
        u = 1.0 / (16^(k-d) * (8*k + j))

    return frac_part(r)

def digpi(d):
    r = 4*sigma(d, 1) - 2*sigma(d, 4) - sigma(d, 5) - sigma(d, 6)
    return floor(16*frac_part(r))

The $10\,001$-th hexadecimal digit of $\pi$.

digpi(10000)

8

Let us measure the elapsed time.

%%timeit
digpi(10000)

1 loop, best of 3: 2.3 s per loop

Is it possible to make this run faster?

Second version (changes)

def ratio(d, k, j):
    return power_mod(16, d - k, 8*k + j) / (8.0*k + j)

def frac_part(x):
    return x - floor(x)

def sigma(d, j):
    r = 0.
    for k in xrange(0, d + 1):
        r += ratio(d, k, j)
        if r > 1.:
            r = frac_part(r)

    kk  = 8*(d + 1) + j
    u   = 1. / 16
    rat = u / kk

    while rat > RR.epsilon():
        r   += rat
        kk  += 8
        u   /= 16
        rat  = u / k

    return frac_part(r)

def digpi(d):
    r = 4*sigma(d, 1) - 2*sigma(d, 4) - sigma(d, 5) - sigma(d, 6)
    return floor(16*frac_part(r))

Speed test.

%%timeit
digpi(10000)

1 loop, best of 3: 2.26 s per loop

Cython version

Load and compile the Cython code (inspect the file).

load("bbp.spyx")

Does it work?

cdigpi(999999)

2

How fast is it?

%%timeit
cdigpi(10000)

10 loops, best of 3: 23.9 ms per loop

Compare this with the previous result! This is faster by factor 10.

Does it really compute what we think it computes?

sum(cdigpi(k) * 16.^(-k-1) for k in xrange(10))

0.141592653589214