Faster square test for integers

Question

This question was asked already in Jan '12 and the most recent answer is from Oct '12, so it's several Mathematica versions out of date.

What is a faster test for whether an integer is a perfect square?

• For very large integers?
• For common-sized integers?

My code below plots the legs of Pythagorean triplets, but it's VERY slow because its calculating the Sqrt[].

isSq = Compile[{{n, _Integer}}, 
  With[{test = Sqrt[n]}, Floor[test] == test]];
t = {};
Do[
  bmax = Floor[a/Sqrt[2]];
  Do[
    If[GCD[a, b] == 1,
    If[
     isSq[a^2 + b^2],
     t = Append[t, {a, b}]]], {b, 1, bmax}], {a, 1, 20000}];
ListPlot[t]

Answer 1

It seems you are calculating legs of Pythagorean triples, $\{a,b\}$, for $b<a/\sqrt{2}$. I used your isSq function, and added a different test for GCD[a,b]==1.

RelativePrimesA[n_Integer] :=
   Block[{max = Floor[n/Sqrt[2]]},
      Complement[Range[max - 1], Apply[Sequence, 
         Map[Range[#, max - 1, #] &, FactorInteger[n][[All, 1]]]]]
   ]

The following code is about 6 times faster than your original code, for an upper limit on a of 10000. Using a ParallelTable with 8 kernels, results in a solution about 30 times faster.

ListPlot[
   Flatten[DeleteCases[
      Table[
         Thread[{a, Pick[#, Map[isSq, a^2 + #^2]] &[RelativePrimesA[a]]}],
         {a, 2, 10000}],
      {}], 1]]

Another method, which approaches the problem slightly differently, is as follows. Find a primitive sum of two squares equalling a square n, then impose the criterion that $b<a/\sqrt{2}$.

PrimitiveSumTwoSquares[n_] :=
   Block[{r, a, s},
      If[(r = PowerModList[-1, 1/2, n]) == {}, {},
      Table[
         a = n;
         s = r[[i]];
         While[a^2 >= n, {s, a} = {a, Mod[s, a]}]; (* do the GCD *)
         {Mod[s, a], a}, (* one last iteration gives solution {a,b}, *)
         {i, Length[r]/2}] (* conjecture half the length *)
      ]]

Now test a range of squares n.

Block[{p},
   ListPlot[
      Map[Reverse, Flatten[DeleteCases[
         Table[
            Pick[p=PrimitiveSumTwoSquares[n], Map[#[[2]]/Sqrt[2] > #[[1]]&, p]],
            {n, Range[2, 100000]^2}],
         {}], 1]]]]

This plot took about 4 seconds.