BourweinPlouffePi_nth.py

# Python program with numba's njit calculates only nth hex
# Computing nth hex digit of Pi with BBP formula.
# Bailey-Borwein-Plouffe_formula for Pi
# Uses numba compiler & njit directive for speed
# Juhani Kaukoranta 2.8.2018
# muokattu 4.8.2019 poistamalla 10 hex numeron tulostus

import math
import time
#from numba import njit

 # Convert Pi to Hex-strings
 # (float number's decimals to hex digits)

def convHex(x) :
    hx = "0123456789ABCDEF"
    y = math.fabs(x)
    chx =''
    for i in range(0,16):
        y = 16.0 * (y - math.floor(y))
        chx += hx[int(y)]
    return chx

 # Compute Modular Exponentiation
@njit
def compModExp(b, e, m):
    if (e == 0):
        return 1
    ans = compModExp(b, e // 2, m)
    ans = (ans * ans) % m
    if ((e % 2) == 1):
        ans = (ans * b) % m
    return ans

# Compute Sums
@njit
def Sums(j):
    s = 0.0
#     t              // Each term of right summation
#     r              // Denominator
#     k              // Loop index
# EPS = Loop-exit accuration of the right summation

    EPS = 1.0e-17
    # Left Sum (0 ... n)
    k = 0
    for k in range(0,n+1):
        r = 8 * k + j
        t = compModExp(16, n - k, r)
        t /= r
        s += t - int(t)
        s -= int(s)

    #Right sum (n + 1 ...)
    k = n+1
    while (1) :
        r = 8 * k + j
        t = 16**(n - k)
        t /= r
        if (t < EPS):
            break
        s += t
        s -= int(s)
        k += 1

    return s # left + right

# paaohjelma kysyy halutun hex-muotoisen desimaalin järjestysnumeroa
# tulostaa sen ja laskennassa käytetyn desimaalimuotoisen fraktion

n = int(input("kuinka mones Piin desimaali halutaan laskea?"))
print("**** PI Computation ( hex digit number: ", n,", waiting..." )
# Compute PI
n = n - 1
t0 = time.perf_counter() # time start
pi = 4 * Sums(1) - 2 * Sums(4) - Sums(5) - Sums(6)
pi = pi - int(pi) + 1
pi_hex = convHex(pi)
print("Hex digit: ", pi_hex[0])
print("FRACTION  : ",   pi)
t1 = time.perf_counter() # time end
print("Aikaa kului ", t1 - t0," sekuntia")