【例題 1】以下の問題を解け
\begin{align}
&nが2 \le n \le Nの整数の時 \ \ (N=10^7) \\
&f(n) = n^2+1 が素数になる個数を求めよ \\
\end{align}
pythonではsympyを使えば一行で書け速度もそこそこですがわずかに1分をオーバーしました。今回はこれを2次ふるい法と呼ばれるアルゴリズムを使って解いてみたいと思います。
from sympy import isprime
N = 10**7
print(sum([isprime(n**2 + 1) for n in range(N+1)]))
# 456361
2次ふるい法 （Quadratic sieving）
ここで紹介する「2次ふるい法」は一般に素因数分解で使われるものとは異なり、2次式で「エラトステネスのふるい」のように合成数を除いて行って素数だけが残るようにするアルゴリズムです。詳細は以下のサイトにありますが長文なので、そのさわりだけを紹介します。
ここで使われているアルゴリズムは以下の性質を利用しています。
- もし$f(n)=n^2+1$が素数ならばこれを$p$とすると、$f(pk \pm n) (kは任意の整数)$は$p$で割り切れる合成数である
\begin{align}
&【証明】 \\
f(pk \pm n) &= (pk \pm n)^2 + 1 \\
&= p^2k^2 \pm 2pkn + n^2 + 1 \\
&\equiv n^2 + 1 = p \equiv 0 \pmod{p} \\
\end{align}
これで一気に合成数をふるい落とすことができます。例えば、
\begin{align}
&f(2)=2^2+1 = 5 なので、n=2, p=5, \\
&f(5k \pm 2) はすべて5で割れる合成数となり \\
&これらを5で割れるだけ割ります \\
&下の例で、\\
&k=1の時は5-2=3と5+2=7を \\
&k=2の時は10-2=8と10+2=12をnとしたf(n)を \\
&5で割れるだけ割ります。この時これらをは \\
&素数でないと分かるので素数フラグをFalseにします。\\
\end{align}
f( 1 ) = 2
f( 2 ) = 5
f( 3 ) = 10 --> 2 : not prime
f( 4 ) = 17
*f( 5 ) = 26
f( 6 ) = 37
f( 7 ) = 50 --> 2 : not prime
f( 8 ) = 65 --> 13 : not prime
f( 9 ) = 82
*f( 10 ) = 101
f( 11 ) = 122
f( 12 ) = 145 --> 29 : not prime
実はこの処理は一石二鳥で以下の2つを同時に行っています。
- f(n)が合成数であることが分かる
- f(n)を既知の素数で割ることによって、新たな素数が出てくる（f(8)の13, f(12)の29）ので、これらをこの後の処理に利用できる
この$f(pk \pm n)$は$p$で割る処理を$pk \pm n \le N$で行い、これを$n \le N/2$まで行います。素数フラグがFalseにならなかった$f(n)$が素数であることが分かります。
これをpythonで実装したのが以下になります。実行時間はisprimeを使った場合の半分程度ですが、numbaを使えばその1/10程度になります。
from numba import njit
@njit
def solve(N):
f = [i**2+1 for i in range(N+1)]
primes = [1]*(N+1)
for n in range(1,N//2+1): # enough to check b< 1/2*N
if f[n] == 1: continue
p = pk = f[n]
#print(p, b)
while pk-n <= N:
while f[pk-n] % p == 0: f[pk-n] //= p
primes[pk-n] = 0
if pk+n <= N:
while f[pk+n] % p == 0: f[pk+n] //= p
primes[pk+n] = 0
pk += p
return sum(primes[2:])
N = 10**7
print(solve(N))
# 456361
【応用】n**2+1の素因数分解を求める
このプログラムをちょっと変更するだけで$n^2+1$の素因数分解を求めることが出来ます。
# Return factors of [k**2+1 for k=[0..N]]
def quadsieveFact(N):
f = [k**2+1 for k in range(N+1)]
fcts = [dict() for _ in range(N+1)]
for n in range(1,N+1):
if f[n] == 1: continue
p = pk = f[n]
fcts[n][p] = 1 # prime number
while True:
eflag = True
for d in [-n,n]:
if pk+d > N: continue
i, eflag = 0, False
while f[pk+d] % p == 0:
f[pk+d] //= p
i += 1
if i > 0:
fcts[pk+d][p] = i # prime factors
if eflag: break
pk += p
return fcts
print(quadsieveFact(10))
# [{}, {2: 1}, {5: 1}, {2: 1, 5: 1}, {17: 1}, {2: 1, 13: 1}, {37: 1}, {2: 1, 5: 2}, {5: 1, 13: 1}, {2: 1, 41: 1}, {101: 1}]
他の2次式への応用
このアルゴリズムは以下の2次式にもほぼそのまま利用することが出来ます。
\begin{align}
f(n)&=4n^2+1 \\
f(n)&=2n^2+1 \\
f(n)&=2n^2-1
\end{align}
（開発環境：Google Colab）
このアルゴリズムは以下のProject Eulerの問題を解くのに役に立ちます。
- Problem 216: Investigating the primality of numbers of the form $2n^2-1$
- Problem 659: Largest prime
- Problem 446: Retractions B
（開発環境：Google Colab）