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 = {};
  bmax = Floor[a/Sqrt[2]];
    If[GCD[a, b] == 1,
     isSq[a^2 + b^2],
     t = Append[t, {a, b}]]], {b, 1, bmax}], {a, 1, 20000}];


Jerry Guern

Posted 2015-12-16T00:03:10.713

Reputation: 4 197

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[]. – J. M.'s ennui – 2015-12-16T00:17:14.537

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... – Jerry Guern – 2015-12-16T00:22:11.477

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: NumberTheory`IntegerSqrt[n]^2 == n – J. M.'s ennui – 2015-12-16T00:28:24.397

You don't need to load that context yourself; try executing NumberTheory`IntegerSqrt[10] in a fresh kernel if you don't believe me. – J. M.'s ennui – 2015-12-16T01:19:00.330

@J.M. It's not that I didn't believe you, I just didn't know about Packages. LOL. I got it working. Thanks! – Jerry Guern – 2015-12-16T01:28:58.660



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.

         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]) == {}, {},
         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.

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

This plot took about 4 seconds.

Pythagorean Legs


Posted 2015-12-16T00:03:10.713

Reputation: 14 269

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. – Jerry Guern – 2015-12-16T21:18:25.670