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<