# We keep only primitive pythagorean triples
pyth := n ->
  Filtered(Cartesian([1 .. n], [1 .. n], [1 .. n]),
  u -> u[3]^2 = u[1]^2 + u[2]^2 and u[1] < u[2]
  and GcdInt(u[1], u[2]) = 1);

pyth(100);
# [ [ 3, 4, 5 ], [ 5, 12, 13 ], [ 7, 24, 25 ], [ 8, 15, 17 ], [ 9, 40, 41 ], [ 11, 60, 61 ], [ 12, 35, 37 ],
#   [ 13, 84, 85 ], [ 16, 63, 65 ], [ 20, 21, 29 ], [ 28, 45, 53 ], [ 33, 56, 65 ], [ 36, 77, 85 ], [ 39, 80, 89 ],
#   [ 48, 55, 73 ], [ 65, 72, 97 ] ]
