Pi calculator using Bailey-Borwein-Plouffe formula

pi.py

# PI calculator
#
# - Optimized for 8-bit microcontroller
# - To be translated to C-language
# - Tested up to 500 digits
#
# Original S function taken from
# http://pythonfiddle.com/bailey-borwein-plouffe-formula/
# Heavily modified afterwards

counter = 0 # total addition counter

def add(a, b):
    global counter
    counter += 1
    return a + b

# subtraction
def sub(a, b):
    return add(a, -b)

# integer multiplication
# b can be float
# a cannot be under 1
# a cannot be float
def imul(a, b):
    result = 0
    i = 0
    while i < a:
        result = add(result, b)
        i += 1
    return result

# negative integer modulo
def imod(a, b):
    a = -a
    while a > 0: 
        a = sub(a, b)
    return -a

# power
def pwr(a, b):
    if b == 0: return 1
    result = a
    i = 1
    while i < b:
        result = imul(a, result)
        i += 1
    return result

# integer division
def idiv(a, b):
    quotient = 0
    while a >= b:
        a = sub(a, b)
        quotient = add(quotient, 1)
    return quotient

# Force python scientific notation to good old normal.
# Not to be implemented in C-code obviously.
def StoN(a):
    v = str(a).split('e')
    if len(v) < 2:
        return "".join(str(a).replace('.', ''))
    coef = float(v[0])
    exp = int(v[1])
    result = ""
    if exp > 0:
        result += str(coef).replace('.', '')
        result += "".join(['0' for i in range(0, abs(exp - len(str(coef).split('.')[1])))])
    elif exp < 0:
        result += "0"
        result += "".join(['0' for i in range(0, abs(exp) - 1)])
        result += str(coef).replace('.', '')
    return result

def longdiv(a, b):
    b = int(abs(b))
    length = len(str(int(a)))
    c = StoN(a)    
    c += "0" * (24 + length)
    result = str(idiv(int(c[0]), b))
    rem = c[0]
    i = 0
    while i < 24 + length:
        rem = str(sub(int(rem), imul(int(result[i]), b)))
        rem += c[i + 1]
        result += str(idiv(int(rem), b))
        i += 1
    return float("." + result[length:])

def S(j, n):
    # Left sum
    s = 0
    k = 0
    while k <= n:
        r = add((sub(n, k) << 3), j)
        s = add(s, longdiv(pwr(16, k), r))
        k += 1
    # Right sum
    t = 0
    while 1:
        r = add((k << 3), j)
        nt = add(t, longdiv(longdiv(1, pwr(16, sub(k, n))), r))
        # Iterate until t no longer changes
        if t == nt:
            break
        t = nt
        k += 1
    return add(s, t)

def main(n):
    a = imul(4, S(1, n))
    b = imul(2, S(4, n))
    c = S(5, n)
    d = S(6, n)
    x = sub(a, b)
    x = sub(x, c)
    x = sub(x, d)
    x = imod(x, 1)
    print("position = " + str(n))
    print("fraction = " + str(x))
    print("hex digits = %014x" % int(x * 16**14))

print("Processing PI hex digits...")
main(3)
main(500)
print("Ready. This many additions have been made: " + str(counter))