Challenge: find the first value where two functions differ

[Ian Kelly]

On Fri, Aug 4, 2017 at 8:51 AM, Steve D'Aprano
<steve+python at pearwood.info> wrote:
> This is a challenge for which I don't have a complete answer, only a partial
> answer.
>
> Here are two functions for calculating the integer square root of a non-negative
> int argument. The first is known to be exact but may be a bit slow:
>
> def isqrt_newton(n):
>     """Integer sqrt using Newton's Method."""
>     if n == 0:
>         return 0
>     bits = n.bit_length()
>     a, b = divmod(bits, 2)
>     x = 2**(a+b)
>     while True:
>         y = (x + n//x)//2
>         if y >= x:
>             return x
>         x = y
>
>
> The second is only exact for some values of n, and for sufficiently large values
> of n, is will fail altogether:
>
>
> import math
>
> def isqrt_float(n):
>     """Integer square root using floating point sqrt."""
>     return int(math.sqrt(n))
>
>
>
> We know that:
>
> - for n <= 2**53, isqrt_float(n) is exact;
>
> - for n >= 2**1024, isqrt_float(n) will raise OverflowError;
>
> - between those two values, 2**53 < n < 2**1024, isqrt_float(n)
>   will sometimes be exact, and sometimes not exact;
>
> - there is some value, let's call it M, which is the smallest
>   integer where isqrt_float is not exact.
>
>
> Your mission, should you choose to accept it, is to find M.
>
> Hint: a linear search starting at 2**53 will find it -- eventually. But it might
> take a long time. Personally I gave up after five minutes, but for all I know
> if I had a faster computer I'd already have the answer.
>
> (But probably not.)
>
>
> Another hint: if you run this code:
>
>
> for i in range(53, 1024):
>     n = 2**i
>     if isqrt_newton(n) != isqrt_float(n):
>         print(n)
>         break
>
>
> you can find a much better upper bound for M:
>
>     2**53 < M <= 2**105
>
> which is considerable smaller that the earlier upper bound of 2**1024. But even
> so, that's still 40564819207303331840695247831040 values to be tested.
>
> On the assumption that we want a solution before the return of Halley's Comet,
> how would you go about finding M?

Here's a much smaller upper bound:

>>> n = 2 ** 53
>>> s = isqrt_newton(n)
>>> n
9007199254740992
>>> s
94906265
>>> m = (s+1)**2 - 1
>>> m
9007199326062755
>>> isqrt_newton(m) == isqrt_float(m)
False