OK, here's the last version I had. Preconditions are that d > 0, n > 0, and n % d == 0.

This version tries to use the narrowest possible integers on each step. The lowermost `good_bits` of dinv at the start of the loop are correct already.

Taking out all the modular stuff, the body of the loop boils down to just

    dinv *= 2 - dinv * d

For insight, if

    dinv * d = 1 + k*2**i

for some k and i (IOW, if dinv * d = 1 modulo 2**i), then

    2 - dinv * d = 1 - k*2**i

and so dinv times that equals 1 - k**2 * 2**(2*i). Or, IOW, the next value of dinv is such that d * dinv = 1 modulo 2**(2*i) - it's good to twice as many bits.

    def ediv(n, d):
        assert d

        def makemask(n):
            return (1 << n) - 1

        if d & 1 == 0:
            ntz = (d & -d).bit_length() - 1
            n >>= ntz
            d >>= ntz
        bits_needed = n.bit_length() - d.bit_length() + 1
        good_bits = 3
        dinv = d & 7
        while good_bits < bits_needed:
            twice = min(2 * good_bits, bits_needed)
            twomask = makemask(twice)
            fac2 = dinv * (d & twomask)
            fac2 &= twomask
            fac2 = (2 - fac2) & twomask
            dinv = (dinv * fac2) & twomask
            good_bits = twice
        goodmask = makemask(bits_needed)
        return ((dinv & goodmask) * (n & goodmask)) & goodmask