素数判定アルゴリズムいろいろ

与えられた$n$について、素数かどうかを判定したい。
アルゴリズムをいろいろ考えてみよう

(修正履歴)
- 2022/03/23 @StrawBerryMoon さん、ご指摘ありがとうございました。
- 2024/02/26 @ringoacidさん、ご指摘ありがとうございました。

# 1. ひたすら割る
$2$以上の数でひたすら割っていき、割れなければ素数、割れれば合成数。
計算量は$O(N)$

```python
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)$

```python
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})$

```python
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を使う

https://qiita.com/ppza53893/items/49e480d22858a449dc69

$\\{2, 3, 5\\}$が既知の素数であるとし、$i=7$から探索を開始する。このとき、$\\{4, 2, 4, 2, 4, 6, 2, 6\\}$の周期で数値を加算しながら素数かどうかを確認できる。
上の数列の場合、$n=7$から始めると、

```math
\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]となる。

[^1]: https://handwiki.org/wiki/Generating_primes

```python
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}$ ([フル](https://raw.githubusercontent.com/ppza53893/wheel_c/master/skip_2_3_5_7_11_13.txt))
$p = \\{2, 3, 5, 7, 11, 13, 17\\} \rightarrow 92160~~\mathrm{cycle}$ ([フル](https://raw.githubusercontent.com/ppza53893/wheel_c/master/skip_2_3_5_7_11_13_17.txt))
$p = \\{2, 3, 5, 7, 11, 13, 17, 19\\} \rightarrow 1658880~~\mathrm{cycle}$ ([フル](https://github.com/ppza53893/wheel_c/blob/master/skip_2_3_5_7_11_13_17_19.txt?raw=true))
$p = \\{2, 3, 5, 7, 11, 13, 17, 19, 23\\} \rightarrow 36495360~~\mathrm{cycle}$ 


# 5. AKS素数判定法を使う
省略。こちらの記事が詳しく解説されています。(**決定的アルゴリズム**)
計算量は$O((\log(N))^{7.5})$。

https://qiita.com/srtk86/items/aaaf7225d9861512a63d

# 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$とかにしてもあり。

```python
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を使わない場合

```python
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
```

<details><summary>もうすこし厳密</summary><div>

```python
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
```
</div></details>

([Wikipedia](https://en.wikipedia.org/wiki/Miller%E2%80%93Rabin_primality_test#Deterministic_variants)をもとに作成。)

# おわりに
こんな感じでいろんなやり方があります。小規模であればwheel factorizationくらいで十分ですが、大きくなってくるとミラーラビンがいかにすごいかがよくわかってきます。