Challenge: find the first value where two functions differ

[Steve D'Aprano]

On Sat, 5 Aug 2017 06:17 am, Terry Reedy wrote:

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

Indeed. But although 2**53 + 1 cannot be represented as a float, there is no
difference between its integer square root and the integer square root of
2**53. Both are 94906265.


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

Indeed. I had forgotten one possibility.


I was thinking that there are only two cases when n is less than 2**53:

(1) n is a perfect square, in which case sqrt(n) is also a whole number which
can be represented exactly, and sqrt(n) is the same as integer sqrt(n);

(2) or n is not a perfect square, in which case sqrt(n) returns a number of the
form a.b, where a is the integer sqrt of n.

But in fact there's a third possibility: if the .b part of a.b is sufficiently
large that it rounds up to (a+1).0. It was this case that I had forgotten
about: the possibility of rounding up to the next whole number.

That's exactly what happens here with Chris' case. The exact mathematically
correct value of sqrt(4503599761588224) is a surd, but we can use Decimal to
get an approximate answer:


py> from decimal import Decimal
py> m = Decimal(4503599761588224)
py> m.sqrt()
Decimal('67108864.99999999254941951410')

Now 67108864.99999999254941951410... is not exactly representable as a binary
float. The two possibilities are:

67108864.99999999  # round down
67108865.00000000  # round up

Rounding down is *closer* to the correct value, which suggests that rounding up
is a bug:

py> m.sqrt() - Decimal('67108864.99999999')  # round down
Decimal('2.54941951410E-9')
py> m.sqrt() - Decimal('67108865.0')  # round up
Decimal('-7.45058048590E-9')


Am I right that math.sqrt(4503599761588224) should round down?


-- 
Steve
“Cheer up,” they said, “things could be worse.” So I cheered up, and sure
enough, things got worse.