N番目の素数を表す式をPythonで

はじめに

N番目の素数を表す式が存在するようです。
出典：https://en.wikipedia.org/wiki/Formula_for_primes

P_{n}=1+
\sum_{i=1}^{2^{n}}\left\lfloor
\left(\frac{n}{\sum_{j=1}^i\left\lfloor
\left(\cos\frac{(j-1)!+1}{j}\pi\right)^{2}
\right\rfloor}
\right)^{1/n}
\right\rfloor\\
P_1 = 2, P_2 = 3, P_3 = 5, P_4 = 7, ...

式内の階乗や三角関数から凄そうなオーラを感じます。
一方でWikipediaには効率的ではない旨の記載もあります。
この式はいったい何なのか、どう非効率なのかをPythonで再現して理解を深めていきます。
私はPythonも数学も得意でないので、文章中に怪しい箇所があれば優しく教えていただけると嬉しいです！

Pythonで書いてみる

数式には慣れていませんが、コードならまだ読める気がします。
N番目の素数を求める処理をPythonで書いていきます。

数式に頼らないコード

まずは数式に頼らない素直なコードを書いてみます。

def is_prime(x):
    # xが素数か判定します
    for t in range(2, x):
        if x % t == 0:
            return False
    return x > 1

def basic_nth_prime(n):
    # 2から順に素数か調べ、n番目の素数を返します
    prime_count = 0
    i = 1
    while prime_count < n:
        i += 1
        if is_prime(i):
            prime_count += 1
    return i

for i in range(1, 5):
    print(f'basic_nth_prime({i}) = {basic_nth_prime(i)}')
# output
# basic_nth_prime(1) = 2
# basic_nth_prime(2) = 3
# basic_nth_prime(3) = 5
# basic_nth_prime(4) = 7

2から順に素数か調べていき、素数であればカウンターを増やし、カウンターがnになれば返却するコードです。

数式を再現したコード

続いて数式を再現してみます。前のセクションで書いたコードは動作確認に使います。

import math

def formula_nth_prime(n):
    # 数式を再現してn番目の素数を返します
    large_sigma_sum = 0
    for i in range(1, 2 ** n + 1):
        small_sigma_sum = 0
        for j in range(1, i + 1):
            small_fraction = (math.factorial(j-1) + 1) / j
            small_sigma_sum += math.floor(math.cos(small_fraction * math.pi) ** 2)
        large_sigma_sum += math.floor(math.pow((n / small_sigma_sum), 1 / n))
    return 1 + large_sigma_sum

def test(n):
    p1 = basic_nth_prime(n)
    p2 = formula_nth_prime(n)
    result = 'success' if p1 == p2 else 'fail'
    print(f'basic_nth_prime({n})={p1}, formula_nth_prime(n)={p2}, {result}')

for i in range(1, 5):
    test(i)
# output
# basic_nth_prime(1)=2, formula_nth_prime(n)=2, success
# basic_nth_prime(2)=3, formula_nth_prime(n)=3, success
# basic_nth_prime(3)=5, formula_nth_prime(n)=5, success
# basic_nth_prime(4)=7, formula_nth_prime(n)=7, success

ちゃんと素数を計算できています。処理内容はよくわかりません。

数式コードを分割する

処理内容を理解していくため、数式コードの一部を関数a, b, cに切り出してみます。

def a(x):
    return (math.factorial(x-1) + 1) / x
def b(x):
    return math.floor(math.cos(x * math.pi) ** 2)
def c(x, y):
    return math.floor(math.pow((y / x), 1 / y))
def formula_nth_prime(n):
    large_sigma_sum = 0
    for i in range(1, 2 ** n + 1):
        small_sigma_sum = 0
        for j in range(1, i + 1):
            small_sigma_sum += b(a(j))
        large_sigma_sum += c(small_sigma_sum, n)
    return 1 + large_sigma_sum

