💨

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

3 min read

素数判定

前提知識

フェルマーの小定理

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+1pの倍数のときa+1 \equiv 0 \pmod pつまりa \equiv -1 \pmod p

a-1pの倍数のとき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 
}