What the proof goes, you did most of the work already. Consider the following:
l = (n.bit_length() - 1)//4
a = isqrt(n >> 2*l)
a = ((a << l) + n//(a << l))//2
return a - (a*a > n)
This computes the square root of the (possibly longer) upper half, applies one Heron's step and a single correction. I think it is functionally equal to what you wrote. Those zeros don't contribute to the quotient so we could instead write:
a = ((a << l) + (n >> l)//a)//2
The problem is that the 3n/n division in the step `(a + n//a)//2` basically recomputes the upper half we already know and so we want to avoid it: instead of 3n/n giving 2n quotient, we want 2n/n giving 1n quotient. If the upper half is correct, the lower half to be taken care of is `n - a**2`:
a = (a << l) + ((n - (a << l)**2) >> l)//a//2
And there is no need to square the zeros either:
a = (a << l) + ((n - (a**2 << 2*l) >> l)//a//2
So I *think* it should be correct, the only thing I'm not sure about is whether the final correction in the original `isqrt` is needed. Perhaps the automated proof of yours could give an answer? |