関数a, b, cを理解していきます。
今回の数式に影響しない要素は考えず、一部の制約や条件を省略します。

関数a

a(x) = \frac{(x-1)!+1}{x}

この関数はxが1か素数のとき整数になり、そうでない場合は整数になりません。（割り切れません）

import math
def a(x):
    return (math.factorial(x-1) + 1) / x
for i in range(1, 9):
    print(f'a({i})={a(i)}')
# output
# a(1)=2.0
# a(2)=1.0
# a(3)=1.0
# a(4)=1.75
# a(5)=5.0
# a(6)=20.166666666666668
# a(7)=103.0
# a(8)=630.125

なぜこうなるの？と不思議に思いますが、ウィルソンの定理が活躍しているようです。
p が素数ならば (p − 1)! ≡ −1 (mod p) が成り立つ。 逆に、整数 p > 1 に対し、(p − 1)! ≡ −1 (mod p) ならば、p は素数である。
Wikipedia：https://ja.wikipedia.org/wiki/%E3%82%A6%E3%82%A3%E3%83%AB%E3%82%BD%E3%83%B3%E3%81%AE%E5%AE%9A%E7%90%86
この定理から、(x-1)!+1がxで割り切れたらxは素数です。
(実はx=1の時も割り切れてしまいます...)

関数b

b(x) = \lfloor(\cos x \pi)^{2}\rfloor

この関数はxが整数の時に1になり、そうでない場合は0になります。

import math
def b(x):
    return math.floor(math.cos(x * math.pi) ** 2)
print(b(1.0), b(1.75), b(5.0), b(20.16))
# 1 0 1 0

まずcos(x pi)を通すことで、整数であれば1か-1になり、二乗することで1にします。整数でなければ二乗して0以上1未満になります。
これをfloor関数にかけることで、0と1に二分します。
関数aとbを組み合わせることで、素数か1なら1, そうでないなら0を出力するようになります！

import math
def a(x):
    return (math.factorial(x-1) + 1) / x
def b(x):
    return math.floor(math.cos(x * math.pi) ** 2)
for i in range(1, 9):
    print(f'b(a({i}))={b(a(i))}')
# output
# b(a(1))=1
# b(a(2))=1
# b(a(3))=1
# b(a(4))=0
# b(a(5))=1
# b(a(6))=0
# b(a(7))=1
# b(a(8))=0

関数c

c(x) = \left\lfloor\left(\frac{y}{x}\right)^{1/y}\right\rfloor

この関数はxがy以下の時に1になり、そうでない場合は0になります。
※xは1以上です

import math
def c(x, y):
    return math.floor(math.pow((y / x), 1 / y))
for i in range(1, 7):
    print(f'c({i}, 3)={c(i, 3)}')
# output
# c(1, 3)=1
# c(2, 3)=1
# c(3, 3)=1
# c(4, 3)=0
# c(5, 3)=0
# c(6, 3)=0

まずyをxで割るので、xがy以下であれば1~yになり、xがyより大きければ0以上1未満になります。
これを1/y乗するので、xがy以下であれば1以上2未満になり、xがyより大きければ0以上1未満になります。
さらにfloor関数にかけることで、0と1に二分します。

	xがy以下	xがyより大きい
y/x	1以上y以下	0以上1未満
(y/x)^(1/y)	1以上2未満	0以上1未満
floor((y/x)^(1/y))	1	0

数式コードを読み直す

関数a, b, cをコメントで置き換えて数式コードを読み直します。

def formula_nth_prime(n):
    large_sigma_sum = 0
    for i in range(1, 2 ** n + 1):
        small_sigma_sum = 0
        for j in range(1, i + 1):
            # 右辺 : jが素数か1なら1, そうでないなら0
            small_sigma_sum += b(a(j))
        # 右辺 : small_sigma_sumがn以下なら1、そうでないなら0
        large_sigma_sum += c(small_sigma_sum, n)
    return 1 + large_sigma_sum

