How to compile the code for generate Pythagorean_triple?

Question

I am finding Pythagorean_triple, it worked slowly. I tried to compile, but it gives some warnings. I also use "Case" or "Do" ,both of them failed.I'm sure my CCompiler has been set correctly. How can I compile the following code?

With[{m = 200},
 Select[Flatten[Table[{x, y, z}, {x, m}, {y, x, m}, {z, y, m}], 
   2], (#1^2 + #2^2 == #3^2 &) @@ # &]
 ]

Answer 1

There are much faster ways to generate Pythagorean triples.

Update: Now twice as fast.

genPTunder[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}
   ];
  Union @@ (Range[lim ~Quotient~ Max@#] ~KroneckerProduct~ {Sort@#} & /@ prim)
 ]

genPTunder[50]

{{3, 4, 5}, {5, 12, 13}, {6, 8, 10}, {7, 24, 25}, {8, 15, 17},
 {9, 12, 15}, {9, 40, 41}, {10, 24, 26}, {12, 16, 20}, {12, 35, 37},
 {14, 48, 50}, {15, 20, 25}, {15, 36, 39}, {16, 30, 34}, {18, 24, 30},
 {20, 21, 29}, {21, 28, 35}, {24, 32, 40}, {27, 36, 45}, {30, 40, 50}}

genPTunder[100000] // Length // Timing

{0.125, 161436}

Over 160,000 triples in an eighth of a second should be serviceable, even without compilation.

Answer 2

Here is one idea; with more time one could think of a better way to generate the table.

ccheck = Compile[{{list, _Integer, 1}},   list[[1]]^2 + list[[2]]^2 == list[[3]]^2, CompilationTarget :> "C"]

pt = Compile[{{m, _Integer}}, 
 Select[Select[Tuples[{Range[m], Range[m], Range[m]}], #[[1]] <= #[[2]] <= #[[3]] &], 
        ccheck[#] &], {{ccheck[_], True | False}}, 
  CompilationTarget :> "C"]

output1 = With[{m = 200}, Select[Flatten[Table[{x, y, z}, {x, m}, {y, x, m}, {z, y, m}], 2], (#1^2 + #2^2 == #3^2 &) @@ # &]] // AbsoluteTiming;

output2 = pt[200] // AbsoluteTiming;

output1[[1]]
(* 8.918892 *)

output2[[1]]
(* 1.250000 *)

output1[[2]]==output2[[2]]
(* True *)

Answer 3

I already knew a fast way, but I didn't know how to let it gives a list.

cpt = Compile[{{m, _Integer}},
   Do[If[i^2 + j^2 == k^2, Print[{i, j, k}]], {i, m}, {j, i, m}, {k, j, m}],
      CompilationTarget -> "C"
   ];

cpt[1000] // Timing

==============================

update, a neat version, but it's not my original:

  ptc = Compile[{{m, _Integer}},
   Do[If[i^2 + j^2 == k^2, Sow@{i, j, k}], {i, m}, {j, i, m}, {k, j, m}], 
     CompilationTarget -> "C", RuntimeOptions -> "Speed"
   ];

 ptc[1000] // Reap // Last // First // Length // Timing
(*{0.359, 881}*)

Answer 4

Here is a compiled version of the for loop:

PTG[m_] := 
  Select[Flatten[Table[{x, y, z}, {x, m}, {y, x, m}, {z, y, m}], 
    2], #[[1]]^2 + #[[2]]^2 == #[[3]]^2 &];
PTGC = Compile[{{m, _Integer}},
   Block[{list = Internal`Bag[Most[{0}]]}, 
    Do[If[i^2 + j^2 == k^2, 
      Internal`StuffBag[list, {i, j, k}, 2]], {i, m}, {j, i, m}, {k, j, m}];
    Internal`BagPart[list, All]
    ], CompilationTarget -> C, RuntimeOptions -> "Speed"];

Here are some timings (the last one is @b.gatessucks solution):

PTG[200] // Length // Timing
Partition[PTGC[200], 3] // Length // Timing
pt[200] // Length // Timing

(*
{5.741, 127}
{0., 127}
{0.78, 127}
*)

(But of course I'd go with mr.W solution) (I was to lazy to add RuntimeOptions -> "Speed", but in this case it really helps)

Answer 5

This generates all primitive pythagorean triples:

pythT[triple_] := 
triple.# & /@ {{{1, 2, 2}, {-2, -1, -2}, {2, 2, 3}}, {{1, 2, 2}, 
{2, 1, 2}, {2, 2, 3}}, {{-1, -2, -2}, {2, 1, 2}, {2, 2, 3}}}
pythT2[triples_] := Join[Flatten[pythT@# & /@ triples, 1], triples]
pythN[n_] := Join[{{3, 4, 5}}, DeleteDuplicates[Sort@Nest[pythT2, pythT@
{3, 4, 5}, n]]]

Last@pythN@10

(*{927538921, 927538920, 1311738121}*)

but it grows like pythNlength[n_] := (3^(n + 2) - 1)/2, so only very low n is needed.

Answer 6

A variation of Mr.Wizard, a little faster:

Clear["`*"];
getPrimitiveTriple = Compile[{{n, _Integer}},
   Module[{bag = Internal`Bag[Most@{0}]},
    Do[
     If[GCD[i, j] == 1,
      Internal`StuffBag[bag, Sort@{j^2 - i^2, 2 i j, i^2 + j^2}, 1]],
     {i, (Sqrt[2 n - 1] - 1)/2}, {j, i + 1, Sqrt[n - i^2], 2}];
    #[[Ordering[#[[All, 3]]]]] &@Partition[Internal`BagPart[bag, All], 3]
    ]
   ];

getAllTriple[n_Integer] := SortBy[Join @@ (Outer[Times, Range[n/Max@#], #] & /@ 
  getPrimitiveTriple[n]), Last];

getAllTriple[50]

{{3, 4, 5}, {6, 8, 10}, {5, 12, 13}, {9, 12, 15}, {8, 15, 17}, {12, 16, 20}, {7, 24, 25}, {15, 20, 25}, {10, 24, 26}, {20, 21, 29}, {18, 24, 30}, {16, 30, 34}, {21, 28, 35}, {12, 35, 37}, {15, 36, 39}, {24, 32, 40}, {9, 40, 41}, {27, 36, 45}, {14, 48, 50}, {30, 40, 50}}

getAllTriple[10^6] // Length // AbsoluteTiming

{1.08712,1980642}