ミラーラビン素数判定法についてメモ

素数判定

前提知識

フェルマーの小定理

$$pを素数とすると \\ a^{p-1} \equiv 1 \pmod p \\$$

証明略

もう一つ

$$pを素数とすると \\ a^2 \equiv 1 \pmod p \\ が成り立つaは1と－1 \pmod pしかない$$

証明

式を変形すると

$$a^2 - 1 \equiv 0 \pmod p \\ (a + 1)(a - 1) \equiv 0 \pmod p$$

$(a+1)(a-1)$がpの倍数になる

pは素数のため2つの数の積で表されるとしたら$1 \times p$の場合のみである

a = 0の場合に$a+1$が1になるが、$a^{2} \equiv 1$を満たさないため不適

よって$a+1$または$a-1$のいずれかがpの倍数になる

$a+1$が$p$の倍数のとき$a+1 \equiv 0 \pmod p$つまり$a \equiv -1 \pmod p$

$a-1$が$p$の倍数のとき$a+1 \equiv 0 \pmod p$つまり$a \equiv 1 \pmod p$

フェルマーテスト

$n$が素数がどうか判定したい

適当に$a$(nと互いに素)をとってきて$a^{p-1} \equiv 1 \pmod n$が成り立つが調べる

成り立っていなければ確実に素数でないといえる

なぜならフェルマーの小定理の対偶をとると、

$$a^{p-1} \not \equiv 1 \pmod p \\ ならばpは素数でない$$

となるため

ただし、成り立っていても確実に素数かどうかはわからない

たくさん$a$をとってきてすべて成り立っていればよしとする

うまくいかないaを取ってきてしまう確率は基本的に$1/2$くらい

フェルマーテストの問題点

いくらaをたくさんとってきても

全部$a^{p-1} \equiv 1 \pmod n$が成り立ってしまうようなnが無数にある(カーマイケル数と呼ばれる)

561, 1105など

ミラーラビンテスト

pを奇数の素数とする

するとp-1は偶数なので

$$p-1 = 2^{s}t$$

と置くことができる

このときに以下が成り立つ

$$pが素数のとき \\ a^t \equiv 1 \pmod p \\ またはあるr(0 \le r \le s-1)について \\ a^{2^{r}t} \equiv -1 \pmod p \\ が成り立つ$$

証明

フェルマーの小定理より

$$a^{2^{s}t} \equiv 1 \pmod p$$

この平方根つまり$a^{2^{s-1}t}$は1または-1(mod p)になる

-1の場合は、(7)の後者の式が成り立つ

1の場合はさらに平方根を考えていく。

そうして考えていったときにr=0になるまでひとつも-1になるものが存在しなかったとする。(全部1)

このとき

$$a^{2^{0}t} = a^{t} \equiv 1 \pmod n$$

これは(5)の前者の式を満たす

ミラーラビンテストは(5)の対偶を使う

$$a^t \not \equiv 1 \pmod p \\ かつすべてのr(0 \le r \le s-1)について \\ a^{2^{r}t} \not \equiv -1 \pmod p \\ ならばpは合成数$$

これもこの条件を満たしていれば確実に合成数だと言えるが、

満たしていなければ確実に素数だといえるわけではない

これもたくさんaを取ってきてすべて満たしていないか調べる

うまくいかないaを選んでしまう確率は最大で$1/4$

決定的アルゴリズム

求めたい数nが$n \le 4,759,123,141$ならば

aは${2, 7, 61}$だけ調べればいい

$n \le 2^{64}$ならば

$2, 325, 9375, 28178, 450775, 9780504, 1795265022$だけ調べればいい

参考

http://miller-rabin.appspot.com/

実装

$1 \le n \le 2^{64}$で正しく動作する

fn modpow(base: i64, mut exp: i64, m: i64) -> i64 {
    let mut ret = 1;
    let mut pow = base;
    while exp != 0 {
        if exp & 1 != 0 {
            ret = ((ret as i128 * pow as i128) % m as i128) as i64;
        }
        pow = ((pow as i128 * pow as i128) % m as i128) as i64;
        exp >>= 1;
    }
    ret
}

fn rabin_miller(n: i64) -> bool {
    if n <= 2 {
        return n == 2
    }
    if n % 2 == 0 {
        return false;
    }
    let mut s = 0;
    let mut t = n-1;
    while t % 2 == 0 {
        s += 1;
        t /= 2;
    }
    let is_prime = |a| {
        let mut x = modpow(a, t, n);
        if x == 1 || x == n-1 {
            return true;
        }
        for _ in 0..s-1 {
            x = ((x as i128 * x as i128) % n as i128) as i64;
            if x == n-1 {
                return true;
            }
        }
        false
    };
    for &a in &[2, 325, 9375, 28178, 450775, 9780504, 1795265022] {
        if a >= n {
            break;
        }
        if !is_prime(a) {
            return false;
        }
    }
    true 
}