# Faster square test for integers

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]

• I'd say the test based on the Jacobi symbol is still the best one of the lot, but of course it can't be embedded in Compile[]. Dec 16 '15 at 0:17
• What is the proper MSE protocol for situations like this, where I need an updated Answer to an existing but very old Question? I hesitated to knowingly post a Duplicate, but... Dec 16 '15 at 0:22
• There's always the possibility of using a bounty to ask for a needed update. Nevertheless, did you already test the ones using the Jacobi symbol? Otherwise, try this: NumberTheoryIntegerSqrt[n]^2 == n Dec 16 '15 at 0:28
• You don't need to load that context yourself; try executing NumberTheoryIntegerSqrt[10] in a fresh kernel if you don't believe me. Dec 16 '15 at 1:19
• @J.M. It's not that I didn't believe you, I just didn't know about Packages. LOL. I got it working. Thanks! Dec 16 '15 at 1:28

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.

• Wow, thanks for all this. As you probably gathered from my code, I'm on the MMa learning curve, and detailed answers like this really help with the learning. Your different approach to the math was very clever too. Dec 16 '15 at 21:18