与えられた$n$について、素数かどうかを判定したい。
アルゴリズムをいろいろ考えてみよう
(修正履歴)
- 2022/03/23 @StrawBerryMoon さん、ご指摘ありがとうございました。
- 2024/02/26 @ringoacidさん、ご指摘ありがとうございました。
1. ひたすら割る
$2$以上の数でひたすら割っていき、割れなければ素数、割れれば合成数。
計算量は$O(N)$
def trial(n: int):
assert n > 1
i = 2
while i < n:
if n%i == 0:
return False
i+=1
return True
2. 偶数は飛ばす
$n\gt2$であれば、素数は全て奇数になる。
1よりも早く判定ができるが、計算量は$O(N)$
def trial_odd(n: int):
assert n > 1
if n == 2:
return True
elif n % 2 == 0:
return False
i = 3
while i < n:
if n%i == 0:
return False
i += 2
return True
3. 探索範囲をsqrt(n)に絞る
実は$\sqrt{n}$まで調べれば、素数かどうかは判別できる($n$の合成数は全て$\sqrt{n}$以下である。)
$n=2^{32}-1$だとしても2と組み合わせてだいたい$32767$回で済める
計算量は$O(N^{1/2})$
def trial_odd_sqrt(n: int):
assert n > 1
if n == 2:
return True
elif n % 2 == 0:
return False
i = 3
while i <= n**0.5:
if n%i == 0:
return False
i+=2
return True
4. Wheel factorizationを使う
$\{2, 3, 5\}$が既知の素数であるとし、$i=7$から探索を開始する。このとき、$\{4, 2, 4, 2, 4, 6, 2, 6\}$の周期で数値を加算しながら素数かどうかを確認できる。
上の数列の場合、$n=7$から始めると、
\begin{align}
&7+4=11\\
&11+2=13\\
&13+4=17\\
&17+2=19\\
&19+4=23\\
&23+6=29\\
&29+2=31\\
&31+6=37\\
&\quad\quad\ldots
\end{align}
となり、見事に2,3,5の倍数が飛ばされている。
これは偶数飛ばしよりもさらに飛ばすことができ、計算量は$O(N/(\log\log N))$1となる。
def trial_wheel_factorize(n: int):
assert n > 1
if n in [2, 3, 5]:
return True
elif n%2 == 0 or n%3 == 0 or n%5 == 0:
return False
i = 7
c = 0
wheel = [4, 2, 4, 2, 4, 6, 2, 6]
while i < n**0.5+1:
if n%i == 0:
return False
i += wheel[c%8]
c += 1
return True
ただし、既知の素数を増やしても使えるというわけではない。周期の数は1つ既知の素数を増やすごとにほぼ指数的に周期が増大する、つまりメモリ的なところで問題が起こってくる。
$p = \{2, 3, 5\} \rightarrow \mathrm{wheel}=\{4, 2, 4, 2, 4, 6, 2, 6\}$
$p = \{2, 3, 5, 7\} \rightarrow \mathrm{wheel}=\{2, 4, 2, 4, 6, 2, 6, 4, 2, 4, 6, 6, 2, 6, 4, 2, 6, 4, 6, 8, 4, 2, 4, 2, 4, 14, 4, 6, 2, 10, 2, 6, 6, 4, 6, 6, 2, 10, 2, 4, 2, 12\}$
$p = \{2, 3, 5, 7, 11\} \rightarrow \mathrm{wheel}=\{4, 2, 4, 6, 2, 6, 4, 2, 4, 6, 6, 2, 6, 4, 2, 6, 4, 6, 8, 4, 2, 4, 2, 4, 14, 4, 6, 2, 10, 2, 6, 6, 4, 2, 4, 6, 2, 10, 2, 4, 2, 12, 10, 2, 4, 2, 4, 6, 2, 6, 4, 6, 6, 6, 2, 6, 4, 2, 6, 4, 6, 8, 4, 2, 4, 6, 8, 6, 10, 2, 4, 6, 2, 6, 6, 4, 2, 4, 6, 2, 6, 4, 2, 6, 10, 2, 10, 2, 4, 2, 4, 6, 8, 4, 2, 4, 12, 2, 6, 4, 2, 6, 4, 6, 12, 2, 4, 2, 4, 8, 6, 4, 6, 2, 4, 6, 2, 6, 10, 2, 4, 6, 2, 6, 4, 2, 4, 2, 10, 2, 10, 2, 4, 6, 6, 2, 6, 6, 4, 6, 6, 2, 6, 4, 2, 6, 4, 6, 8, 4, 2, 6, 4, 8, 6, 4, 6, 2, 4, 6, 8, 6, 4, 2, 10, 2, 6, 4, 2, 4, 2, 10, 2, 10, 2, 4, 2, 4, 8, 6, 4, 2, 4, 6, 6, 2, 6, 4, 8, 4, 6, 8, 4, 2, 4, 2, 4, 8, 6, 4, 6, 6, 6, 2, 6, 6, 4, 2, 4, 6, 2, 6, 4, 2, 4, 2, 10, 2, 10, 2, 6, 4, 6, 2, 6, 4, 2, 4, 6, 6, 8, 4, 2, 6, 10, 8, 4, 2, 4, 2, 4, 8, 10, 6, 2, 4, 8, 6, 6, 4, 2, 4, 6, 2, 6, 4, 6, 2, 10, 2, 10, 2, 4, 2, 4, 6, 2, 6, 4, 2, 4, 6, 6, 2, 6, 6, 6, 4, 6, 8, 4, 2, 4, 2, 4, 8, 6, 4, 8, 4, 6, 2, 6, 6, 4, 2, 4, 6, 8, 4, 2, 4, 2, 10, 2, 10, 2, 4, 2, 4, 6, 2, 10, 2, 4, 6, 8, 6, 4, 2, 6, 4, 6, 8, 4, 6, 2, 4, 8, 6, 4, 6, 2, 4, 6, 2, 6, 6, 4, 6, 6, 2, 6, 6, 4, 2, 10, 2, 10, 2, 4, 2, 4, 6, 2, 6, 4, 2, 10, 6, 2, 6, 4, 2, 6, 4, 6, 8, 4, 2, 4, 2, 12, 6, 4, 6, 2, 4, 6, 2, 12, 4, 2, 4, 8, 6, 4, 2, 4, 2, 10, 2, 10, 6, 2, 4, 6, 2, 6, 4, 2, 4, 6, 6, 2, 6, 4, 2, 10, 6, 8, 6, 4, 2, 4, 8, 6, 4, 6, 2, 4, 6, 2, 6, 6, 6, 4, 6, 2, 6, 4, 2, 4, 2, 10, 12, 2, 4, 2, 10, 2, 6, 4, 2, 4, 6, 6, 2, 10, 2, 6, 4, 14, 4, 2, 4, 2, 4, 8, 6, 4, 6, 2, 4, 6, 2, 6, 6, 4, 2, 4, 6, 2, 6, 4, 2, 4, 12, 2, 12\}$
$p = \{2, 3, 5, 7, 11, 13\} \rightarrow 5760~~\mathrm{cycle}$ (フル)
$p = \{2, 3, 5, 7, 11, 13, 17\} \rightarrow 92160~~\mathrm{cycle}$ (フル)
$p = \{2, 3, 5, 7, 11, 13, 17, 19\} \rightarrow 1658880~~\mathrm{cycle}$ (フル)
$p = \{2, 3, 5, 7, 11, 13, 17, 19, 23\} \rightarrow 36495360~~\mathrm{cycle}$
5. AKS素数判定法を使う
省略。こちらの記事が詳しく解説されています。(決定的アルゴリズム)
計算量は$O((\log(N))^{7.5})$。
5. ミラーラビン素数判定法を使う
確率的アルゴリズム。100%ではないが素数かどうかを判定できる。(kを使わない場合)約$n\lt 3.317\times10^{24}$までなら決定的なアルゴリズムにできることが分かっている。
計算量は$O(k\log^3{(N)})$。決定的になると変わるかもしれない。
kを使う場合
誤差は$4^{-k}$。$k=20$であれば$\sim9.094\times 10^{-13}$。ほぼ正解にはなるけど、念をもって$k=50$とかにしてもあり。
import random
def mr_test_normal(n: int):
assert n > 1
if n == 2:
return True
elif n%2 == 0:
return False
d = n - 1
s = 0
while d & 1 == 0:
s += 1; d >>= 1
k = 20
for _ in range(k):
a = random.randint(2, n-1)
x = pow(a, d, n) # a^d mod n
if x == 1 or x == n-1:
continue
for _ in range(s):
x = pow(x, 2, n)
if x == n-1: # a^(d*2^r) mod n
break
else: return False
return True
もっと誤差を低減したい場合、その手法について研究がなされていた大阪市立大の修士の方がいますので、そちらを参考にしてください(直接リンクは参照できませんが、Googleで「ミラーラビン素数判定法」と調べればヒットします)。
kを使わない場合
def mr_test(n: int):
assert n > 1
if n == 2:
return True
elif n%2 == 0:
return False
assert n < int(3.317 * 10**24)
d = n - 1
s = 0
while d & 1 == 0:
s += 1; d >>= 1
it = (2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41)
for a in it:
if a == n:
return False
a %= n
x = pow(a, d, n) # a^d mod n
if x == 1 or x == n-1:
continue
for _ in range(s):
x = pow(x, 2, n)
if x == n-1: # a^(d*2^r) mod n
break
else: return False
return True
もうすこし厳密
def select_base(n: int) -> tuple[int, ...]:
bases: tuple[int, ...]
if n < 2047:
bases = (2,)
elif n < 1373653:
bases = (2, 3)
elif n < 9080191:
bases = (31, 73)
elif n < 25326001:
bases = (2, 3, 5)
elif n < 3215031751:
bases = (2, 3, 5, 7)
elif n < 350269456337:
bases = (
4230279247111683200, 14694767155120705706, 16641139526367750375
)
elif n < 55245642489451:
bases = (
2, 141889084524735, 1199124725622454117, 11096072698276303650
)
elif n < 7999252175582851:
bases = (
2, 4130806001517, 149795463772692060, 186635894390467037,
3967304179347715805
)
elif n < 585226005592931977:
bases = (
2, 123635709730000, 9233062284813009, 43835965440333360,
761179012939631437, 1263739024124850375
)
elif n < 2**64:
bases = (2, 325, 9375, 28178, 450775, 9780504, 1795265022)
elif n < 3317044064679887385961981:
bases = (2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37)
else:
raise OverflowError
return bases
def mr_test(n: int):
assert n > 1
if n == 2:
return True
elif n % 2 == 0:
return False
d = n - 1
s = 0
while d & 1 == 0:
s += 1
d >>= 1
mr_bases = select_base(n)
for a in mr_bases:
if a == n:
return False
a %= n
x = pow(a, d, n) # a^d mod n
if x == 1 or x == n-1:
continue
for _ in range(s):
x = pow(x, 2, n)
if x == n-1: # a^(d*2^r) mod n
break
else:
return False
return True
# mr_test(15485863) # -> True
(Wikipediaをもとに作成。)
おわりに
こんな感じでいろんなやり方があります。小規模であればwheel factorizationくらいで十分ですが、大きくなってくるとミラーラビンがいかにすごいかがよくわかってきます。