Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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 &) @@ # &]
 ]
share|improve this question
    
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
add comment

4 Answers

up vote 10 down vote accepted

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.

share|improve this answer
    
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
3  
@xzczd perhaps, but why fix what will remain broken (slow) when there is a better alternative? :^) –  Mr.Wizard Dec 7 '12 at 10:40
add comment

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 *)
share|improve this answer
    
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
add comment

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}*)
share|improve this answer
1  
Have a look at Internal`Bag 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
add comment

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)

share|improve this answer
    
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
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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