How to compile the code for generate Pythagorean_triple?

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 &) @@ # &]
]


5 Answers

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.

• I think OP is trying to figure out the mistake in his code…still +1 for the unbelievable fast code :D. Dec 7, 2012 at 10:37
• @xzczd perhaps, but why fix what will remain broken (slow) when there is a better alternative? :^) Dec 7, 2012 at 10:40

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 *)

• You have to check your compiled code. It still calls the kernel. Solve this by (1) using With in combination with CompilationOptions->{"InlineCompiledFunctions"->True} to really inline ccheck. (2) Don't use Tuples inside Compile. Use it outside and give the built tuples to the compiled function. This saves some tenth seconds. Dec 7, 2012 at 18:35

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}*)

• Have a look at InternalBag here mathematica.stackexchange.com/questions/845/… Dec 7, 2012 at 11:42
• +1 for learning to use Sow and Reap. Dec 7, 2012 at 14:53
• @Mr.Wizard Needs["CompiledFunctionTools"]; CompilePrint@ptc ... Strange that the calls to MainEvalute don't seem to slow things down too much... Dec 7, 2012 at 15:41

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 = InternalBag[Most[{0}]]},
Do[If[i^2 + j^2 == k^2,
InternalStuffBag[list, {i, j, k}, 2]], {i, m}, {j, i, m}, {k, j, m}];
InternalBagPart[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)

• Hmm, strange, If I run it like this or in it's own cell I get different timings. I think I'll have to improve the benchmarking for this case, but don't have time now Dec 7, 2012 at 15:39
• Re: benchmarking, search for timeAvg Dec 7, 2012 at 18:13

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.