このブログ記事は、2023年8月に部内向けに書いたものを外部公開向けに変更したものです。

はじめに

こんにちは。koukawa_ppです。

さて、螺旋はお好きですか？
ここでいう螺旋とは、二次元平面における螺旋です。

こういうタイプのもです。
こういった螺旋について、今回は深く掘り下げていこうというお話です。

なお、今回取り上げる螺旋には、以下のようなものがあります。

• アルキメデスの螺旋
• サックスの螺旋
• ウラムの螺旋

これらについて、特にPythonのコードも交えながら、色々螺旋を描いてみましたので、共有したいと思います。

なお、今回のコードや画像ファイルは全てこちらのGitHubリポジトリに入れてあります。
ぜひご覧ください。

実行環境

• Windows 10 Home 64bit
• Python 3.9.13
• Visual Studio Code

アルキメデスの螺旋を描く

そもそも、なぜ螺旋を描きたいとなったのかと言いますと、まずどこかで、「素数螺旋」というものがある、ということを耳にしたためです。
とは言えその時は螺旋について詳しく知らなかったので、とりあえず何も考えずにオーソドックスな螺旋を描こうという話になりました。

ここでいうオーソドックスな螺旋こそが、アルキメデスの螺旋です。
すなわち、極座標表示すれば、$r=a\theta$ と表される曲線です。
最初は無知故適当に見えればよいかという形で描画していたのですが、あとから調べてみるとこちら や数Ⅲの教科書などに載っていたので、これに則って描画することにしました。
この参考のページには、コンパスと定規で描く近似法を示していただいています。

次の問題です。
$r=a\theta$ という式を得られたわけですから、
あとはこれを $x=r\cos x,y=r\sin x$ という形で変数変換すれば、
これを描画することはお手の物です。
しかし、このまま何も考えずに $\frac{\pi }{180}$ ラジアンずつなどと描画すると、こうなります。

そうです、直径が大きくなればなるほど、
もっと言ってしまえば $\theta$ が大きくなればなるほど、
点と次の点の間があまりにも空きすぎてしまい、
綺麗な螺旋を描けなくなってしまうのです。

これを改善するために、螺旋上に等間隔に点を打つ必要が出てきました。
これを実現する方法は簡単です。

1. ある基準の長さ $d$ を設定する。
2. ある時点で最後に打った点と基準線のなす角度を $\alpha$ として、次に打つべき点と基準線のなす角度を $\beta$ とするとき、螺旋の $\alpha \leqq \theta \leqq \beta$ の部分における曲線の長さが、$d$ と一致すればよい。
3. 2.から、適切な角度 $\beta$ を求める。

発想としてはそこまで難しくありません。
しかし、行うは難しです。
螺旋の長さなんてそう簡単に求められるわけがありません。
……いや、求めます。

調べてみると、極座標における曲線の長さを求める公式がありました。
$\theta$ が ${\theta }_0$ から ${\theta }_1$ までにおける曲線の長さ $L$ は、

$$L={\int }_{{\theta }_0}^{{\theta }_1}\sqrt{r^2+{\left(\frac{dr}{d\theta }\right)}^{2}d\theta }$$

で与えられるとのことです。詳しくは こちら をご覧ください。いつもお世話になっております。
なのでこの曲線の長さは、

$$L={\int }_{\alpha }^{\beta }\sqrt{r^2+{\left(\frac{dr}{d\theta }\right)}^2}d\theta ={\int }_{\alpha }^{\beta }\sqrt{(a\theta {)}^2+a^2}d\theta =a{\int }_{\alpha }^{\beta }\sqrt{{\theta }^2+1}d\theta$$

と表されます。
次に こちら によれば

$$\int \sqrt{x^2+a^2}dx=\frac{1}{2}(x\sqrt{x^2+a^2}+a^2\log x+\sqrt{x^2+a^2})$$

とのことですので、つまり

