Compile関数の威力

code: https://colab.research.google.com/drive/1woAe7wRlAS4LeebCssP45DMv9jq8X8cM

オリジナルコード

壁時間で約10秒かかりました。

In [ ]:
newtonOriginal[z_Complex, n_Integer] :=
  Module[{zn = z},
     Do[zn = (2 zn + 1/zn^2)/3, {n}]; 
     If[Re[zn] > 0, 1, If[Im[zn] > 0, 2, 3]]]
In [ ]:
AbsoluteTiming[
   no = Table[
       newtonOriginal[x + I y, 25], 
       {y, -1, 1, 2./399}, {x, -1, 1, 2./399}];]
Out[ ]:
Output
In [ ]:
ArrayPlot[no]
Out[ ]:
Output

コンパイル

同じコードをコンパイルすることで約20倍速くなります。

In [ ]:
newt = Compile[{{z, _Complex}, {n, _Integer}}, 
   Module[{zn = z},
      Do[zn = (2 zn + 1/zn^2)/3, {n}]; 
      If[Re[zn] > 0, 1, If[Im[zn] > 0, 2, 3]]]];
    
AbsoluteTiming[
     nc = Table[newt[x + I y, 25], {y, -1, 1, 2./399}, {x, -1, 1, 2./399}];]
Out[ ]:
Output

Listable

期待したほどの性能向上は残念ながら見られませんでした。

In [ ]:
newt = Compile[{{z, _Complex}, {n, _Integer}}, 
   Module[{zn = z},
      Do[zn = (2 zn + 1/zn^2)/3, {n}]; 
      If[Re[zn] > 0, 1, If[Im[zn] > 0, 2, 3]]], 
   RuntimeAttributes -> Listable];
  
AbsoluteTiming[
   nl = Table[newt[x + I y, 25], {y, -1, 1, 2./399}, {x, -1, 1, 2./399}];]
Out[ ]:
Output

並列化

コンパイル対象コードに If があるせいか、性能向上は若干に留まっています。

In [ ]:
newtp = Compile[{{z, _Complex}, {n, _Integer}}, 
   Module[{zn = z},
      Do[zn = (2 zn + 1/zn^2)/3, {n}]; 
      If[Re[zn] > 0, 1, If[Im[zn] > 0, 2, 3]]], 
   RuntimeAttributes -> Listable, Parallelization -> True];
   
AbsoluteTiming[
 np = newtp[Table[x + I y, {y, -1, 1, 2./399}, {x, -1, 1, 2./399}], 
                       25];]
Out[ ]:
Output
In [ ]:
ArrayPlot[np]
Out[ ]:
Output