2013/03/19

When the inexact square root of an integer is exact

In Exact sqrt entry, I wrote that, for a natural number smaller than 2^53, we could use double-precision floating point sqrt and take the answer exact if the result was an integer. As Mark H Weaver pointed out, it was wrong.

Suppose we have nonnegative integers N, M, m, and non-zero real number e, where M = m^2 < 2^53, N = (m*(1+e))^2 < 2^53.

If the absolute value of e isn't greater than 2^-53, square root of N, which is m*(1+e), can be rounded to m in the double-precision sqrt calcluation. The result becomes an integer, but not exact. That's the case we want to exclude.

N - M = (2e+e^2)M. The maximum bound of (2e+e^2) is (2^-52 + 2^-106) when e = 2^-53, so if M is greater than 2^52, N - M can exceed 1 and we might incorrectly recognize N as a square number of m.

The greatest square number below 2^52 is (2^26-1)^2 = 4503599493152769, and (2^-52 + 2^-106)*4503599493152769 is 81129635996755064984528658366465/81129638414606681695789005144064, which is smaller than 1.

2^52 = 4503599627370496 is also a square number, and calculating (sqrt (inexact (+ (expt 2 52) 1))) yields a rounded integer 67108864.0. But the lower bound, (-(2^52) + 2^-106)*4503599627370496 is -18014398509481983/18014398509481984, which is grater than -1, so applying inexact sqrt on 2^52-1 won't lead us to the rounded integer. (Indeed, (sqrt (inexact (- (expt 2 52) 1))) yields 67108863.99999999.)

So we can use inexact square root when the input is smaller than 2^52.

Tags: Flonums, sqrt, exact-integer-sqrt