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

-
I guess it is good. Thanks for the Accept. – Mr.Wizard Dec 28 '12 at 10:43
Thanks for your help. – chyaong Dec 29 '12 at 13:49

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. – xzczd Dec 7 '12 at 10:37
@xzczd perhaps, but why fix what will remain broken (slow) when there is a better alternative? :^) – Mr.Wizard Dec 7 '12 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. – halirutan Dec 7 '12 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/… – Ajasja Dec 7 '12 at 11:42
+1 for learning to use Sow and Reap. – Mr.Wizard Dec 7 '12 at 14:53
@Mr.Wizard Needs["CompiledFunctionTools"]; CompilePrint@ptc ... Strange that the calls to MainEvalute don't seem to slow things down too much... – Ajasja Dec 7 '12 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 – Ajasja Dec 7 '12 at 15:39
Re: benchmarking, search for timeAvg – Mr.Wizard Dec 7 '12 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.

-