• 内側のjループでは、1からiの範囲で[素数か1]の整数を数えています。
• 外側のiループでは、1から2^nについて、ある数以下の[素数か1]の個数がn以下になる整数を数えています。
• 関数全体では、ある数以下の[素数か1]の個数がn個になる最大の整数 + 1を計算しています。
  • 言い換えると、ある数以下の素数の個数がn-1個になる最大の整数 + 1を計算しています。

ある数	ある数以下の素数の数
2	1
3	2
4	2
5	3
6	3
7	4

n=3の場合は、「ある数以下の素数の数がn-1=2個になる最大の整数」が4なので、これに1を足して5が計算されます。3番目の素数は5なので期待通りの結果です。
関数aの分類にて1が素数グループに含まれているため、素数グループの要素数が1ずれています。これを最後の1加算で調整しているようです。

数式を読み直す

P_{n}=1+
\sum_{i=1}^{2^{n}}\left\lfloor
\left(\frac{n}{\sum_{j=1}^i\left\lfloor
\left(\cos\frac{(j-1)!+1}{j}\pi\right)^{2}
\right\rfloor}
\right)^{1/n}
\right\rfloor\\
P_1 = 2, P_2 = 3, P_3 = 5, P_4 = 7, ...

Pythonで理解を深めてきたので、この数式も読めるようになっています。
当初感じていた階乗や三角関数のオーラも、その役割がわかります。
階乗は素数判定のため、三角関数は素数判定の結果を0と1に2分するためだったわけですね。

効率について理解を深める

ここからは効率の観点から理解を深めていきます。

数値計算の観点

コード内に巨大な数値や実数を使う箇所が複数あります。
これらの数値は通常の整数に比べて数値計算のコストや精度で課題が出てくるはずです。
nに合わせて影響が大きくなるので、テスト範囲を広げて試してみます。

for i in range(1, 11):
    test(i)

basic_nth_prime(1)=2, formula_nth_prime(n)=2, success
basic_nth_prime(2)=3, formula_nth_prime(n)=3, success
basic_nth_prime(3)=5, formula_nth_prime(n)=5, success
basic_nth_prime(4)=7, formula_nth_prime(n)=7, success
basic_nth_prime(5)=11, formula_nth_prime(n)=11, success
basic_nth_prime(6)=13, formula_nth_prime(n)=13, success
basic_nth_prime(7)=17, formula_nth_prime(n)=129, fail
Traceback (most recent call last):
    (省略)
    return (math.factorial(x-1) + 1) / x
OverflowError: integer division result too large for a float

7番目の素数を求める際に129を出力してテストに失敗しており、8番目の素数を求める際にerrorが起きています。
原因を探ってみます。

数値計算の課題1

7番目で129が出ている原因は、素数判定がうまく動かず2**7+1=129が計算されているからだと考えられます。
数式の素数判定だけを取り出して素直な実装と比べてみます。

print([i for i in range(2, 30) if is_prime(i)])
print([i for i in range(2, 30) if b(a(i)) == 1])
# [2, 3, 5, 7, 11, 13, 17, 19, 23, 29]
# [2, 3, 5, 7, 11, 13]

どうやら7番目以降の素数が正しく計算できていないようです。
関数a, bの動作を確認します。

for i in range(11, 20):
    print(i, a(i), b(a(i)))
# 11 329891.0 1
# 12 3326400.0833333335 0
# 13 36846277.0 1
# 14 444787200.0714286 0
# 15 5811886080.066667 0
# 16 81729648000.0625 0
# 17 1230752346353.0 0
# 18 19760412672000.055 0
# 19 336967037143579.0 0

関数aの結果は出力の中央列です。関数aに素数を与えて整数の出力が期待されます。検証の範囲では問題なさそうです。
関数bの結果が出力の右列です。関数bに整数を与えて1の出力が期待されます。i=17, 19の時に1を出力できていません。
関数bの中身を模してcosにa(i)*piを与えてみます。

