Hi James,

The code you provided formatted correctly.

prim = With[{max = 100}, 
  Map[Last, 
   List @@ (Reduce[
       x^2 + y^2 == z^2 && 0 < x < y < max && 0 < z < max && 
        GCD[x, y, z] == 1, {x, y, z}, Integers, 
       Backsubstitution -> True] /. And :> List), {2}]]

It is hundreds of times slower than genPTunder that I referred to in my previous reply.

It looks like the triples were filtered in some way to generate the image you provided. So this does not generate an identical image. It does match the image on Wikipedia.

genPTunder[6000] // Map[Most] // Join[#, Map[Reverse, #]] & //
  ListPlot[#, AspectRatio -> 1, ImageSize -> 600, PlotStyle -> Blue] &

Missed that. Easy to modify

coprimePT = genPTunder[10000] // Map[Most] // Select[CoprimeQ @@ # &];
coprimePT // 
 ListPlot[{#, Map[Reverse, #]}, AspectRatio -> 1, ImageSize -> 600, 
   PlotStyle -> {Red, Black}] &

You may also find my method of generating PPTs useful. It generates all ppt's {a, b, c}, with the short leg as the input, where a<b<c. It is easily modified if you want all ppts with any leg equal to the input value by changing this part # > x to # > 0.

pyths[x_Integer] := (t = Select[Divisors[x^2], # <= x &]; 
  m = Sort[Select[(x^2 - t^2)/(2 t), IntegerQ[#] && # > x &]]; 
  Table[{x, m[[i]], Sqrt[x^2 + m[[i]]^2]}, {i, 1, Length[m]}])

or primitive ppt's

primpyths[x_Integer] := (t = Select[Divisors[x^2], # <= x &]; 
  m = Sort[Select[(x^2 - t^2)/(2 t), 
     IntegerQ[#] && # > x && GCD[x, #] == 1 &]]; 
  Table[{x, m[[i]], Sqrt[x^2 + m[[i]]^2]}, {i, 1, Length[m]}])

pyths[1000]

{{1000, 1050, 1450}, {1000, 1875, 2125}, {1000, 2400, 2600}, {1000, 
  3045, 3205}, {1000, 4950, 5050}, {1000, 6210, 6290}, {1000, 9975, 
  10025}, {1000, 12480, 12520}, {1000, 15609, 15641}, {1000, 24990, 
  25010}, {1000, 31242, 31258}, {1000, 49995, 50005}, {1000, 62496, 
  62504}, {1000, 124998, 125002}, {1000, 249999, 250001}}

primpyths[1000]

{{1000, 15609, 15641}, {1000, 249999, 250001}}

James, Thank you for the link, however I can't seem to get it to run, there seems to be some information missing so I can't see the end result. I did notice that you were using Reduce to get the triples which in my view is very slow, so I have taken the liberty of amending my method of producing PPTs to produce what I think was the desired output of your notebook.

pyths[x_Integer, max_Integer] := (t = Select[Divisors[x^2], # <= x &];
   m = Sort[
    Select[(x^2 - t^2)/(2 t), 
     IntegerQ[#] && # > x && Sqrt[x^2 + #^2] <= max && 
       GCD[x, #] == 1 &]]; 
  Table[{x, m[[i]], Sqrt[x^2 + m[[i]]^2]}, {i, 1, Length[m]}])

t1 = Flatten[Table[pyths[i, 1720], {i, 3, 1720}], 1];

t1 // Map[Most] // 
 ListPlot[{#, Map[Reverse, #]}, AspectRatio -> 1, ImageSize -> 600, 
   PlotStyle -> {Red, Black}] &

On my pc if they were all computed together runs in 270ms. Maybe this helps you.

Paul, I forgot to mention, I've already done the hard work on this notebook, just run this one, the long version, it only takes a few minutes at most: https://www.notebookarchive.org/ppt-graph-nb--2021-06-6y23baq/

JP

Thanks, JP