5
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ A quick test on my computer with $MachinePrecision = 15.9546 show that 10^14+1.0 == 10^14. $\endgroup$ – John McGee Dec 25 '15 at 14:54
  • $\begingroup$ @JohnMcGee But doesn't isSq2 work with integers rather than reals? $\endgroup$ – rogerl Dec 25 '15 at 15:19
  • 1
    $\begingroup$ I believe that Compile of Sqrt yeilds only floating point result even for square integers, e.g sqrt2=Compile[{{n,_Integer},Sqrt@n]; $\endgroup$ – John McGee Dec 25 '15 at 15:51
5
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ But the problem is not with Sqrt, it's with ==: if you define isSq := IntegerQ@Sqrt@# &, that works properly with $2^{52}+1$ (that is, it returns False). Also, $2^{52}\approx 4\times 10^{15}$, which is much bigger than where I started to see failures of isSq2. So is there documentation anywhere of what tolerances == uses? (I don't really care what those tolerances are, but it would be nice to know in what range you could count on "exactness"). $\endgroup$ – rogerl Dec 25 '15 at 20:38
  • $\begingroup$ @rogerl Oops, I forgot I made a slight change to isSq2 -- I'll update. $\endgroup$ – Michael E2 Dec 25 '15 at 20:49
  • $\begingroup$ You are right about IntegerQ; I was just confused. But why should I start to see errors at $10^{14}$, which is over $45$ times smaller (over five bits) less than $2^{52}$? Again, my goal here is just to determine the range in which I can use this kind of integer function. $\endgroup$ – rogerl Dec 25 '15 at 20:49
  • $\begingroup$ @rogerl It follows that for the version of isSq2 in the OP, the borderline for the reliability of Floor[#] == # is around Floor@Sqrt[1/$MachineEpsilon/10^Internal`$EqualTolerance] or 3.51844*10^13. $\endgroup$ – Michael E2 Dec 25 '15 at 21:47
  • $\begingroup$ Thanks for the very detailed analysis. It clears things up for me quite a bit. One small comment: you say $MachineEpsilon == 2^-53; in fact it appears to be 2^-52. $\endgroup$ – rogerl Dec 26 '15 at 15:26

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.