Miller–Rabin 素数判定法

前提知識

$p$ を奇素数とする。フェルマーの小定理より、 $p$ と互いに素な正整数 $a$ について

a^{p-1} \equiv 1 \pmod p \qquad (1)

である。

$p-1$ を $2$ で割れるだけ割って $2^{s} d$ と表す。ただし $s$ は正整数、 $d$ は奇数である。上の式は

a^{2^s d} \equiv 1 \pmod p \qquad (2)

と表せる。

このとき、次の定理が成り立つ。

定理①

各変数を上のように定義したとき、

\exists r\in \mathbb{Z}(0 \le r \le s-1) \ \ a^{2^r d} \equiv {-1} \pmod p \qquad \mathrm{(a)}

または

a^d \equiv 1 \pmod p \qquad \mathrm{(b)}

が成り立つ。

以下これを示す。

証明

$a^{2^{s-1} d} \pmod p$ の値について考える。

$a^{2^s d}$ の平方根なので、 $(2)$ よりこれは $1$ または $-1$ と合同である。

$-1$ と合同なら $\mathrm{(a)}$ が成立。

$1$ と合同なら平方根 $a^{2^{s-2} d}$ を再びとり、同じように考える。

$0 \le r \le s-1$ の間で $-1$ と合同なものがあれば $\mathrm{(a)}$ が成立する。

そういうものがない場合、 $a^{2^0 d} \equiv a^d \equiv 1 \pmod p$ となるので $\mathrm{(b)}$ が成立。

証明終わり

これの対偶を取ると、

定理②

$n$ を正整数、 $a$ を $n$ と互いに素な正整数として、

$n-1 = 2^{s} d$ 、ただし $d$ は正奇数、 $s$ は正整数としたとき、

\forall r \in \mathbb{Z}(0\le r\le s-1)\ \ a^{2^r d} \not\equiv {-1} \pmod n \qquad \mathrm{(c)}

かつ

a^d \not\equiv 1 \pmod n \quad \mathrm{(d)}

ならば、 $n$ は合成数

がいえる。

判定法

定理②を用いて確率的な素数判定をするものが Miller–Rabin 素数判定法である。

具体的には、 n を素数判定したい正整数とすると、

$n-1$ を $2$ の冪 $2^s$ と奇数 $d$ の積に分解する
以下の試行を何回か繰り返す
1. 正整数 $0 < a \le n-1$ をランダムに選ぶ
2. $a$ と $n$ が互いに素でないなら $n$ は合成数。終了
3. 上の条件 $\mathrm{(c)}, \mathrm{(d)}$ をともに満たす場合、 $n$ は合成数。終了
4. i. に戻る
試行に引っかからなかった場合、 $n$ は確率的素数である。終了

ここで、定理②の逆は成り立たない。つまり $n$ が合成数でも上の2つの条件 $\mathrm{(c)}, \mathrm{(d)}$ を満たさない場合がある。

しかし $a$ を様々に変えて試行すれば、いつかは $\mathrm{(c)}, \mathrm{(d)}$ を満たす $a$ が見つかる可能性が上がる。つまり $n$ を合成数と判定できる確率が上がる。

条件付き決定的判定

また、特定の底 $a$ について判定を行えば、ある範囲では素数であると決定的に判定できるような $a$ と範囲がいくつか発見されている。

例えば、 $2^{64}$ 未満の数は、 $a=2,325,9375,28178,450775,9780504,1795265022$ の7つの底について判定を行えば、素数であると決定的に判定できることが分かっている。

実装例

TypeScript


export const millerRabin = (n: bigint) => {
	if (n < 0n) throw Error('引数は正の整数でなければなりません');

	// 2, 2の倍数 の場合の判定
	if (n === 2n) return true;
	if (n === 1n || n % 2n === 0n) return false;

	const bit_num = n.toString(2).length;


	const s = BigInt((n - 1n).toString(2).match(/0+$/g)?.[0].length ?? 0);
	const d = (n - 1n) >> s;

	if (n < 2n ** 64n) {
		/** n が 2^64 未満の時、決定的に判定できる 参考: https://miller-rabin.appspot.com/#bases7 */
		const bases_under_64 = [2n, 325n, 9375n, 28178n, 450775n, 9780504n, 1795265022n] as const;

		challenge: for (const b_ of bases_under_64) {
			const base = b_ >= n ? b_ % n : b_;

			if (base === 0n) continue challenge;

			if (exEuclidean(base, n).gcd != 1n) return false;

			let y = modPow(base, d, n);

			if (y === 1n) continue challenge;

			for (let i = 0n; i < s; i++) {
				if (y === n - 1n) continue challenge;

				y = (y * y) % n;
			}
			return false;
		}
		return true;
	} else {
		/** 試行回数 */
		const max_rot = 40;

		challenge2: for (let i = 0; i < max_rot; i++) {
			let b_ = 0n;

			while (b_ < 2n || b_ >= n) {
				b_ = getRandBIByBitLength(bit_num); // bit_num ビット以下の乱数を出力
			}

			if (exEuclidean(b_, n).gcd !== 1n) return false;

			const base = b_;

			let y = modPow(base, d, n);

			if (y === 1n) continue challenge2;

			for (let i = 0n; i < s; i++) {
				if (y === n - 1n) continue challenge2;
				y = (y * y) % n;
			}
			return false;
		}
		return true;
	}
};