from math import log

def rshift(x, n):
    if n >= 0: return x >> n
    else:      return x << (-n)

def lshift(x, n):
    if n >= 0: return x << n
    else:      return x >> (-n)

def giant_steps(start, target):
    """Return a list of integers ~= [~2*start, ..., target/2, target],
    describing suitable precision steps for quadratically convergent
    algorithms (e.g. Newton's method).

        >>> giant_steps(53, 1000)
        [64, 126, 251, 501, 1000]

    """
    L = [target]
    while L[-1] > start*2:
        L = L + [L[-1]//2 + 1]
    return L[::-1]

def div_newton(p, q, guard_bits=0):
    """
    Asymptotically fast division of positive integers p and q.
    """
    if not q:
        raise ZeroDivisionError
    if not p:
        return 0
    szp = int(log(p,2))+1
    szq = int(log(q,2))+1
    szr = szp - szq
    initial_prec = 15
    if min(szp, szq, szr) < 2*initial_prec:
        return p//q
    origp = p
    if guard_bits:
        p <<= guard_bits
        szp += guard_bits
        szr += guard_bits
    r = (1 << (2*initial_prec)) // (q >> (szq - initial_prec))
    last_prec = initial_prec
    for prec in giant_steps(initial_prec, szr):
        a = lshift(r, prec-last_prec+1)
        b = rshift(r**2 * rshift(q, szq-prec), 2*last_prec)
        r = a - b
        last_prec = prec
    return ((p >> szq) * r) >> (szr + guard_bits)

def idivmod(p, q):
    """
    Fast divmod for large integers.
    """
    psign = p < 0
    qsign = q < 0
    if psign and qsign:
        r, m = idivmod(-p, -q)
        return r, -m
    elif psign:
        r, m = idivmod(-p, q)
        if m:
            return -r-1, q-m
        return -r, m
    elif qsign:
        r, m = idivmod(p, -q)
        if m:
            return -r-1, q+m
        return -r, m
    r = div_newton(p, q, 30)
    m = p - r*q
    # The worst case would be for div_newton to return complete garbage.
    # In that case the following divmod ensures that the result is
    # correct (it will just be slow). Typically, however, r will be
    # off by 0 or 1 units, so the builtin divmod only needs to perform
    # a single-limb division here.
    r2, m = divmod(m, q)
    return r + r2, m

def idiv(p, q):
    """
    Fast division for large integers.
    """
    return idivmod(p, q)[0]

def check(a, b):
    q1 = divmod(a, b)
    q2 = idivmod(a, b)
    assert q1 == q2

def test_division(size=100000, N=100):
    from random import seed, randint, choice
    print size, N
    for i in xrange(1, N):
        seed(i)
        e1 = randint(10,size)
        e2 = randint(10,size)
        a = randint(0,1<<e1) * choice([1, -1])
        b = randint(1,1<<e2) * choice([1, -1])
        #print i, e1, e2
        check(a, b)
        ab = a*b
        check(ab, b)
        check(ab+1, b)
        check(ab-1, b)

def numeral(n):
    if n < 0:
        return "-" + numeral(-n)
    if not n:
        return "0"
    # Fast enough to do directly
    size = int(log(n,10)+1)
    if size < 250:
        return str(n)
    # Divide in half
    half = (size // 2) + (size & 1)
    A, B = idivmod(n, 10**half)
    ad = numeral(A)
    bd = numeral(B).rjust(half, "0")
    return ad + bd

from timeit import default_timer as clock

def time_division(INT, twiceness=2.0, upto=18):
    fmt = "%10s %12f %12f %12f %8s"
    print "%10s %12s %12s %12s %8s" % ("size", "old time", "new time",
        "faster", "correct")
    print "-"*78
    for i in range(4,upto):
        n = 2**i
        Q = INT(10)**n
        P = INT(10)**int(twiceness*n)
        t1 = clock()
        R1 = divmod(P,Q)
        t2 = clock()
        R2 = idivmod(P,Q)
        t3 = clock()
        size, old_time, new_time = n, t2-t1, t3-t2
        faster, correct = old_time/new_time, R2==R1
        print fmt % (size, old_time, new_time, faster, correct)

def time_str(INT, upto=18):
    fmt = "%10s %12f %12f %12f %8s"
    print "%10s %12s %12s %12s %8s" % ("size", "old time", "new time",
        "faster", "correct")
    print "-"*78
    from random import randint
    for i in range(4,upto):
        n = 2**i
        Q = randint(INT(10)**(n-1), INT(10)**n)
        t1 = clock()
        R1 = str(Q)
        t2 = clock()
        R2 = numeral(Q)
        t3 = clock()
        size, old_time, new_time = n, t2-t1, t3-t2
        faster, correct = old_time/new_time, R2==R1
        print fmt % (size, old_time, new_time, faster, correct)

def testit():
    test_division(size=1000, N=2000)
    test_division(size=10000, N=200)
    test_division(size=100000, N=20)
    test_division(size=1000000, N=2)

if __name__ == "__main__":
    testit()