Challenge: find the first value where two functions differ

[MRAB]

On 2017-08-04 15:51, Steve D'Aprano 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?
> 
Why would isqrt_float not give the correct answer? Probably because of 
truncation (or rounding up?) of the floating point. I'd expect it to 
fail first near a square.