fast_div.py

from mpmath import mp, mpf
from fractions import Fraction
from gmpy2 import bit_scan1 as ctz

mp.dps = 999
print(mp)


def ceil(x):
    return int(mp.ceil(x))


def log2(x):
    return mp.log(x, 2)


def pow2(x):
    assert isinstance(x, int)
    return 2 ** x


def choose(N, d, prec):
    l = ceil(log2(d))
    sh_post = l

    # m_low = floor( 2^(N+l) / d)
    m_low = pow2(N+l) // d

    # m_high = floor( (2^(N+l) + 2^(N+l-prec))/d )
    m_high = (pow2(N+l) + pow2(N+l-prec)) // d

    while m_low//2 < m_high//2 and sh_post > 0:
        m_low = m_low//2
        m_high = m_high//2
        sh_post = sh_post-1
    return (m_high, sh_post, l)


def recip(N, d, fast=True):
    if fast:
        sh_pre = 0
        (m, sh_post, l) = choose(N, d, N)
        if m >= pow2(N) and d % 2 == 0:
            e = ctz(d)
            d_odd = d >> e
            sh_pre = e
            (m, sh_post, _) = choose(N, d_odd, N-e)
        if d == pow2(l):
            return (l, 0, 1, "power of two")
        elif m >= pow2(N):
            assert sh_pre == 0
            return (1, sh_post-1, m - 2**N, "large")
        else:
            return (sh_pre, sh_post, m, "default")
    else:
        l = ceil(log2(d))
        m = (pow2(N+l) // d) - pow2(N) + 1
        return (1, l-1, m)


assert recip(32, 10) == (0, 3, 3435973837, "default")

list = []
for i in range(0, 20):
    v = recip(64, 10**i)
    list.append(v)

for i, v in enumerate(list):
    print(f"{v}, // 10^{i}")