6
$\begingroup$

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]

PythagoreanTripletLegs

$\endgroup$
  • $\begingroup$ 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[]. $\endgroup$ – J. M. is away Dec 16 '15 at 0:17
  • $\begingroup$ 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... $\endgroup$ – Jerry Guern Dec 16 '15 at 0:22
  • $\begingroup$ 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 $\endgroup$ – J. M. is away Dec 16 '15 at 0:28
  • $\begingroup$ You don't need to load that context yourself; try executing NumberTheory`IntegerSqrt[10] in a fresh kernel if you don't believe me. $\endgroup$ – J. M. is away Dec 16 '15 at 1:19
  • $\begingroup$ @J.M. It's not that I didn't believe you, I just didn't know about Packages. LOL. I got it working. Thanks! $\endgroup$ – Jerry Guern Dec 16 '15 at 1:28
7
$\begingroup$

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.

Pythagorean Legs

$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – Jerry Guern Dec 16 '15 at 21:18

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.