「ウラムの螺旋」を描く [日常のあれやこれや]
ここんとこ新しい光学素子の設計の最適化や組み立て設備のシーケンサがわりのRaspberry Piやなんかのソフトウェアを書く仕事をいくつか抱えている。どうしても大きなものになりがちで、一人てやってると煮詰まってしまいやすい。とくにここ数年歳のせいかスループットが落ちてるような気がする。
そんなことでついつい、ちょこちょこっとした全然関係ないものを気晴らしに書いてしまって、そっちが面白くなってしまうことがある。もっと真面目にやらないといけないと思うんだけど、集中力が落ちてしまってどうしようもない。買えるものなら集中力や思考力を買いたいものである。それが手に入るなら金に糸目はつけないんだけど。
それはまあいいとして（よくないけど）、そんな気晴らしの中でこないだのAnglePathを使えばすぐ描けると思って、「ウラムの螺旋」を描くMathematicaのコードを書いた。実は、ちょっと風邪をひいてしまって昨日今日と会社をお休みしてた。その間、退屈だったのでついつい遊んでしまっていた。
今回のMathematicaコードはRaspberry Piのバージョンでも動くのでRaspberry Pi3なんかを持っている人は試してみてほしい。しかし「ウラムの螺旋」って不思議だなあ．．．．
格子点をぐるぐる螺旋を描いてまわるAnglePathはすぐできる。ウラムの螺旋の描き方では1回転するごとにふたつ平方数が現れ、その次が曲がり角になる。曲がり角になる点$n$は一つ前$m=n-1$に関して$\lfloor \sqrt{m} \rfloor$が$\sqrt{m}$と一致する点だということになるなので、例えばまず、なんばんめの単位長さになるかという$n$を与えてそこが螺旋の曲がり角かどうかを返す関数は
MacBook Pro Late2013で約2秒かかった（平方数の判定やDiskを描くのとか最適化すべきところはあるけど、そこで頑張っても倍も速くはならない）。うーん、なんだこの斜め線は。奇数は45°方向に並ぶので斜めに並んでいるように見えるのは当たり前なのかもしれない。それとは別に縦横に空いた線（素数がないところ。原点の近くから左や下に伸びてる隙間とか）があったりする。縦横には偶数と奇数がかわりばんこに並ぶのでその方向に素数がないというのもなんとも不思議だ。
人間の目（と脳）はパターンを見つけたがる傾向があるので、無意識のうちに過大評価してる可能性もある。しかしリーマン予想が正しいとすれば素数の分布はランダムだ、ということになるんじゃなかったっけ？
素数って特別な数のように思ってたけどずいぶんたくさんあって、しかもこのくらいの範囲では数が大きくなっても全然減らないな。MathematicaでLogIntegral（素数定理による分布の近似関数${\rm li}(z)=\int_0^zdt/\log t$のMathematica実装）を数値評価したら$2^{170}$ぐらいでやっと100個にひとつにまで減る程度なので、$512^2=2^{18}$個ぐらいでは密度が下がったようには見えないんだな。いまさらだけど。
ところでよく考えたら、このやりかただと螺旋の全部の格子点に対して計算してる。ウラムの螺旋を描きたいだけなら、素数を与えたときにその螺旋上の座標を返すような関数を作ればもっと速いはずである。
でもその関数がどうなるのかすぐにはピンとこない。最初に書いた螺旋の曲がり角の点を手掛かりにすると
もちろんさっきのと分布は同じ（のはず）。こっちでは約0.2秒ほどでさっきより一桁速くなった（${\rm li}(512^2)/512^2 \approx0.1$ぐらいなのでもっともらしい）。ちなみにPrimePiは引数以下にある素数の（LogIntegralと違って厳密な）個数を返す組み込み関数である。これだとMathematicaのフロントエンドがDiskオブジェクトを描画する時間の方がずっと長くなる。
しかしそれ以前に座標を求める関数coord[]のIf文の入れ子がカッコ悪い。こういうぐちゃっとしたのって一晩経つだけでなにやってるのかわからなくなる。もうちょっとスッキリした書き方はないものか。螺旋上の値$n$に対して次の$m$ \begin{align} m &= \frac{\sqrt{4n-3}+1}{4} \nonumber \\ m &= \frac{\sqrt{n-1}}{2} \nonumber \\ m &= \frac{\sqrt{4n-3}-1}{4} \nonumber \\ m &= \frac{\sqrt{n}-1}{2} \nonumber \end{align} が整数になるとき、その$n$が曲がり角（4方向の対角線上）なのでもうすこし簡単に書けるはず。しかしこれ、式の形が微妙に対称性悪いな。もともとの螺旋の書き方が厳密には対称じゃないからか。しょうがないのか。
これで2k回転（$4096 \times4096$、素数の個数$\approx$1M個）までの素数を描いてみたけど、端っこの方でも斜め線のパターンが見える。不思議だ。速くなったのはいいんだけど、AnglePathは関係なくなってしまったな。
ところで、Wikipediaの最後にはもっと変というか不思議な絵がある。重複を数えた素因数の数に円の大きさを比例させると斜めのパターンがもっとはっきり現れる。さっきのprimeDisk[]のかわりに
これを最初に絵を描いたやりかたで描く。つまり
小さな正方形の格子パターンが見えるけど、それはビットマップとのモアレのせいで、ちゃんと合わせないといけなかった。でももうめんどくさい。
こっちは45°方向にはっきりと並んでいるように見える。しかも上下左右の三角の中身がパッとみて違うパターンに見える。しかもところどころに20°ぐらいとか70°ぐらいの方向の直線もあるように見える。不思議だ。
今の数論ではこんな簡単なパターンの説明もつけられないらしい。友愛数や双子素数なんていうような個別問題にかかずらってなくて、もっと整数全体の構造の研究に力を入れないといけないのではないか。
僕にはできないけど。
そんなことでついつい、ちょこちょこっとした全然関係ないものを気晴らしに書いてしまって、そっちが面白くなってしまうことがある。もっと真面目にやらないといけないと思うんだけど、集中力が落ちてしまってどうしようもない。買えるものなら集中力や思考力を買いたいものである。それが手に入るなら金に糸目はつけないんだけど。
それはまあいいとして（よくないけど）、そんな気晴らしの中でこないだのAnglePathを使えばすぐ描けると思って、「ウラムの螺旋」を描くMathematicaのコードを書いた。実は、ちょっと風邪をひいてしまって昨日今日と会社をお休みしてた。その間、退屈だったのでついつい遊んでしまっていた。
今回のMathematicaコードはRaspberry Piのバージョンでも動くのでRaspberry Pi3なんかを持っている人は試してみてほしい。しかし「ウラムの螺旋」って不思議だなあ．．．．
格子点をぐるぐる螺旋を描いてまわるAnglePathはすぐできる。ウラムの螺旋の描き方では1回転するごとにふたつ平方数が現れ、その次が曲がり角になる。曲がり角になる点$n$は一つ前$m=n-1$に関して$\lfloor \sqrt{m} \rfloor$が$\sqrt{m}$と一致する点だということになるなので、例えばまず、なんばんめの単位長さになるかという$n$を与えてそこが螺旋の曲がり角かどうかを返す関数は
cornerQ[n_] = (Floor[Sqrt[n]]^2 + 1 == n)
|| (Floor[Sqrt[n]]^2 + Floor[Sqrt[n]] + 1 == n);
などと書ける。これを使って90°曲がるか、あるいはまっすぐいくかを角度で返すような
angleAt[n_?cornerQ] := Pi/2 angleAt[_] := 0として、AnglePathの引数にしてやって、$n$が素数だったら丸印を返してそうでなければ何もしないという関数
primeDisk[{n_, pos_}] := Disk[pos] /; PrimeQ[n]
primeDisk[_] := Sequence[]
に渡すようにして、こんなふうに絵に描く。
Timing[Graphics[{Pink,
Map[primeDisk,
With[{mx = 512^2},
Transpose[{Range[mx], AnglePath[ParallelArray[angleAt, mx - 1]]}]]]},
ImageSize -> Large]
]
これで$512 \times 512$角の中を描いたのが
人間の目（と脳）はパターンを見つけたがる傾向があるので、無意識のうちに過大評価してる可能性もある。しかしリーマン予想が正しいとすれば素数の分布はランダムだ、ということになるんじゃなかったっけ？
素数って特別な数のように思ってたけどずいぶんたくさんあって、しかもこのくらいの範囲では数が大きくなっても全然減らないな。MathematicaでLogIntegral（素数定理による分布の近似関数${\rm li}(z)=\int_0^zdt/\log t$のMathematica実装）を数値評価したら$2^{170}$ぐらいでやっと100個にひとつにまで減る程度なので、$512^2=2^{18}$個ぐらいでは密度が下がったようには見えないんだな。いまさらだけど。
ところでよく考えたら、このやりかただと螺旋の全部の格子点に対して計算してる。ウラムの螺旋を描きたいだけなら、素数を与えたときにその螺旋上の座標を返すような関数を作ればもっと速いはずである。
でもその関数がどうなるのかすぐにはピンとこない。最初に書いた螺旋の曲がり角の点を手掛かりにすると
anchor[n_] = Floor[Sqrt[n - 1]];として、座標を返す関数
coord[n_] := Module[{anc, pos, rem},
anc = anchor[n];
rem = n - anc^2 - 1;
If[OddQ[anc], pos = (anc + 1)/2;
If[rem <= anc, {pos, 1 - pos + rem}, {pos - rem + anc, 1 - pos + anc}],
pos = anc/2;
If[rem <= anc, {-pos, pos - rem}, {rem - anc - pos, pos - anc}]]
]
を試行錯誤的にむりやり書いて、
Timing[Graphics[{Orange,
ParallelMap[Disk[coord[#]] &, Array[Prime, PrimePi[512^2]]]},
ImageSize -> Large]
]
として描いたのが
しかしそれ以前に座標を求める関数coord[]のIf文の入れ子がカッコ悪い。こういうぐちゃっとしたのって一晩経つだけでなにやってるのかわからなくなる。もうちょっとスッキリした書き方はないものか。螺旋上の値$n$に対して次の$m$ \begin{align} m &= \frac{\sqrt{4n-3}+1}{4} \nonumber \\ m &= \frac{\sqrt{n-1}}{2} \nonumber \\ m &= \frac{\sqrt{4n-3}-1}{4} \nonumber \\ m &= \frac{\sqrt{n}-1}{2} \nonumber \end{align} が整数になるとき、その$n$が曲がり角（4方向の対角線上）なのでもうすこし簡単に書けるはず。しかしこれ、式の形が微妙に対称性悪いな。もともとの螺旋の書き方が厳密には対称じゃないからか。しょうがないのか。
これで2k回転（$4096 \times4096$、素数の個数$\approx$1M個）までの素数を描いてみたけど、端っこの方でも斜め線のパターンが見える。不思議だ。速くなったのはいいんだけど、AnglePathは関係なくなってしまったな。
ところで、Wikipediaの最後にはもっと変というか不思議な絵がある。重複を数えた素因数の数に円の大きさを比例させると斜めのパターンがもっとはっきり現れる。さっきのprimeDisk[]のかわりに
Clear[primeAndCompositeDisk];
primeAndCompositeDisk[{n_, pos_}] := {Pink, Disk[pos, 0.3]} /; PrimeQ[n]
primeAndCompositeDisk[{1, _}] := Sequence[];
primeAndCompositeDisk[{n_, pos_}] := {Blue, Disk[pos, 0.12 (PrimeOmega[n] - 2)]}
というようなのを定義する。1は素数でも合成数でもない、となっているので無視してる。Diskの半径は描きながら潰れないように調整したので係数に意味はない。ここでPrimeOmega[]は重複を数えた素因数の数を返す組み込み関数である。ちなみに異なる素因数の（重複を同一視した）数を返す関数はPrimeNu[]である。ちなみにPrimeOmegaもPrimeNuも1を約数とはみなさない（1を含めるWikipediaの素数の記述とは違う。でも重複を数えた場合、つまりいくつの素数の積になるか、というときに1を含めるのっておかしいよな。だって1は何回重複しても1だし。だからMathematicaのほうがconsistentだわな。素因数分解を半群とみるとかそういう代数構造を考えるからかな）。これを最初に絵を描いたやりかたで描く。つまり
Timing[Graphics[
Map[primeAndCompositeDisk,
With[{mx = 300^2},
Transpose[{Range[mx], AnglePath[ParallelArray[angleAt, mx - 1]]}]]],
ImageSize -> Large]
]
を実行したのがこれ（Wikipediaにある絵の面積4倍で描いた）
こっちは45°方向にはっきりと並んでいるように見える。しかも上下左右の三角の中身がパッとみて違うパターンに見える。しかもところどころに20°ぐらいとか70°ぐらいの方向の直線もあるように見える。不思議だ。
今の数論ではこんな簡単なパターンの説明もつけられないらしい。友愛数や双子素数なんていうような個別問題にかかずらってなくて、もっと整数全体の構造の研究に力を入れないといけないのではないか。
僕にはできないけど。
2017-06-20 21:32
nice!(0)
コメント(0)
トラックバック(0)
コメント 0