今回はウラムの螺旋についてです。 
ドラマ「危険なビーナス」を見ていて
気になってしまったので
これを作るプログラムを
作成してみようと思います。

解説動画はこちら



まずウラムの螺旋とは
数学者のウラムさんが発見と言うか作った
素数をとある法則に従って並べた図形です。


1から渦巻状に数値を配置する際に
素数だけ残すと、対角線状に素数が残ります。

まずは1から数値を並べるところを
考えていきましょう。

def uram(N=7):
    # 箱を用意
    data = [ [0 for n1 in range(N)] for n2 in range(N)]

    # 入れる値を生成
    nums,tmp,c = [],[],1
    for num in range(1,N*N+1):
        tmp.append(num)
        if num==c*c:
            nums.append(tmp)
            tmp,c=[],c+1
            #c+=1

    # 起点位置決め
    c_y,c_x = N//2,N//2 if N%2==0 else N//2+1

    # 値の代入
    for i,tmp in enumerate(nums,start=1):
        for j,n in enumerate(tmp):
            sub = 1 if i%2==0 else -1
            if j==0:
                c_y,c_x=c_y,c_x + sub*1
            elif j>len(tmp)//2:
                c_y,c_x=c_y,c_x - sub*1
            else:
                c_y,c_x=c_y - sub*1,c_x
            data[c_y][c_x]=n
    return data

data = uram(5)
for row in data:
    for t in row:
        print('{0:}'.format(t),end='\t')
    print()

最初にN*Nの箱を作ります。

次に入れる数値の生成ですね。
こんな感じでN*Nの数値のリストを作ります。

N = 5
nums,tmp,c = [],[],1
for num in range(1,N*N+1):
    tmp.append(num)
    if num==c*c:
        nums.append(tmp)
        tmp,c=[],c+1
print(nums)

[[1],
[2, 3, 4],
[5, 6, 7, 8, 9],
[10, 11, 12, 13, 14, 15, 16],
[17, 18, 19, 20, 21, 22, 23, 24, 25]]

ここでは螺旋状に数値を入れる際に
折り返しがあるので
入れる数値のリストの半分のところで
折り返すようにしています。

Nが偶数の時はその前の数値から
右に一つずれ、そこから上
真ん中のところで左に折り返し

Nが奇数の時はその前の数値から
左に一つずれ、そこから下
真ん中のところで右に折り返す
ように数値を入れていくと
螺旋状に数値が並んだ
リストが生成できます。

次に素数だけを
残すようにしてみましょう。

sympyのisprimeと言うメソッドで
素数判定ができます。

これを使って素数の部分だけを
残すように変えてみます。

from sympy import isprime

def uram2(N=7):
    # 箱を用意
    data = [ [0 for n1 in range(N)] for n2 in range(N)]

    # 入れる値を生成
    nums,tmp,c = [],[],1
    for num in range(1,N*N+1):
        tmp.append(num)
        if num==c*c:
            nums.append(tmp)
            tmp=[]
            c+=1

    # 起点位置決め
    c_y,c_x = N//2,N//2 if N%2==0 else N//2+1

    # 値の代入
    for i,tmp in enumerate(nums,start=1):
        for j,n in enumerate(tmp):
            text = ' | '
            sub = 1 if i%2==0 else -1
            if j==0:
                c_y,c_x=c_y,c_x + sub*1
            elif  j>len(tmp)//2:
                text = ' - '
                c_y,c_x=c_y,c_x - sub*1
            else:
                c_y,c_x=c_y - sub*1,c_x
            # 素数判定
            if isprime(n):
                data[c_y][c_x]=n
            else:
                data[c_y][c_x]=text
    return data

data = uram2(5)
for row in data:
    for t in row:
        print('{0:}'.format(t),end='\t')
    print()

はいこれで素数だけ残りましたね。

次は塗りつぶしですが
数値を文字列に変えるだけでできます。


これだと使いづらいので
画像にするようにしてみましょう。

matplotlibではリストを投入すると
そのまま画像化できます。

素数を0,他は255を入れた
リストを作ってplotするだけです。

せっかくなので1000x1000
100万までの素数で
ウラムの螺旋を作ってみましょう。

from sympy import isprime
import matplotlib.pyplot as plt
%matplotlib inline

def uram4(N=7):
    # 箱を用意
    data = [ [255 for n1 in range(N)] for n2 in range(N)]

    # 入れる値を生成
    nums,tmp,c = [],[],1
    for num in range(1,N*N+1):
        tmp.append(num)
        if num==c*c:
            nums.append(tmp)
            tmp=[]
            c+=1

    # 起点位置決め
    c_y,c_x = N//2,N//2 if N%2==0 else N//2+1

    # 値の代入
    for i,tmp in enumerate(nums,start=1):
        for j,n in enumerate(tmp):
            sub = 1 if i%2==0 else -1
            if j==0:
                c_y,c_x=c_y,c_x + sub*1
            elif  j>len(tmp)//2:
                c_y,c_x=c_y,c_x - sub*1
            else:
                c_y,c_x=c_y - sub*1,c_x
            if isprime(n):
                data[c_y][c_x]=0

    return data

data = uram4(1000)
plt.figure(figsize=(10,10))
plt.imshow(data)
plt.axis("off")
plt.show()

メモリとCPUが強力であれば
もっと大きな螺旋も描けると思います。

これ以上はなかなかに痺れますね。

素数の法則は見つかっていないと
思いますが、なんかロマンがある図形ですね。

お暇な方は試してみてください。
それでは。

タグ ：

#プログラミング

#Python

#可視化

#ウラムの螺旋

#matplotlib