仮想計算機構

IT業界と無縁な派遣社員のブログ

Twitterでみかけた素数の一般項を計算してみる

Python

面白そうだったので計算してみました。

プログラム

愚直に実装しました。count は for 文の中の計算が何回行われたかを記録します。

from math import floor,cos,factorial,modf,pi

def p(n):
    count = 0
    a = 1
    for m in range(1,2**n+1):
        s = 0
        for k in range(1,m+1):
            s += floor(cos((factorial(k-1)+1)/k*pi)**2)
            count += 1
        a += floor((n/s)**(1/n))
    return [a,count]

for n in range(1,10):
    pn,c = p(n)
    print("{}:{}".format(pn,c))

実行結果

素数が順番に出力されると思いきや13から129に飛びました。えー

> python .\main.py
2:3
3:10
5:36
7:136
11:528
13:2080
129:8256
Traceback (most recent call last):
  File ".\main.py", line 19, in <module>
    pn,c = p(n)
  File ".\main.py", line 13, in p
    s += floor(cos((factorial(k-1)+1)/k*pi)**2)
OverflowError: integer division result too large for a float

階乗の値がでかすぎると怒られています。整数同士の加減計算であれば問題ないのですが割り算やコサインの部分が面倒です。というかコサインの計算がアレなせいで誤った出力になったのではないかと思われます。上述のプログラムのうち、コサインの部分だけ以下で置き換えてみます。

s += floor(cos(modf((factorial(k-1)+1)/k/2)[0]*2*pi)**2)

実行結果

> python .\main.py
2:3
3:10
5:36
7:136
11:528
13:2080
17:8256
Traceback (most recent call last):
  File ".\main.py", line 22, in <module>
    pn,c = p(n)
  File ".\main.py", line 11, in p
    s += floor(cos(modf((factorial(k-1)+1)/k/2)[0]*2*pi)**2)
OverflowError: integer division result too large for a float

13の次が17となりましたが、依然としてエラーは解決できていません。

プログラム

最初のプログラムの反省点を生かし、コサインが周期関数であることを利用してプログラムを書き換えます。剰余記号を使わず剰余を計算するところがキモでしょうか。

from math import floor,cos,factorial,pi

def p(n):
    count = 0
    a = 1
    for m in range(1,2**n+1):
        s = 0
        for k in range(1,m+1):
            b = factorial(k-1)+1
            r = (b-k*(b//k))/k
            s += floor(cos(r*2*pi)**2)
            count += 1
        a += floor((n/s)**(1/n))
    return [a,count]

for n in range(1,10):
    pn,c = p(n)
    print("{}:{}".format(pn,c))

実行結果

> python .\main.py
2:3
5:36
7:136
11:528
13:2080
17:8256
23:131328
29:524800
31:2098176

問題なさそうです。計算量の面から見て実用的ではないですが、ツイート主さんのおっしゃるようにロマンがあります。もちろんそれを計算機でやったところでロマンなんか雲散霧消してしまうわけですが。

おまけ

ウィルソンの定理 - Wikipedia
こちらの定理が元ネタのようです。

プログラム
from math import factorial

def p(n):
    return factorial(n-1)%n==n-1

for n in range(2,100):
    if p(n):
        print(n)
実行結果
> python .\wilson.py
2
3
5
7
11
13
17
19
23
29
31
37
41
43
47
53
59
61
67
71
73
79
83
89
97