Challenge: find the first value where two functions differ

[Terry Reedy]

On 8/4/2017 11:44 AM, Chris Angelico wrote:
> On Sat, Aug 5, 2017 at 12:51 AM, Steve D'Aprano
> <steve+python at pearwood.info> wrote:
>> def isqrt_float(n):
>>      """Integer square root using floating point sqrt."""
>>      return int(math.sqrt(n))

The operations in the integer version are well-defined and so the answer 
should be the same on any Python 3 and indeed in any language that does 
proper integer arithmetic (but most do not).

In CPython, Math.sqrt wraps the C compiler's double sqrt function.  To 
make it deterministic, we must specify the float representation, 
rounding rule, and accuracy of the sqrt function.  Let's assume that the 
float answer is IEEE binary64, with 53 binary digits, with rounding rule 
'round to nearest, ties to even' and that the algorithm is as exact as 
possible, which is to say, the exact mathematical answer rounded as 
specified.

With 53 binary digits, all counts from 0 to 2**53 - 1 are exactly 
represented with a exponent of 0,  2**53 = 2**52 * 2, so it is exactly 
represented with an exponent of 1. Many other higher counts can be 
exactly represented with various exponents, but 2**53 + 1 requires 54 
digits.  So it is the first count that does not have an exact binary64 
representation.

Under these assumptions, how can isqrt_float(n) fail?

Suppose 0 <= m <= 2**26.5, so that m*m <= 2**53.  Under the assumptions 
that m.sqrt is as exact as possible, sqrt(m*m) == float(m) and 
int(float(m)) == m.  The true square root of m*m -1 will be m - e (e < 
1) and the integer square root is m-1.  However, if e is less than (and 
I believe equal) to 1/2 of the interval between n and the next lowest 
float, m - e will be rounded up to float(m) and ifloat(m*m-1) will 
incorrectly return m instead of m-1.

As m increases, e decreases while the interval between floats increases, 
so we should expect an eventual error

In this range of m, sqrt(m*m +1) will always be greater than m, so 
rounding down cannot lead to an error.  For larger m, where float(m*m) 
is inaccurate enough, it might.

As m increases, e decreases while the interval between floats increases, 
so we should expect an eventual error.  Whether it occurs for m <= 2* or 
not might be determined by careful analysis, but Chris Angelico already 
gave us the brute force answer.

>> We know that:
>>
>> - for n <= 2**53, isqrt_float(n) is exact;

but it is not ;-).  In this range, the only possibility is 
math.sqrt(m*m-1) being m instead of m - delta and that is the case for
Chris' answer 4503599761588224 = 67108865**2 - 1. 
int(math.float(4503599761588224) is 67108865 instead of 67108864.

>> - 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.


> Hmm. Thinking aloud a bit here. We know that isqrt_float(n) is not
> technically *exact* for any n that is not a square. So what I'd do is
> assume that for (n*n+1) to (n+1)*(n+1)-1, it's going to return the
> same value. I would then probe every perfect square, and the values
> one above and one below it. Now, that's still probably too many to
> check, but it's going to be a lot faster; for starters, we don't
> actually need to use isqrt_newton to compare against it.
> 
> for root in range(2**26, 2**53):
>      square = root * root
>      if isqrt_float(square - 1) != root - 1:
>          print(square - 1)
>          break
>      if isqrt_float(square) != root:
>          print(square)
>          break
>      if isqrt_float(square + 1) != root:
>          print(square + 1)
>          break
> 
> That gave me this result almost instantaneously:
> 
> 4503599761588224

Your limit was barely low enough ;-).

 >>> math.log2(math.sqrt(4503599761588224))
26.000000021497833

If you had started with, say int(2**26.5)-1, you would have missed this.
I reran with range(2, 2**26) and there is nothing lower.

 >>> math.log2(4503599761588224)
52.00000004299566

 >>> import decimal as d
 >>> I = d.Decimal(4503599761588224)
 >>> I.sqrt()
Decimal('67108864.99999999254941951410')
 >>> float(I.sqrt())
67108865.0

> which has been rounded up instead of down. I don't know if that counts
> as sufficiently wrong?

The isqrt_float answer is definitely wrong, and so is the claim for n < 
2**53.

>>>> isqrt_float(4503599761588224)
> 67108865
>>>> isqrt_newton(4503599761588224)
> 67108864
>>>> 67108865 ** 2
> 4503599761588225

-- 
Terry Jan Reedy