[qntm]

I had the same idea. The problem is that shifting the bits to the right results in decimal points falling off the end, which is actually a big deal because they add up to something substantial (i.e. more than 1) by the end of the sum.

If N = 15692343, the result should be 5230781. Using "N /= 2", you get:

    0
     + 3923085.75 = 3923085.75
     + 980771.4375 = 4903857.1875
     + 245192.859375 = 5149050.046875
     + 61298.21484375 = 5210348.26171875
     + 15324.5537109375 = 5225672.81542969
     + 3831.13842773438 = 5229503.95385742
     + 957.784606933594 = 5230461.73846436
     + 239.446151733398 = 5230701.18461609
     + 59.8615379333496 = 5230761.04615402
     + 14.9653844833374 = 5230776.01153851
     + 3.74134612083435 = 5230779.75288463
     + 0.935336530208588 = 5230780.68822116
     + 0.233834132552147 = 5230780.92205529
     + 0.0584585331380367 = 5230780.98051382
     + 0.0146146332845092 = 5230780.99512846
     + 0.0036536583211273 = 5230780.99878211
     + 0.000913414580281824 = 5230780.99969553
     + 0.000228353645070456 = 5230780.99992388
    ...

which eventually reaches the correct answer for any desired degree of accuracy. But using "N >>= 2", you get:

    0
     + 3923085 = 3923085
     + 980771 = 4903856
     + 245192 = 5149048
     + 61298 = 5210346
     + 15324 = 5225670
     + 3831 = 5229501
     + 957 = 5230458
     + 239 = 5230697
     + 59 = 5230756
     + 14 = 5230770
     + 3 = 5230773
     + 0 = 5230773

which is out by 8.

[tzs]

This can be fixed by keeping track of some of the spillage. Here's a solution using this approach that works for 32-bit unsigned integers:

    def badd(A, C):
        while C != 0:
            t = A & C
            A = A ^ C
            C = (t << 1) & 0xFFFFFFFF
        return A

    def div3(ah):
        qh = 0
        ql = 0
        al = 0
        while ah != 0:
            al = (al >> 2) & 0x0000FFFF
            al = badd(al, (ah & 0x3) << 14)
            ah = ah >> 2
            qh = badd(qh, ah)
            ql = badd(ql, al)
            if ql & 0xFFFF0000:
                qh = badd(qh, ql >> 16)
                ql = ql & 0xFFFF
        if ql > 0x8000:
            qh = badd(qh, 1)
        return qh