$$a{\int }_{\alpha }^{\beta }\sqrt{{\theta }^2+1}d\theta =\frac{a}{2}(\beta \sqrt{{\beta }^2+1}+\log (\beta +\sqrt{{\beta }^2+1})-\alpha \sqrt{{\alpha }^2+1}-\log (\alpha +\sqrt{{\alpha }^2+1})$$

というわけです。これが $d$ と等しくなる $\beta$ を見つけられれば、晴れて螺旋上に等間隔で点を打つことができます。


では、ある $a,\alpha$ と $d$ が与えられたとき、どのようにして $\beta$ を求めれば良いのでしょうか。
直感から、以下の方程式を解けばいいことが分かります。

$$\frac{a}{2}(\beta \sqrt{{\beta }^2+1}+\log (\beta +\sqrt{{\beta }^2+1})-\alpha \sqrt{{\alpha }^2+1}-\log (\alpha +\sqrt{{\alpha }^2+1})=d$$

ここは、文明の利器、scipy を頼ることにします。
scipy.optimize.fsolve() 関数は、第一引数に $f(x)=0$ の解を求めたい関数 $f(x)$、第二引数にその解の値に近い $x$ の値を与えることで、その方程式の近似解を求めることができるという優れものです。

ここで、$a$ と $d$ については、draw_spiral() 関数の引数の radius と distance をそれぞれ与えればよく、また $\alpha$ は、前回の $\theta$ の値をそのまま使えばよいです。
では、fsolve() 関数の第二引数には何を与えるかということですが、$\beta$ に近い値を与えればよいということなので、ここは素直に $\theta$ の値を与えればよいということになります。

これらを実現するために使っているコードというのが、

draw_spiral.pyfrom scipy import optimize

def _length(x):
    return a * (x * np.sqrt(x * x + 1) - theta * np.sqrt(theta * theta + 1) + np.log(x + np.sqrt(x * x + 1)) - np.log(theta + np.sqrt(theta * theta + 1))) / 2 - d

def draw_spiral(
        draw_index: set[int] = None,
        width: int = 1000,
        height: int = 1000,
        radius: float = 5.0,
        distance: float = 1.0,
        num_of_cycle: int = 10,
        file_name: str = None
    ) -> None:
    # 中略
    while theta < 2 * math.pi * num_of_cycle:
        # 中略
        theta = optimize.fsolve(_length, theta)[0]

というわけです。
また、どうも optimize.fsolve の第一引数にとる関数は、引数に ndarray を取る必要があるらしいので、_length() 関数にはやむなく numpy の関数を用いて、optimize.fsolve() 関数の0番目を theta に代入するという形になっています。
（最初 math を使ってたらwarningが出てきた……）
これを使って、等間隔に点を打つと以下のような感じになります。

点と次の点の間が等しくなったような気がしませんか？
これでうまく機能しました。


次に、draw_spiral() 関数に draw_index という引数があると思いますが、ここは初期値として None になっています。
ここに set を渡すと、指定されたインデックスだけの点が表示されるようになります。
たとえば、draw_index に 0, 4, 6, 9 と代入されているとすると、0, 4, 6, 9番目の点だけが表示されるようになります。
これでお分かりのように、点は0番目から開始します。
set にランダムに数字を追加すると、以下のようなランダムな螺旋が描けます。

なお、最後に ImageOps.flip() を実行していると思います。
これは、僕が左上が原点だと思っていたところ、pillow によると左下が原点であるということらしいので、つじつまを合わせるために上下反転させたということです。
GitHub の方に API リファレンスを載せておきますので、draw_spiral.py でいろいろ遊びたくなりましたらそちらをご覧ください。


さて、これで螺旋をランダムに描くと、以下のような感じになります。

まあ綺麗っちゃ綺麗かもしれませんが、
ランダムさが生んだ綺麗さですよね。

無理数も無理やり螺旋にしてみましょう。
たとえば、 $\sqrt{2}$ は、1.41421356...と続きますが、0から始めて各位の数字を足していくという数列を考えてみましょう。
つまり、$\sqrt{2}$ の小数第 $n$ 位を $b_n$ として、今回 set に追加する数字の内小さい方から $n$ 番目を $a_n$ とするなら、

$$a_0=0,a_1=1,a_{n+1}=a_n+b_n$$

とすれば良いわけです。
こうすると、$a_0=0,a_1=1,a_2=5,a_3=6$ と求められていきます。
時々 $b_n=0$ となりますが、それは $a_{n+1}=a_n$ として処理することにして、そうすると無理数から作る螺旋も一意に決定することができます。
こうして $\sqrt{2}$ から作った螺旋は以下のようになります。

$\sqrt{3}$ は以下の通り。

$\sqrt{5}$ は以下の通り。

良さげですね。

では、待ちに待った素数螺旋を描いてみましょう。
素数のところだけ表示させた螺旋は以下のようになります。

ん？規則、的？

……まあ確かに途切れてるところは似ているし？
周期的に途切れてるし？
まあ綺麗に見えるし？
素数螺旋ってこういうものなのかぁ～（棒）

ウラムの螺旋

koukawa_ppはこんらんした！
そすうの　らせんが　そこまで　うまく　みえない！

さて、こんなところで終わっていられません。
こういった模様になってくると、人によって評価が分かれるところ。
本来大多数の人が規則的だと思えなければ、あとは定量的に測るなどして、（無理やり）規則性を発見できなければ、それは単なる見た目だろ、という話になってしまいます。
多くの人に認められた黄金比でさえも、確かに美しくは見えるがなぜかはわからない、といったように、あるものが美しい理由というのが科学的に解明されたようなものは決して多いわけではありませんが、今回のこれについては大多数が美しいと感じるかと言われると、多分ハテナマークを浮かべる人が一定数いると思ったわけです。

そこでネットをいろいろ漁ってみると、「ウラムの螺旋」というものを見つけました。

こちら によると、ウラムの螺旋は、中心から螺旋を描くように描画されます。
数字は格子点上にあるようにも解釈できます。
このことから、ウラムの螺旋もPythonで描いてみようということになりました。

そうして作ったのが、draw_ulam_spiral.py です。
これは、アルゴリズムの部分が若干ややこしくなっています。

このように螺旋を描画していくのが、ウラムの螺旋です。
まず1のところを座標(0,0)として、2のところを(1,0), 4のところを(0,1)とします。
このとき、ある方向に進んで、曲がって、進んで、曲がってを繰り返すと思いますが、
この曲がる点というのは、上下左右に4つあります。

そして、各々の点について、次に曲がる点は青字のところになります。

つまり、

としておいて、これらの座標に到達したときは、
右下の点の座標は右下に、
右上の点の座標は右上に、
左上の点の座標は左上に、
左下の点の座標は左下にずらせばよいわけです。

つまり2(1, 1)にたどり着いた時は、

とすればよいのです。

さて、これらのことを使ってアルゴリズムを見ていきましょう。
まず、いくつか変数を定義しています。

draw_ulam_spiral.pyx, y = 0, 0
lrx, lry = 1, 0
urx, ury = 1, 1
ulx, uly = -1, 1
llx, lly = -1, -1
dx, dy = 1, 0

x, yは、現在いる座標を示しています。
lrx, lryは右下の座標、urx, uryは右上の座標、
ulx, ulyは左上の座標、llx, llyは左下の座標です。
また、dx, dyは次x方向およびy方向にどれだけ進むかを表しています。

(x, y)を描画したら、次の座標を求めるために、こうします。

draw_ulam_spiral.pyx += dx
y += dy

次に、今いる座標が右下、右上、左上、左下のいずれかの場合は、dxおよびdy、また右下、右上、左上、左下のうち現在いる場所の座標を次のものに変更させる必要があります。

まず右下にいる場合、次は上の方向に進む必要があるので、

draw_ulam_spiral.pydx = 0
dy = 1

とすればよいことが分かります。
次に、右下の座標を一つ右下にずらすので、

draw_ulam_spiral.pylrx += 1
lry -= 1

とすればよいことが分かります。これが、

draw_ulam_spiral.pyif x == lrx and y == lry:
    dx = 0
    dy = 1
    lrx += 1
    lry -= 1

となっている理由です。
これを他にも右上、左上、左下についてelifでくっつければOKです。
本当はdictとかで取得すればよかったと、今さらながら気がつきました。

あとこのまま行くと、点と点の間の長さを決められなくなるので、実際に描画する際は、xおよびyにdistanceを掛けることで最終的な座標を求めています。

これを使って実際にウラムの螺旋を描くと、こうなります。

……まあ、螺旋も何もないですよね。
なぜならこれは若干厄介な方法で格子点を描画しているのに過ぎないのですから。
この螺旋の本領を発揮するのは、先ほどアルキメデスの螺旋でもやったように、部分的に描画しない点を設定することです。

この描画する点を指定するのは、先ほど同様 draw_index です。
先ほどと同じように指定できます。
ただここで一つ違うのが、先ほど点は0番目からカウントと言っていましたが、ここでは base_val 番目からカウントを行います。
それゆえ、draw_spiral() 関数の引数で指定をしない限り、1番目からのカウントになることに注意してください。
ではこれで素数螺旋を描画すると、こうなります。

おっ、きれいですね。
何かx状に模様が入っています。
ただこれは、ウラムの螺旋だからこうなっただけかもしれません。
ランダムに螺旋を作ってみましょう。

確かに、先ほどのように整った螺旋には見えません。
$\sqrt{2},\sqrt{3},\sqrt{5}$ で、さっきと同じルールで作った螺旋も見てみましょう。

特段規則性は見られませんね。
素数でないと、今のところはああいった模様は描けないみたいです。

あともう一つ、先ほどの記事によれば、こんなこともあるらしいですね。

そしたらやってみるしかないですよね。
その実、そのために base_val を設定したのもあります。
というわけで、41で見てみると、以下のようになります。

確かに中心部にかなりはっきりと斜め線が見えている……。
これはオイラー素数の話に繋がってくるみたいですが、
ここではその話はしません。

ウラムの螺旋と最小公倍数

ウラムの螺旋については、これともう一つ面白いものを見つけたので、そちらも紹介していきましょう。

最初のセッティング通り、base_val を初期値1のままにして、以下のようなプログラムを実行しました。

draw_rule_ulam_spiral.pyfrom draw_ulam_spiral import draw_spiral

hmp = 150000
for rule in range(1, 1001):
    draw_index = set([i for i in range(hmp) if i % rule == 0])
    draw_spiral(draw_index=draw_index, num_all_points=hmp, file_name='rule_ulam_spiral_{}.jpg'.format(rule))

ここでのdraw_ulam_spiralは、draw_ulam_spiral.pyのことです。
こうすると、rule個おきに点を打ったウラムの螺旋の画像ファイルをゲットできます。
なお、$1\leqq rule\leqq 1000$ です。
これについて、ruleが2, 3, 6の場合を見てみましょう。

同じ雰囲気を感じませんか？
そうです、思い付けば「当然だよね」となりますが、よく考えてみると、このruleが6の場合のウラムの螺旋は、
ruleが2の場合の螺旋とruleが3の場合の螺旋を比べて、どちらにも描画されている点からなる螺旋であるということです。
これは、6が2と3の最小公倍数であることから求まります。

24の場合も同様で、8と3の螺旋のうち、重なっている点のところだけが描画されます。
下の図は、上からruleが3, 8, 24のウラムの螺旋です。

8や24の場合、右上に斜め線が見えるという形なんですよね。
実は16や32の場合は、左下に斜め線が見えます。
下の図です。

かなり不思議です。
ウラムの螺旋は今回初めて知りましたが、面白かったです。

サックスの螺旋

しかし、このまま行くと、アルキメデスの螺旋において、等間隔に点を打つということを実現するため、難しい積分の話までしていたことが無駄になってしまいます。
何とかしてアルキメデスの螺旋を使った、綺麗な素数螺旋を描くことはできないものかと調べてみました。
先ほどの記事では同時に、「サックスの螺旋」というものがあることも紹介されていました。
書き方としては、こちらのページ が詳しいと思います。
中心には0を置き、あとは周回ごとに平方数を配置するというものです。

アルゴリズム的にはアルキメデスの螺旋とさほど変わりませんが、やることは極端に複雑になっていきます。

まず、素数の処理を簡単にするために、0と1については最初に配置しておくことにします。
これを行っているのが、以下の部分です。

draw_sacks_spiral.pyd = _local_length(2 * math.pi, 4 * math.pi) / 3
if all_show:
    draw.point((width / 2, height / 2), fill=(255, 255, 255))
    num_of_points += 1

    r = radius * 2 * math.pi
    draw.point((width / 2 + r * math.cos(2 * math.pi), height / 2 + r * math.sin(2 * math.pi)), fill=(255, 255, 255))
    theta = optimize.fsolve(_length, theta)[0]

これはall_showがTrueの場合にのみ実行されるので、素数のみを見せる場合、すなわちall_showがFalseの場合は実行されません。

まず最初のdraw.pointが、0の点を描画する処理です。
その次のdraw.pointの周辺が、1の点を描画する処理です。

次に、1, 4, 9, 16, ...を配置するのは、それぞれthetaが2π, 4π, 6π, 8π, ...となるときです。
そして、2πから4πの間に2と3、4πから6πの間に5から8、6πから8πの間に10から15……といったように、それぞれ等間隔に配置する必要があります。

これは言いかえると、
2πから4πの間を3等分、
4πから6πの間を5等分、
6πから8πの間を7等分、
……という感じだと言えます。
これを一般化すると、
2nπから2(n+1)πの間を2n+1等分と言えます。
このことを使って、thetaを求めていきます。

draw_sacks_spiral.pynext_i = 2
while i < num_of_cycle * num_of_cycle + 1:
    # 中略
    if i == next_i * next_i - 1:
        theta = 2 * next_i * math.pi
    elif i == next_i * next_i:
        d = _local_length(2 * next_i * math.pi, 2 * (next_i + 1) * math.pi) / (2 * next_i + 1)
        next_i += 1
        theta = optimize.fsolve(_length, theta)[0]
    else:
        theta = optimize.fsolve(_length, theta)[0]

まず、next_iとは、現在のiより大きい最小の平方数の正の平方根を与えます。
それゆえ、iを2で初期化するときに、next_iも2で初期化されています。
（2より大きい最小の平方数は4であり、この正の平方根は2であるため）

まず、iがnext_i * next_iより1小さいときは、次のiがnext_i * next_iであることが分かります。
それゆえ、誤差を無くすために、thetaに2 * next_i * math.piを代入しています。
なぜ2 * next_i * math.piなのかということですが、これは具体例を考えればよいです。
たとえば現在のiが3のとき、次の平方数は4であり、next_iは2であるので、iはnext_i * next_iより1小さいことが分かります。
4はthetaが4πのところに配置されるべきなので、2 * next_i * math.piは成り立ちます。
または、現在のiが8のとき、次の平方数は9であり、next_iは3であるので、iはnext_i * next_iより1小さいことが分かります。
9はthetaが6πのところに配置されるべきなので、2 * next_i * math.piは成り立ちます。

次に、iがnext_i * next_iと等しい場合を考えます。
このとき、iは平方数です。これより、次のいくつかについて、どの間隔で点を配置すべきかを求める必要があります。
先ほどの話より、現在のthetaは2 * next_i * math.piであることが分かります。
つまり、

先ほど記述したこのルールに則ると、
2 * next_i * math.pi から 2 * (next_i + 1) * math.pi の間を 2 * next_i + 1 等分した値を、dに代入する必要があるということになります。

これは_local_length 関数において、α=2 * next_i * math.pi, β=2 * (next_i + 1) * math.pi を代入し、これを 2 * next_i + 1 で割ったものと一致します。
これより、dは、

draw_sacks_spiral.pyd = _local_length(2 * next_i * math.pi, 2 * (next_i + 1) * math.pi) / (2 * next_i + 1)

という更新式が得られます。
次に、next_i * next_i より大きい最小の平方数は、next_i + 1 * next_i + 1 ですから、この正の平方根は、next_i + 1 です。
ゆえに、next_i を1増やす必要があります。
そのあとはいつものように、thetaを更新すればよいです。

iがこのいずれの条件も満たさない場合は、thetaを普通に更新すればよいです。

これらのことをまとめると、先ほどのようなコードが得られます。


他のことについては、アルキメデスの螺旋でやったこととほぼ同じです。
ただし、draw_indexを無くし、その代わりにall_show引数を追加しました。
ここでall_showをTrueにすると、全ての点が表示されます。
反面、all_showをFalseにすると、素数の点だけが表示されます。

さて、all_showをTrueにして実行してみると、以下のようになります。

ここで、素数に切り替えてみます。すると、

なんか素数が連なっている部分もありそうですね。
radiusを小さくして、より大きな範囲を表示すると、以下のようになります。

なんと！という感じです。
やっぱり、なんか綺麗ですよね。
よくできてます。

おわりに

本日は、様々な螺旋をご紹介しました。
何かお気に入りの螺旋(?)は見つかりましたか？
また、付属のプログラムもいろいろいじっていただければ幸いです。
最後までお読みいただき、ありがとうございました。