与えられた数 \(n\) が素数かどうかを 誤り率 \(4^{-k}\)、時間 \(O(\log(n) k)\) で判定する、ミラー・ラビン素数判定 (Miller-Rabin primality test) を紹介します（計算量は \(n\) が word size に収まるとした場合）。
フェルマーテスト
\(n\) が素数のとき、フェルマーの小定理より任意の \(x \neq 0\) について \(x^{n-1} = 1 \bmod n\) が成り立ちます。従って、ランダムな \(x \in [1,n)\) について、この等式が必ず成り立つならば \(n\) が素数であると高い確率で言えそうです。また、この等式を満たさない \(x\) が存在すれば、\(n\) は合成数です。この素数判定法をフェルマーテストと言います。
しかし、残念ながら合成数であるにも関わらず、フェルマーテストが高い確率で失敗する数が存在します。\(n\) と互いに素な任意の \(a\) について \(x^{n-1} = 1 \bmod n\) となるような \(n\) が無限に存在し、カーマイケル数といいます。
最小のカーマイケル数は \(n=561=3 \cdot 11 \cdot 17\) です。カーマイケル数であることは、中国剰余定理より各素因数に対して
- \(561=2 \cdot 280 + 1 = 1\bmod 2\)
- \(561=10 \cdot 56 + 1 = 1\bmod 10\)
- \(561=16 \cdot 35 + 1 = 1\bmod 16\)
が成り立つことから従います。
\(n\) と互いに素でない数 \(a\) をフェルマーテストに掛けると、必ず \(n\) が合成数だと判明します。例えば、\(n=p_1p_2p_3\) がカーマイケル数ならば \((1-1/p_1)(1-1/p_2)(1-1/p_3)\) の確率で合成数だと分かりません。従って、素数と判明するまでには \(p_1\) 回程度試行する必要があり \(n\) の多項式時間になってしまいます。そこで、ミラー・ラビンの素数判定法へと改良します。
ミラー・ラビン素数判定法
素数 \(p\) ならば、\(x^2=1 \bmod p\) は高々 2 個しか解を持たず、\(1\) の平方根は \(\pm 1\) のみです。一方 \(n\) が合成数の場合は複数の平方根を持ち得ます。例えば、カーマイケル数 \(n=561\) を mod とする 1 の平方根は \(x=1,67,188,254,307,373,494,560\) の 8 つあります。
\(n-1=q2^Q\) (\(q\) : 奇数) とします。フェルマーテストで合成数と確定しなかったランダムな数 \(a\) （つまり \(a^{n-1}=1\)）に対して、
$$a^{q}, a^{q2^1}, a^{q2^2}\ldots, a^{q2^{Q-1}}$$に \(\pm 1\) 以外の \(1\) の平方根が含まれていないか調べるのがミラー・ラビン素数判定法です。
ミラー・ラビン素数判定法ではランダムな \(a\) に対して、\(1-4^{-1}\) の確率で \(n\) が合成数だと確定させることができます（証明参照）。\(k\) 個のランダムな数に対して行えば、確率 \(1-4^{-k}\) です。計算量は一回当たり \(O(\log(n))\) です。
誤り率の証明
片側誤り率 1/4 以下であることを証明します。
偶数の素数は \(2\) のみですから、\(n=\prod_{i=1}^m p_i^{e_i}\) が奇数と仮定しても一般性を失いません。
\(a \not \in (\mathbb{Z}/n\mathbb{Z})^\times\) は \(\gcd(a,n) \neq 1\) より \(a^{n-1} \neq 1 \bmod n\) です。従って、\(n\) が合成数だと判定できます。
\((\mathbb{Z}/n\mathbb{Z})^\times\) の元を考えましょう。
\(x^{n-1}=1 \pmod n\) を満たす \(x\) からなる部分集合 \(H\) は \((\mathbb{Z}/n\mathbb{Z})^\times\) の部分群をなします。ラグランジュの定理より、\((\mathbb{Z}/n\mathbb{Z})^\times\) は互いに共通部分を持たない \(H,g_1H,\ldots,g_sH\) の和集合$$(\mathbb{Z}/n\mathbb{Z})^\times = H \cup g_1 H \cup \cdots \cup g_s H$$ で表せます。\(a \not \in H\) だと、\(a^{n-1} \neq 1\) より \(n\) が合成数だと判定できます。フェルマーテストの成功率は \(s/(s+1)\) です。
(i) \(n\) が平方因子 \(p^2\) を持つ場合：
\(a_k=1+kp\) （\(k \neq 0\)：原始根）とすると $$\begin{align}a_k^{n-1} &= 1 + k(n-1)p \bmod p^2\\ &=1-kp\bmod p^2\end{align}$$ です。\(k=1,2,\ldots,p-1\) それぞれに対して、\(a_k\) が異なる \(g_iH\) に属するので、\(s \geq p-1\) です。\(p=3\) のとき、テストをすり抜けるのは \(x=1\) のみなので失敗確率 \(1/8\) です。\(p \geq 5\) のとき、\(s \geq 4\) より失敗確率 \(1/5\) 以下です。
(ii) 相異なる素数の積 \(n = pq\) の場合（\(p < q\)）：
\(a^{n-1}=a^{pq-1}=a^{q-1}\bmod p\) ですが、\(x^{q-1} = 1 \bmod p\) は解を \(\gcd(q-1,p-1)\) 個しか持たないので \(s \geq 1\) です。よって失敗確率 \(1/2\) 以下です。\(1/4\) 以下であることを示すには、フェルマーテストが失敗する場合に当たる \(H\) の構造を調べる必要があります。
\(n-1=q2^Q\) (\(q\) : 奇数) とします。\(x^q=1\) を満たす \(x\) からなる部分集合 \(F \subset H\) も部分群をなし、互いに共通部分を持たない \(F,a_1F,\ldots,a_tF\) の和集合 $$ H = F \cup a_1 F \cup \cdots \cup a_t F$$ で表せます。\(q\) が奇数より \(b \in F \Rightarrow -b \not\in F\) にもかかわらず \(-b \in H\) です。中国剰余定理で素冪の mod にバラシて適用することで、\(t \geq 2^m\) が得られます（\(n=\prod_{i=1}^m p_i^{e_i}\)）。つまり \(|H| \geq 2^m|F|\) です。
\(H\) を次の３つに分割します。
- \(F = \{x : x^q=1\}\)
- \(F_1 = \{x : \) \(x^{q2^s}=-1, x^{q2^{s+1}}=1\) なる \(s\) が存在\(\}\)
- \(F_2 = \{x : \) \(x^{q2^s}\neq \pm1, x^{q2^{s+1}}=1\) なる \(s\) が存在\(\}\)
\(F,F_1,F_2\) のうち \(a \in F_2\) の場合だけ \(a\) が合成数だと判定できます。\(F,F_1\) の全体に占める割合を知りたいです。中国剰余定理より、素因数分解 \(n=\prod_{i=1}^m p_i^{e_i}\) の各素冪 \(p_i^{e_i}\) の mod に分解できます。つまり、\(x \in F_1\) ならば $$x^{q2^s}=-1 \bmod n \Leftrightarrow \begin{cases}
x^{q2^s}=-1 \bmod p_1^{e_1} \\
x^{q2^s}=-1 \bmod p_2^{e_2} \\
\vdots \\
x^{q2^s}=-1 \bmod p_m^{e_m} \\
\end{cases}$$ です。各素冪に対する mod を独立に 1 にすることができるので、\(|F_2| \geq (2^m-1)|F_1| \) です。\(2^m|F|\leq |H|\) と併せて
$$\begin{align}
|F|+|F_1| \leq& |F| + (|H|-|F|)/(2^m-1) \\
\leq& |H|/(2^m-1)+|F|(2^m-2)/(2^m-1) \\
=& \frac{2-2^{-(m-1)}}{2^m-1} |H|\\
\end{align}$$を得ます。\(m=2\) つまり (ii) のとき \(|F|+|F_1| \leq 2^{-1}|H|\) ですから、失敗確率は 1/4 以下です。
(iii) \(n\) が 3 つ以上の素因数を持つ場合：
このとき、\(m\geq3\) より \(|F|+|F_1| \leq 4^{-1}|H|\) なので失敗確率は 1/4 以下です。
実装
C++の実装は次の通りです：
#include <iostream>
#include <random>
using namespace std;
long long pow(long long a, long long n, long long p) {
return n!=0?pow(a*a%p,n/2,p)*(n%2==1?a:1)%p:1;
}
// n が素数か判定
bool miller_rabin(long long n) {
if (n % 2 == 0) return n == 2;
long long q = n, Q = 0;
while (q % 2 == 0) {
q/=2;
++Q;
}
random_device seed_gen;
default_random_engine engine(seed_gen());
uniform_int_distribution<> dist(1, n - 1);
for (int i = 0; i < 20; ++i) {
long long a = dist(engine);
if (pow(a, n - 1, n) != 1) {
return false;
}
a = pow(a, q, n);
for (int j = 0; j < Q; ++j) {
long long na = a * a % n;
if (na == 1 && a == n - 1) return false;
a = na;
}
}
return true;
}
int main() {
for (int i=2;i<100;++i)
cout<<i<<" "<<miller_rabin(i)<<endl;
}
Comments