9
$\begingroup$

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 &) @@ # &]
 ]
Improve this question
$\endgroup$
0

5 Answers 5

Reset to default
22
$\begingroup$

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.

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

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 *)
Improve this answer
$\endgroup$
1
  • $\begingroup$ 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. $\endgroup$
    – halirutan ♦
    Dec 7, 2012 at 18:35
6
$\begingroup$

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}*)
Improve this answer
$\endgroup$
3
  • 3
    $\begingroup$ Have a look at Internal`Bag here mathematica.stackexchange.com/questions/845/… $\endgroup$
    – Ajasja
    Dec 7, 2012 at 11:42
  • $\begingroup$ +1 for learning to use Sow and Reap. $\endgroup$
    – Mr.Wizard
    Dec 7, 2012 at 14:53
  • $\begingroup$ @Mr.Wizard Needs["CompiledFunctionTools`"]; CompilePrint@ptc ... Strange that the calls to MainEvalute don't seem to slow things down too much... $\endgroup$
    – Ajasja
    Dec 7, 2012 at 15:41
5
$\begingroup$

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)

Improve this answer
$\endgroup$
2
  • $\begingroup$ 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 $\endgroup$
    – Ajasja
    Dec 7, 2012 at 15:39
  • $\begingroup$ Re: benchmarking, search for timeAvg $\endgroup$
    – Mr.Wizard
    Dec 7, 2012 at 18:13
1
$\begingroup$

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.

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.