Reducing integer division to integer multiplication.

division.py

# Reducing integer division to integer multiplication.
#
# Let M(n) denote the time it takes to multiply two n-bit
# integers. Say we want to find the quotient of the division a/b for
# two integers a and b. The simplest reduction is a binary search for
# the answer in the range [0, a], leading to a O(M(n)*n) algorithm for
# division. How can we do better?
#
# Here's an idea. Suppose we find 1 / b = 0.uvwxyz... accurately up to
# some precision (in, say, base 10) then we can just multiply
# sufficiently many digits of the inverse of `b` with `a` to get the
# answer. We can find inverses using Newton's method, by setting up
# our problem as finding a root of f(y) = b - 10^k / y, where the
# factor 10^k shifts unity so that the root is located at 10^k / b,
# for some k to be determined. The successive approximants to the root
# are given by
#
#  y_{n+1} = y_n * (2 - y_n * b / 10^k) = 2 * y_n - y_n * y_n * b / 10^k.
#
# Notice that each iteration involves three integer multiplications
# (division in the chosen base is free). As the speed of convergence
# of Newton's method to the root is quadratic (doubling precision at
# each step), we are going to hit the value floor(10^k / b) in O(\log
# k) iterations. For k roughly \log n chosen such that a < 10^k, we
# then get a division algorithm that should run in O(\log n * M(n)).

from random import randint

# global counter for the number of multiplication steps used.
muls = 0

# Multiplication algorithm.
def M(a, b):
    global muls
    muls += 1
    return a*b

# Return the ceiling of a / b
def intceil(a, b):
    return a // b if (a % b == 0) else a // b + 1

# Returns an integer y such that y*y / (1 + y) < p / b <= y, enclosing
# the quotient of the division in an interval of size y / (1 + y),
# which is less than 1. In other words, y = ceil(p / b).
def inv(b, p):
    print('Computing inverse ceil(%d / %d).' % (p, b))
    y_old, y_new = 0, 1
    iters = 0
    while y_old != y_new:
        y_old = y_new
        y_new = M(2, y_new) - M(y_new, M(y_new, b)) // p
        # Ensure that successive approximants are increasing in value.
        y_new = max(y_new, y_old)
        iters += 1
    # Assert the fixed point.
    assert y_new == intceil(p, b)
    print('Determined inverse to be %d in %d Newton iterations.' % (y_new, iters))
    return y_new

# The division algorithm that returns a pair (q, r) such that:
# a = q*b + r and 0 <= abs(r) < abs(b).
def D(a, b):
    print('Dividing %d by %d.' % (a, b))
    global muls
    muls = 0
    if b == 0: raise ZeroDivisionError
    neg_a = True if a < 0 else False
    neg_b = True if b < 0 else False
    a = -a if a < 0 else a
    b = -b if b < 0 else b
    q = D_san(a, b) if a >=b else 0
    r = a - M(q, b)
    assert(0 <= r < b)
    if neg_a and neg_b:
        r *= -1
    elif neg_a:
        q *= -1
        r *= -1
    elif neg_b:
        q *= -1
    print('Determined that %d = %d * %d + %d using %d multiplications.' %
          (a * (-1 if neg_a else 1), q, b * (-1 if neg_b else 1), r, muls))
    return (q, r)

# Assume a >= b > 0. We want to compute floor(a / b) but we only have
# floor(a * ceil(p / b) / p). We have the inequalities a/b <= a *
# ceil(p / b) / p < a/b + 1 and therefore, floor(a / b) must be either
# floor((a * ceil(p / b)) / p) or floor((a * ceil(p / b)) / p) - 1.
def D_san(a, b, base=10):
    print('Computing in base %d' % (base,))
    k, p = 1, base
    while a >= p:
        p *= base
        k += 1
    q = M(a, inv(b, p)) // p
    # Find the correct quotient using one more multiplication.
    q = q - 1 if (a - M(q, b) < 0) else q
    return q

if __name__ == '__main__':
    a, b = 8997126187981209870390813274098, 187961918702646
    print('Computed: %s (versus %s)' % (D(a, b), (a // b, a % b)))

    # some tests
    N1, N2 = 10 ** 50, 10 ** 25
    for i in range(100):
        a, b = randint(-N1, N1), randint(-N2, N2)
        if b == 0: continue
        q, r = D(a, b)
        assert(a == q*b + r)
        assert(0 <= abs(r) < abs(b))
        print('%d / %d = %s.\n' % (a, b, (q, r)))