ウラムの螺旋(Ulam Spiral)をPythonとSympyで500万まで描いてみる。

ウラムの螺旋というのがある。

https://ja.wikipedia.org/wiki/%E3%82%A6%E3%83%A9%E3%83%A0%E3%81%AE%E8%9E%BA%E6%97%8B

数字をらせん状に置いていき、素数のみ塗っていくと構造（水平線、垂直線、対角線）が見えるというもの。

昔から知っているがそういや自分で描いたことなかった。理由は螺旋に置くいい感じのアルゴリズムが思いつかなかったから…

いい感じじゃなくて愚直にやるなら簡単か、と思ってやってみた。しかもPythonならSympyがあるので素数の判定も簡単。

コードはこんな感じ。

import numpy as np

import matplotlib.pyplot as plt

from sympy import isprime

n = 2302

a =np.zeros((n, n - 1),dtype=np.int64)

#中央

i = n//2 -1

j = n//2 -1

# 1, 2だけあらかじめ入力

a[i, j] = 0

i += 1

a[i, j] = 1

# 螺旋を描く

num = 3

s = 1

for t in range(n//2-1):

    for k in range(s):

        j -= 1

        a[i, j] = isprime(num)

        num += 1

    s += 1

    for k in range(s):

        i -= 1

        a[i, j] = isprime(num)

        num += 1

    for k in range(s):

        j += 1

        a[i, j] = isprime(num)

        num += 1

    s += 1

    for k in range(s):

        i += 1

        a[i, j] = isprime(num)

        num += 1

# 最後の部分

s -= 1

for k in range(s):

    j -= 1

    a[i, j] = isprime(num)

    num += 1

plt.figure(figsize=(50,50))

plt.imshow(a.T, cmap="Greys")

print(num - 1)

では100万まで描いたとき。

500万まで描いたとき

対角線がよく見える。実は1億でも普通にできるが、ファイルサイズが大きくなりすぎるので今回はここまで。