円周率の近似値を手元に置くこと

B00264YQ80

円周率の近似値が印刷されたマグカップを手元に置くのが一つの手です。

487310002X

円周率の近似値が印刷された書籍、牧野貴樹『円周率1000000桁表』を手元に置くのも一つの手です。

必要な時にすぐ計算できるように、スーパーπをダウンロードしておくのも一つの手です。スーパーπはかなり前からある有名なプログラムで、ベンチマークなどでもよく使われています。本稿執筆時点で、419万桁の計算時間を集めたサイトスーパーπランキングのトップは、Intel i7 2600K@5.454GHz (Sandy Bridge)の35秒とのこと。私のマシン(Core i7 930)では121秒でした。

Mathematicaを使うのも一つの手です。スーパーπのベンチマークと同じだけ計算したい場合は、「pi = ToString@N[Pi, 4194307];」とします。有効な結果はStringDrop[pi, -2]です(最後の2桁を捨てる)。私のマシン(Core i7 930)では16秒と、スーパーπよりかなり速く計算できました(もちろん計算結果は同じです)。

結果をそのままファイルに書き出しても4MB程度なので、「記録しておくより、必要な時に計算したほういい」という状況ではありませんが。

ライフゲーム上の万能チューリングマシン

『ライフゲイムの宇宙』という本を以前紹介しました。そこで、

ところで、ライフゲイムで万能チューリング・マシンを作るために必要なセルは10^13程度と見積もられている。これはクラスタ計算機なら到達できるところだろう。視覚化の問題はあるが、誰かやらないだろうか。

などと言っていたのですが、クラスタを使うまでもなく、2002年に実現しているということを教えてもらいました[Paul Chapman]。

万能チューリングマシンとは何か。チューリングマシンと万能チューリングマシンの違いは、電卓とパソコンの違いのようなものだと考えてください。電卓は、あらかじめ組み込まれたタスクしかこなせませんが、パソコンは、プログラムが与えられればどんなタスクでもこなせます。万能チューリングマシンは、適切な入力が与えられればどんなチューリングマシンにもなれます。チューリングマシンでできることが「計算可能なこと」だと考えられているので、万能チューリングマシンは、計算可能なものならどんなものでも計算できるのです。計算がすぐに終わるかどうかは別問題です。

ライフゲームが興味深い理由の一つに、ライフゲーム上に万能チューリングマシンを構築できることがあるのは間違いないでしょう。しかし、そのための空間はとても大きく、『ライフゲイムの宇宙』当時には難しかったわけです。

しかし、計算技術の進歩によって、ふつうのパソコンで実現できるようになったというのは先に述べたとおりです。2011年になっても改良が加えられ、Golly用のファイルも配布されたりしているので、簡単に動かして試せるようになっています[Paul Rendell]。

万能チューリングマシンをシミュレートするためのライフゲームの世界はとても大きいです。1:16くらいで表示するとこんな感じですが、

Golly Scale=1:16

全体を表示できるようにしたときの縮尺はなんと2^8:1です。ステップが進むにつれて大きくなります。

Golly Scale=1:16

こういうのは、先日紹介したMathematicaによる手軽な方法ではまったく歯が立たないでしょうね。

4150112908『ライフゲイムの宇宙』は原書が1985年、翻訳が1990年です。発行から時間が経っているので、このように、内容が古くなっている部分があります。他にもたとえば、p.43に載っているエデンの園(Garden of Eden)は9×33のパターンですが、2009年には11×11のものが発見されています(参考:LifeWiki)。

「エデンの園配置とはなにか、知っていますか?」

「ええ、もちろん。セル・オートマトン理論において、あるシステムの状態が、それ以前のいかなる状態の結果でもありえないことをいうの。ほかのいかなるセルのパターンも、その状態を生じさせられない。エデンの園配置を望むなら、それそのものからはじめなくてはならない—それをシステムの最初の状態として、自ら設定する必要があるわけ」(イーガン『順列都市』(下p.20))

バットマン方程式

バットマン方程式というのがあるそうです(Is this Batman equation for real?)。なんでも、この式を満たす(x, y)をプロットするとバットマンのマークになるのだとか。

batman equation

さっそくMathematicaのContourPlotで描こうとしたら、うまくいきませんでした。

仕方がないので、こんな感じにx方向に走査しながら方程式を満たす点を探して、一応描けました。