for i in [11, 13, 17]:
    print(i, a(i), math.cos(a(i) * math.pi))
# 11 329891.0 - 1.0
# 13 36846277.0 - 1.0
# 17 1230752346353.0 - 0.9999999789777547

i=11, 13の時は-1.0を計算できていましたが、i=17の結果は-0.9999...となっています。
数値計算の精度に影響が出たようです。

数値計算の課題2

8番目の素数を求める際のerrorもみてみます。

Traceback (most recent call last):
    (省略)
    return (math.factorial(x-1) + 1) / x
OverflowError: integer division result too large for a float

階乗由来の巨大な数値を割って素数判定する際に、floatの範囲を超えて数値計算に影響が出たようです。

数値計算の影響から理解を深める

計算補助のモジュールを導入することで影響を軽減できますが、いっそ数式に従わず通常の整数で完結するように書き換えても良さそうです。
例えば、a, b関数を素直な実装のis_prime関数に乗り換える。
他にも怪しい箇所はありますが、ここまでの確認で「計算過程で巨大な数値や実数を扱うため非効率」とは言えそうです。

数値計算を改善(追記)

関数a内の除算でFractionを使い、関数bの整数判定を書き換えたところ、10番目の素数までは正しく求められたとコメントをいただきました！
https://qiita.com/kagasan/items/cf8dfd6b10448c356795#comment-1300198ba0651775b998
Fractionを使うことで素数判定b(a(x))内の実数処理が消えるため、精度の問題が解決されます。

# 除算にFractionを使う
from fractions import Fraction
def a(x):
    # return (math.factorial(x-1) + 1) / x
    return Fraction(math.factorial(x-1) + 1) / x
def b(x):
    # return math.floor(math.cos(x * math.pi) ** 2)
    return 1 if x == math.floor(x) else 0

関数bが元の姿から結構変わった点には一旦目をつむり、何番目の素数まで正しく計算できるか手元で試してみました。
後ろの素数を求めるほど計算時間がかかるようになり、12番目の素数までは正しく求められましたが、13番目の計算の途中で中断しました。
数値計算の精度だけでなく、計算時間にも難がありますね。

アルゴリズムの観点

数の扱いだけでなくアルゴリズムの観点でも見てみます。
数式を再現したコードを眺めると3つの処理が入れ子になっているように見えます。

• 素数の候補の区間を調整する処理（数式で言うところの1つ目のシグマ, コードで言うところのiループ）
  • ループ回数が2のn乗
• 区間内の素数の数を調べる処理（数式で言うところの2つ目のシグマ, コードで言うところのjループ）
  • ループ回数は最大で2のn乗
• 素数判定の処理(数式で言うところの２つ目のシグマの中身, コードで言うところのa, b関数)
  • 階乗では最大で2のn乗-1個の要素を掛け合わせる

上記の入れ子の処理は何らかの方法で計算コストを下げられそうです。
例えば次のようなアプローチが思い浮かびます。

• 3つ目の素数判定の処理結果を再利用して計算コストを下げる
• 2つ目の区間内を調査する処理で、過去の区間の処理結果を再利用して計算コストを下げる
• 1つ目の区間を調整する処理で、処理2の単調性を活かして計算コストを下げる
• そもそも区間を決め打ちせずに前から素直に素数判定する

他にもアプローチはありますが「数式をそのまま再現するとアルゴリズムが冗長で非効率」と言えそうです。

おわりに

興味深いN番目の素数を表す式についてみてきました。
数式にはウィルソンの定理に由来する素数判定の処理と、素数候補の区間を表現するためのシグマ、条件分けのための三角関数や根号やガウス記号が含まれています。
計算機に解かせることは可能ですが、数値計算の精度やコスト, アルゴリズムの観点からあまり効率的ではないようです。
現在の計算機環境でN番目の素数を計算したい場合は、他のアプローチを選ぶとよさそうです。