Thank you everybody for your help.

Related thread.

{2 m n, m^2 - n^2, m^2 + n^2}

is from this formula on wiki and you can easily verify that $(2mn)^2 + (m^2-n^2)^2 = (m^2 + n^2)^2$. The second part of the code is kind of clever way to generate

• Same pythagorean triples up to common factor
• Eleminate $m^2 + n^2$'s larger than the lim

Remove the last part of genPTunder gives:

    genPTunder2[lim_Integer?Positive] := 
     Module[{prim}, 
      prim = Join @@ 
        Table[If[
          CoprimeQ[m, n], {2 m n, m^2 - n^2, m^2 + n^2}, ## &[]], {m, 2, 
          Floor@Sqrt@lim}, {n, 1 + m~Mod~2, m, 2}]]

Notice that we need to add 2*{3,4,5}, 3*{3,4,5}, 4*{3,4,5} and remove {24,7,25} if we set lim to 20. And we also want to have each tuple sorted numertically:

In[9]:= data = genPTunder2[20]
Out[9]= {{4, 3, 5}, {12, 5, 13}, {8, 15, 17}, {24, 7, 25}}

Quotient function generate the common factors that can be used for multiplication of each triple:

In[11]:= {Range[20~Quotient~Max@#], Sort[#]} & /@ data
Out[11]= {{{1, 2, 3, 4}, {3, 4, 5}}, {{1}, {5, 12, 13}}, {{1}, {8, 15,
   17}}, {{}, {7, 24, 25}}}

$20/5 = 4$ therefore {3, 4, 5} can be multiplied up to four times so that the hypothenus is still no larger than 20. Similiarly $\lfloor 20/13 \rfloor = 1$, {5,12,13} is the only case; $\lfloor 20/25 \rfloor = 0$ so Range[20~Quotient~25] returns empty list.

KroneckerProduct works like map each multiplier to a given triple and the empty list is eaten by Union in the end:

In[13]:= KroneckerProduct[{1, 2, 3, 4}, {3, 4, 5}]
Out[13]= {{3, 4, 5}, {6, 8, 10}, {9, 12, 15}, {12, 16, 20}}
In[14]:= KroneckerProduct[{}, {7, 24, 25}]
Out[14]= {}

Because most of the funtion used here are either optimized for integer calculation, automatic multithreading or linear algebra operation, this SE implementation is very efficient for large number

Frank,

How about going from the other direction by generating only numbers that are squares of two other squared integers.

Original code on 2500 points:

In[1]:= AbsoluteTiming[
 t = Flatten[Table[{{x, y}, Sqrt[x^2 + y^2]}, {x, 2500}, {y, 2500}], 
    1];]

Out[1]= {24.1006, Null}

Generate all possible z's that fit the pattern:

In[2]:= AbsoluteTiming[zs = Map[#^2 &, Range[1, 2500]]; 
 zss = Total /@ Permutations[zs, {2}];]

Out[2]= {5.23768, Null}

Then find the sums of squares that are in the original list of squares:

vals = Intersection[zss, zs];

Now you can sort out which ones you want. Does this help?

Regards,

Neil