f[x_, y_] := ((x/7)^2 Sqrt[
      Abs[Abs[x] - 3]/(Abs[x] - 3)] + (y/3)^2 Sqrt[
      Abs[y + (3 Sqrt[33])/7]/(y + (3 Sqrt[33])/7)] - 
    1) (Abs[x/2] - ((3 Sqrt[33] - 7)/112) x^2 - 3 + 
    Sqrt[1 - (Abs[Abs[x] - 2] - 1)^2] - 
    y) (9 Sqrt[
      Abs[(Abs[x] - 1) (Abs[x] - 3/4)]/((1 - Abs[x]) (Abs[x] - 
           3/4))] - 8 Abs[x] - 
    y) (3 Abs[x] + .75 Sqrt[
      Abs[(Abs[x] - 3/4) (Abs[x] - 1/2)]/((3/4 - Abs[x]) (Abs[x] - 
           1/2))] - 
    y) (9/4 Sqrt[Abs[(x - 1/2) (x + 1/2)]/((1/2 - x) (1/2 + x))] - 
    y) ((6 Sqrt[10])/
     7 + (3/2 - Abs[x]/2) Sqrt[
      Abs[Abs[x] - 1]/(Abs[x] - 1)] - (6 Sqrt[10])/14 Sqrt[
      4 - (Abs[x] - 1)^2] - y)

