Compilation, square roots, and integers

Question

After looking at this question, particularly this answer, I wrote my own performance test, using the two functions

isSq2 = Compile[{{n,_Integer}}, Floor@# == # & @ Sqrt @ n];

SquareQ08[n_] :=
     JacobiSymbol[n, 541] =!= -1 && JacobiSymbol[n, 547] =!= -1 &&
     JacobiSymbol[n, 557] =!= -1 && JacobiSymbol[n, 563] =!= -1 &&
     JacobiSymbol[n, 569] =!= -1 && JacobiSymbol[n, 647] =!= -1 &&
     JacobiSymbol[n, 653] =!= -1 && JacobiSymbol[n, 659] =!= -1 &&
     IntegerQ@Sqrt@n;

I found that indeed isSq2was much faster for "small" numbers. But I also found that isSq2 worked for larger numbers. Looking further, it seems that Compile works on machine-size integers, which on my box range up to $2^{63}-1$. So first question:

Is this true? Does Compile really work on all machine-size integers (which may vary depending on your box)? Is this explicitly documented anywhere other than that statement on the documentation page for Compile?

So assuming that the answer to the above question was true, I compared performace for larger integers, and ran across the following issue:

Table[isSq2[10^14 + i], {i, {-1, 0, 1}}]

(* {False, True, True} *)

Table[SquareQ08[10^14 + i], {i, {-1, 0, 1}}]

(* {False, True, False} *)

Clearly isSq2 erroneously returns True for $10^{14}+1$, which is far smaller than $2^{63}-1$. Further, testing random collections of 10000 numbers or so in that range (including the 10000 integers on either side of $10^{14}$), I can find no other integers for which the two functions return different values.

What is going on here? Is there something wrong with my code? And if so, why does it fail precisely for those two integers?

Edit: After the discussion below with @JohnMcGee, the issue is that compiled Sqrt always returns floats, and for the case I gave, 10^14+1 == 10^14. So the real reason that isSq2 does not work for larger integers is not that Compiled stuff fails for integers larger than $2^{31}$; it is because of roundoff error in floating point computations! So then my question is:

For what range of floating point values can I count on isSq2 doing the right thing? The documentation of == is a little vague in this regard.

Answer 1

I would guess Sqrt is sufficiently accurate up to about 2^53/2=2^52, where 53 is the number of bits representable in a floating-point number (binary64). This can be deduced from the series expansion of Sqrt near a perfect square:

Simplify[Series[Sqrt[x^2 + e], {e, 0, 2}], x > 0]
x + e/(2 x) - e^2/(8 x^3) + O[e]^3

When Abs[e/2] <= x^2 * $MachineEpsilon/2, the term e/(2 x) and higher terms will not change the sum, and the square root of x^2 + e will be (erroneously) given as x. When x and e are integers Floor[x] will equal x and isSq2 will return True. [See note at end.]

Since e is an integer, the smallest it can be is ±1. So up to x^2 < 2^52 then, since $MachineEpsilon == 2^-52, it should be possible to use the floating-point Sqrt to test for perfect squares. Since 2^52 is itself a perfect square, it should work for that value as well. But it should fail on 2^52 + 1. A simple test:

isSq2[2^52 - 1]
isSq2[2^52]
isSq2[2^52 + 1]
(*
  False
  True
  True
*)

The other issue in isSq2 is that == compares inexact numbers with a relative tolerance:

Approximate numbers with machine precision or higher are considered equal if they differ in at most their last seven binary digits (roughly their last two decimal digits). [From docs]

To overcome that, one should test x - y == 0. So the following fixes isSq2 for all integers up to 2^52:

isSq2 = Compile[{{n, _Integer}}, Floor@# - # == 0 &@Sqrt@n]

isSq2[10^14 + 1]
(*  False  *)

Further explanation of the tolerance. The tolerance is stored in Internal`$EqualTolerance. Here is a test of borderline cases:

SeedRandom[0];
Table[
 x = 10^RandomReal[{1, 100}];
 y = x + x*$MachineEpsilon/2*10^Internal`$EqualTolerance;
 x == y,
 {8}]
(*  {True, False, False, False, True, True, True, False}  *)

If we examine the OP's particular case of Sqrt[10^4 + 1], we can see that the difference between it (y) and its Floor (the same as x) is less than the tolerance.

x = Sqrt[10.^14];
y = Sqrt[10.^14 + 1];
x*$MachineEpsilon/2*10^Internal`$EqualTolerance
x - y
x == y
(*
  1.42109*10^-7
  -5.02914*10^-8
  True
*)

---

Note: How much to explain the borderline Abs[e/2] <= x^2 * $MachineEpsilon/2 depends on how well a reader understands the numerics of typical floating-point computations. Briefly, one often thinks roughly in terms of a continuous model; but floating-point computations are discrete, and the continuous model is only an approximation. Since we are analyzing edge cases, reasoning precisely and accurately can become somewhat intricate. Indeed my first argument, while correct in its conclusion, was specious (although perhaps not far from the truth). $MachineEpsilon is defined as "the difference between 1.0 and the next-nearest number representable as a machine-precision number." The factor $MachineEpsilon/2 in the borderline is divided by 2 because of rounding to the nearest floating-point number. Discretization becomes a factor at a power of 2, where the distance between succesively floating-point numbers jumps. Let us consider, then, when x^2, and therefore x, is a power of 2. The distance below a power of 2 to the next representable number is half the distance above, so the threshold at which Sqrt[x^2 + e] is x is different on each side of x. It is smaller below such an x than above. So it is more likely for the root to evaluate to x, when e > 0. Hence we use the greater $MachineEpsilon/2 as the factor in the borderline.