bbp.spyx

##
#
# BBP: Cython
#

cimport cython
from libc.math cimport floor
from libc.float cimport DBL_EPSILON

cdef unsigned long cpower_mod(unsigned long a, unsigned long b, unsigned long n):
    cdef unsigned long r = 1
    while b > 0:
        if b % 2 == 1:
            r = r * a % n
        b /= 2
        a = a * a % n
    return r

@cython.cdivision(True)
cdef double cratio(unsigned long d, unsigned long k, unsigned long j):
    return <double> cpower_mod(16, d - k, 8*k + j) / (8*k + j)

cdef double cfrac_part(double x):
    return x - floor(x)

@cython.cdivision(True)
cdef csigma(unsigned long d, unsigned long j):
    cdef double r = 0.
    cdef unsigned long k

    for k in range(0, d + 1):
        r += cratio(d, k, j)
        if r > 1.:
            r = cfrac_part(r)

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

    while rat > DBL_EPSILON:
        r   += rat
        kk  += 8
        u   /= 16
        rat  = u / kk

    return cfrac_part(r)

def cdigpi(unsigned long d):
    cdef double r = 4*csigma(d, 1) - 2*csigma(d, 4) - csigma(d, 5) - csigma(d, 6)
    return <int> floor(16*cfrac_part(r))