points =
  Flatten[Table[{x, #} & /@ Cases[y /. NSolve[f[x, y] == 0, y], _Real],
    {x, -8, 8, 0.1}], 1];

ListPlot[points]

batman equation1

Mathematicaを持っている人は、pointsの計算の前にDistributeDefinition[f]を実行し、Tableの代わりにをParallelTableを使った方が速いでしょう。

こんなもんかなと思っていたら、Mathematicaできれいに描く方法がPlaying With Mathematicaで紹介されていました(The Batman curve)。パーツに分けてしまえば描けて当然という気もしますが、こっちのほうがきれいですね。

batman equation1

追記:Wolfram Alphaが対応したので、Mathematicaでも簡単に描けるようになりました。

Universal Mathematica Manipulator—Poor Man’s Mathematica

As of December 2017, the method introduced here seems to be difficult to realize.

Newer Version(新しいバージョン)

Universal Mathematica Manipulator

Wolfram CDF Player is a free of charge software to view documents in Computable Document Format (CDF). Users cannot edit codes using the software. For editing, Mathematica is required. Consequently, the relationship between CDF player and Mathematica is similar to the one between Adobe Reader and Acrobat.

Since CDF Player has almost all Mathematica functions to allow rich user experiences, users should take full advantage of those functions.

Although users cannot edit CDF documents with CDF player, they can still input data. If there were a document with a universal manipulator that receives program codes and outputs the evaluated results of them, then the document can perform as poor man’s Mathematica, or a complimentary alternative for Mathematica.

Users cannot use the ordinary Mathematica expressions as inputs, because, according to Wolfram research, CDF does not support non-numeric input fields [FAQ]. However, any expressions can be represented with numbers (e.g., Gödel number), and CDF supports numbers. Consequently, if there were a converter that converted Mathematica expressions to numbers and vice versa, CDF documents could process Mathematica expressions by using the converter.

Then, I developed such a converter. As a result, a CDF document that can process any Mathematica documents with the converter became reality. Therefore, I call the CDF document “Universal Mathematica Manipulator (UMM).”

Yet, one little problem remains. Users cannot paste on CDF documents unless he/she has Mathematica. As numbers that represent Mathematica expressions are very large, it takes a great deal of trouble to input them digit by digit. Luckily, both Mac OS X and Windows have scripting languages that can simulate key strokes: AppleScript for Mac OS X and VBScript for Windows.

The script in AppleScript is given below. It obtains a number from clipboard and emulates key strokes to input it digit by digit.

tell application "System Events"
     set str to get the clipboard
     keystroke "a" using command down
     keystroke str
end tell

This script is created on AppleScript Editor and saved in User Script Folder (~/Library/Scripts). Make sure that Script menu is shown in the menu bar in AppleScript Editor’s preferences. One downside to the script in AppleScript is to have to hit either “Enter” or “Tab” key to resume the script when it stops. Though, the way to eliminate the extra step is unknown.

Also, the script in VBScript is below. It is better than the one in AppleScript, as we do not need to press the enter key after the number is inputted. Nevertheless, this script has a limitation where it cannot handle an input with 5,000 or longer in length.

Set objIE = CreateObject("InternetExplorer.Application")
objIE.Navigate("about:blank")
text = objIE.document.parentwindow.clipboardData.GetData("text")
objIE.Quit

Set objShell = WScript.CreateObject("WScript.Shell")
WScript.Sleep(3000)
objShell.SendKeys("^a0{Enter}")
WScript.Sleep(1000)
objShell.SendKeys(text & "{Enter}")

In executing the script in VBScript, a warning message of Internet Explorer appears. You can disable the message at your own risk. (Tools > Internet Options > Security tab > Custom level… > Scripting > Allow Programmatic clipboard access > Enable)

This system is tested on Wolfram CDF Player 8.0.3.0 Mac OS X x86 and Microsoft Windows (32-bit). Relatively large codes like a solution for Burrau’s problem of three bodies can be executed.

solution for Burrau's problem of three bodies on UMM

Newer Version(新しいバージョン)

CUDAを使ってマンデルブロ集合の描画を3桁速くする方法

追記:Mathematica 10で導入されたMandelbrotSetPlot[]は、GPUを使いませんがかなり速いです。この記事でやっていることは、以下のコードで再現できます。

追記:OpenCLの場合

描画が3桁速くなると、マンデルブロ集合をインタラクティブに楽しめるようになります。

はじめに

以前、Mathematicaでマンデルブロ集合を描く方法を紹介しました(「追悼、マンデルブロ博士」)。

Clear@f;
f[c_] := f[c, 0., 0]
f[_, _, 100] = -1;
f[c_, z_, i_] := With[{x = z^2 + c}, If[2 < Abs@x, i, f[c, x, i + 1]]]
DensityPlot[f[x + y I], {x, -2, 1}, {y, -1.5, 1.5}, PlotPoints -> 100]

描くのは確かに簡単だったのですが、「ここをもうちょっと拡大したい」という要望に応えるためのインタラクティブな仕掛けは、遅すぎて使えませんでした。「追悼記事なのに・・・」という思いがあったので、ちょっと改善することにしました。

目的

計算を高速化して、インタラクティブなマンデルブロ集合描画システムを作ります。

手法

(1)計算自体(上の例のf)高速化と(2)並列化で高速化します。(1)のために、再帰を反復に書き換えます。(2)のために、CUDAを使います。マンデルブロ集合の計算は、いわゆる「embarrassingly parallel」というやつなので、CUDAの練習に向いていると思います。

この記事のようにインタラクティブなものではありませんが、Mathematicaのマニュアルにもマンデルブロ集合の例が掲載されています。

条件

時間を計る条件は次のようにします。「有理数だと遅い」と思われるかもしれませんが、浮動小数点数を使うと、TableとParallelTableの結果が違ってしまう危険があることへの対応です。ちなみに、これはバグだと思うのですが、Wolfram社の回答は「バグでは無い」でした。

xMin = -2;
xMax = 1;
yMin = -3/2;
yMax = 3/2;
dx = (xMax - xMin)/300;
dy = (yMax - yMin)/300;

計算機環境(2013年3月更新)

  • CPU:Core i7 950(HTオフ)
  • 主記憶:12GB
  • GPU:NVIDIA GeForce GTX 580
  • OS:Windows 7 64 bit
  • Mathematica 9.0.1

結果

比較の基準

最初のバージョンをCUDAを使うバージョンで、計算時間を比較してみましょう。

最初のバージョンの実行時間を次のようなコードで計ると11.6秒でした。

CUDA

CUDAを使うバージョンの実行時間を次のようなコードで計ると0.003秒でした。4000倍くらい速くなることがわかります。

まあ、GPU無しの比較対象は最初に思いつくものですから、「まずGPU無しでどこまで速くなるのか追求してから」と思う方は、下の「おまけ」を読んでください。GPU無しでも最初のものより80倍くらいは速くできるので、それと比べれば、「GPUの効果は50倍くらい」というのが妥当なところでしょう。

うまくいかないときは、CUDAのための環境がちゃんとあることを確認してください。Visual StudioとCUDAがインストールされていればいいはずです(OpenCL版の方が簡単かもしれません)。私は、Visual Studio 2008 Professionalのデフォルトではx64のコンパイラがインストールされないために少しはまりました。本稿執筆時点でMathematicaが使うデフォルトのCUDA SDKは3.1で、これがサポートするのはVisual Studio 9 (2008)までだということにも注意が必要です(カスタマイズすればいいのでしょうが)。追記:2011年8月にCUDA SDK 4に対応したので、Visual Studio 2010がそのまま使えるようになりました。

リアルタイム描画

描画には、ArrayPlotを使うことにします。ListDensityPlotのほうが高度ですが、少し遅いのと、座標の扱いに問題があったため使えませんでした。Imageは速そうですが、やはり座標の扱いに問題がありました。

ユーザインターフェースはMathematicaのManipulateで構築します。

realtime

画像上をクリックすると、その場所を中心にして再描画します。Scaleで描画領域の広さを、Log[maxSteps]で反復の上限を変更できます。

ぬるぬる動くというわけにはいきませんが、「ここをもうちょっと拡大したい」という要望には十分応えられると思います。もっと速くしたい場合は、ArrayPlotをどうにかするということになるのでしょうが、もう少しコードを書かなければならないでしょう。

ところで、さすがにこれってCDFにしてもだめなんですよねえ。

おまけ(細かいことなので、これ以降は読まなくていいです)

追悼記事がダメだった理由

追悼記事がダメだった最大の要因は、DensityPlotにあったようです。DensityPlotは、下に示すように場所ごとにメッシュの細かさを自動調整してくれる便利な関数なのですが、この内部が並列化されていないため、あまり速くならないのです。

densityplot

CUDAなしでどこまで行くか

コンパイル+並列化

Mathematicaでは、機械精度の数を使うことに限定すれば、ユーザ定義関数をコンパイルできます。最初のバージョンをコンパイルし、並列化して実行すると30倍くらい速くなります。

fc = Compile[{{c, _Complex}, {z, _Complex}, {i, _Integer}},
   If[i == 100, -1, 
    With[{x = z^2 + c}, If[2 < Abs@x, i, fc[c, x, i + 1]]]], 
   CompilationTarget -> "C"];
DistributeDefinitions[fc];

time = AbsoluteTiming[
   result = ParallelTable[ fc[x + y I, 0, 0],
      {y, yMin, yMax, dy}, {x, xMin, xMax, dx}];][[1]];

time0/time

28.44619
反復

再帰を反復に書き直すともっと速く、80倍くらい速くなります。

hc = Compile[{{z, _Complex}},
   Module[{c = 0. + 0. I, i = 0},
    While[i++ < 100,
     c = c^2 + z;
     If[2 < Abs@c, Return@i]];
    Return@-1]];
DistributeDefinitions[hc];

time = AbsoluteTiming[
   result = ParallelTable[ hc[N@x + N@y I],
      {y, yMin, yMax, dy}, {x, xMin, xMax, dx}];][[1]];

time0/time

77.13116

ちなみに、この例ではコンパイル時に「CompilationTarget -> “C”」を付けてもあまり効果は無いようです。

CUDA無しでのリアルタイム描画

この方法を使ってリアルタイムに描画してみます。スケールや反復の上限もリアルタイムに変更できるように、少し修正します。

hc2 = Compile[{{xIndex, _Integer}, {yIndex, _Integer},
    {maxSteps, _Integer},
    {width, _Integer}, { height, _Integer},
    {xMin, _Real}, {xMax, _Real},
    {yMin, _Real}, {yMax, _Real}},
   With[{z = xMin + (xMax - xMin) xIndex/width
       + (yMin + (yMax - yMin) yIndex/height) I},
    Module[{c = 0. + 0. I, i = 0},
     While[i++ < maxSteps,
      c = c^2 + z;
      If[4 < Abs@c, Return@i]];
     Return@-1]],
   CompilationTarget -> "C"];
DistributeDefinitions[hc2];

width = 300; height = 300;
Manipulate[
 With[{
   xMin = center[[1]] - 10^scale, xMax = center[[1]] + 10^scale,
   yMin = center[[2]] - 10^scale, yMax = center[[2]] + 10^scale,
   maxSteps = IntegerPart[10^logMaxSteps]},
  result = ParallelTable[
    hc2[xIndex, yIndex, maxSteps, width, height, xMin, xMax, yMin, yMax],
    {yIndex, 1, height}, {xIndex, 1, width}];
  ArrayPlot[result, DataRange -> {{xMin, xMax}, {yMin, yMax}},
   DataReversed -> True, ColorFunction -> "Rainbow"]],
 {{center, {-0.5, 0}}, Locator},
 {{scale, 0, "Scale"}, -13, 1},
 {{logMaxSteps, 2, "Log[maxSteps]"}, 1, 3}]

実行結果は上に掲載したものと同じです。速いCPUがある場合はこれで十分なのかもしれませんが、私の環境では、「一応動くがストレスフル」という感じでした。