ミラー・ラビンの素数判定法 (Miller-Rabin 法)

ミラー・ラビンの素数判定法は、その背景にある整数論的考察もめっちゃ面白いので、ぜひそれも味わいましょう！！

1. はじめに

本記事では正の整数 が与えられたときに、 が素数であるか否かを判定する問題を考えます。

よく知られた方法は、 を で順に試し割りしていく方法です。 まで試せばよいことから の計算量となります。競プロでは、 程度の制約下で、1 回だけ素数判定するのが限界というイメージがありますよね。

ところが、本記事で紹介する Miller-Rabin の素数判定法を用いると、2 秒以内の実行制限時間内に、

• 程度の制約下で、
• 回程度の回数の

素数判定が行えます！

例によって、Yosupo Judge Library Checker の問題「Primality Test」を通します。

なお注意点として、Miller-Rabin 法は確率的アルゴリズムであり、通常は、非常に低い確率で「合成数を素数だと誤判定してしまう」ことがあり得ます。しかし、うまくやることで、 以下の範囲では確率 100% で素数判定ができることが示されています。実用的には決定的アルゴリズムといってよいでしょう。

目次

• 1. はじめに
• 2. 素数の性質
• 3. Miller-Rabin 法の原理
  • 数値例
• 4. Miller-Rabin 法
  • 注意：N が小さいときは決定的アルゴリズムに
• 5. Miller-Rabin 法の実装
  • 方法 1：__uint128_t 型を用いる方法
  • 方法 2：モンゴメリ乗算を用いる方法
• 6. おわりに

2. 素数の性質

Miller-Rabin 法を理解する上で欠かせない、素数の重要な性質を 1 つ挙げます。

---

【素数の性質 1】

を整数とし、 を素数とする。

が の倍数であるならば、 のうちの少なくとも 1 つは の倍数である。

---

一見すると当然のようにも思える定理ですが、意外と自明ではなく、奥が深い定理です。有名な「素因数分解の一意性」にも深く関わってくる定理です。ここでは証明はしないですが、関心のある方は高木貞治先生の『初等整数論講義』などを読んでみてください。

もう一つ、Miller-Rabin 法を理解する上で欠かせない素数の定理として、Fermat の小定理があります。

---

【素数の性質 2：Fermat の小定理】

が素数とし、 を で割り切れない整数とする。このとき、

が成立する。

---

Fermat の小定理については、次の記事に詳しく書きました！

qiita.com

Fermat の小定理も、その証明方法も含めて理解することで、整数と仲良くなれる。そんな定理です。

3. Miller-Rabin 法の原理

Miller-Rabin 法の原理となる整数論的性質を深掘りしていきます。ここまでの素数の性質が頭にしっかり染み込んでいれば、Miller-Rabin 法の原理はいたって単純です！

たとえば、 としましょう。このとき、Fermat の小定理より、 を で割り切れない整数として

が成り立ちます。1 を右辺に移行して因数分解してみましょう。

となります。これは が 41 の倍数であることを意味しているので、素数の性質 1 より、これら 4 つの値のうちの少なくとも 1 つが 41 の倍数であることになります。つまり、

• 
• 
• 
• 

のうちのいずれかが成り立つことになります。もしすべて成り立たないとするならば、41 は素数ではないと言い切れることになります。

一般の整数 の場合

続いて、素数かどうかを判定したい一般の整数 に対して、同様のことを考えてみましょう。

が素数であるならば、Fermat の小定理より、 の倍数でない整数 に対して

が成り立ちます。先ほどと同様に、1 を右辺に移して因数分解します。 ( は奇数) と表すことにすると、

となります。したがって、 が素数であるならば、

• 
• 
• 
• 
• 
• 

のうちのいずれかが成り立つことになります。

逆に、これらがどれも成り立たないとき、 は合成数であると言いきれます。ただし、上記の式のうちのいずれかが成り立つ場合に、 が素数だと言い切れるわけではないことには注意しましょう。

ここで、以下の定理が知られています。証明は 37zigen さんの記事にあります。つまり、 をランダムに選ぶのを 回実施してテストすれば、合成数を誤って素数と判定してしまう確率を 以下にできると言えます。

---

が奇数の合成数とし、 を 以上 以下の整数の中からランダムに選ぶとする。

このとき、75% 以上の確率で

• 
• 
• 
• 
• 
• 

が成立する。

---

数値例

・ (素数)

() のとき、 とすると、

となります。 となっています。

・ (合成数)

() のとき、 とすると、

となります。この時点で は合成数だと断定できます。

・ (合成数)

() のとき、 とすると、

となります。これだけでは は素数とも合成数とも言い切れません。しかし とすると、

となります。これによって、 は合成数であると言い切れます。

4. Miller-Rabin 法

Miller-Rabin の素数判定法をまとめます。

---

• のとき：素数である
• が 以外の偶数であるとき：合成数である

以下、 を 以上の奇数とする。 ( は奇数) と表す。

以上 以下の整数 に対して、

• 
• 
• 
• 
• 
• 

のいずれかが成り立つことを「 は のテストに通る」と呼ぶことにする。 のテストに通らない が検出されたならば、直ちに は合成数であると言える。

ランダムに を選択することを 回繰り返して、いずれの もテストに通ることが確かめられた場合、確率 で は素数であると言える。

---

Miller-Rabin 法は、素数を誤って合成数と判定してしまうことは絶対にないですが、合成数を誤って素数と判定してしまうことはあり得ます。しかし、試行回数 を増やせば、誤判定確率を減らせます。

計算量は、

• での整数の乗算 1 回：計算量
• 乗算の回数：

ということで、全体で と評価できます。しかし、競プロ的には、 の範囲における での整数の乗算は定数時間でできると考える方が自然で、実質的な計算量は となります。

注意： が小さいときは決定的アルゴリズムに

このように Miller-Rabin 法は確率的アルゴリズムですが、十分小さい に対しては、テストに用いる を上手に選ぶことで決定的アルゴリズムとなります。

---

• のとき：
• のとき：

に対してテストに通るならば、 は素数であると言い切れる

---

なお、効率的な の選び方について、次のサイトによくまとまっています。

miller-rabin.appspot.com

5. Miller-Rabin 法の実装

Miller-Rabin 法はとてもシンプルなので、速度を気にしなければ実装は容易です。

しかし、「 を で割った余り」を求める部分をどのように実装するかによって、結構な速度差が生じます。Yosupo Judge の Primality Test をベンチマークとしてみましょう。

「 を で割った余り」を計算するときには、どうしても 64 bit 整数の範囲におさまらない値が出てきます。その対策として、次の 2 つの方法を実装してみました。

• 方法 1：単純に __uint128_t 型を用いる方法 (僕の実装で 1783ms)
• 方法 2：モンゴメリ乗算を用いる方法 (僕の実装で 247ms)

方法 1：__uint128_t 型を用いる方法

を計算するとき、 の値は 64 bit 整数におさまりません。そこで、__uint128_t 型を用います。次のコードは、Yosupo Judge で 1783ms でした。

#include <bits/stdc++.h>
using namespace std;

// A^N mod M
template<class T> T pow_mod(T A, T N, T M) {
    T res = 1 % M;
    A %= M;
    while (N) {
        if (N & 1) res = (res * A) % M;
        A = (A * A) % M;
        N >>= 1;
    }
    return res;
}

bool MillerRabin(long long N, vector<long long> A) {
    long long s = 0, d = N - 1;
    while (d % 2 == 0) {
        ++s;
        d >>= 1;
    }
    for (auto a : A) {
        if (N <= a) return true;
        long long t, x = pow_mod<__int128_t>(a, d, N);
        if (x != 1) {
            for (t = 0; t < s; ++t) {
                if (x == N - 1) break;
                x = __int128_t(x) * x % N;
            }
            if (t == s) return false;
        }
    }
    return true;
}

bool is_prime(long long N) {
    if (N <= 1) return false;
    if (N == 2) return true;
    if (N % 2 == 0) return false;
    if (N < 4759123141LL)
        return MillerRabin(N, {2, 7, 61});
    else
        return MillerRabin(N, {2, 325, 9375, 28178, 450775, 9780504, 1795265022});
}

void YosupoPrimarityTest() {
    int N;
    cin >> N;
    for (int i = 0; i < N; ++i) {
        long long A;
        cin >> A;
        if (is_prime(A))
            cout << "Yes" << endl;
        else
            cout << "No" << endl;
    }
}

int main() {
    YosupoPrimarityTest();
}

方法 2：モンゴメリ乗算を用いる方法

モンゴメリ乗算は、時間のかかる除算を実質的に行うことなく、乗算・加算・減算・シフト演算のみで、効率的に整数の積の剰余を求めることのできるアルゴリズムです。

ここでは詳細は省略しますが、モンゴメリ乗算を用いると、32 bit 整数に収まらない整数を法とする効率的な modint を作れます。

次のコードは、Yosupo Judge で 247ms でした。特別なチューニングを施していない単純な実装としては、まずまずの速度だと思います。

#include <bits/stdc++.h>
using namespace std;

// montgomery modint (MOD < 2^62, MOD is odd)
struct MontgomeryModInt64 {
    using mint = MontgomeryModInt64;
    using u64 = uint64_t;
    using u128 = __uint128_t;
    
    // static menber
    static u64 MOD;
    static u64 INV_MOD;  // INV_MOD * MOD ≡ 1 (mod 2^64)
    static u64 T128;  // 2^128 (mod MOD)
    
    // inner value
    u64 val;
    
    // constructor
    MontgomeryModInt64() : val(0) { }
    MontgomeryModInt64(long long v) : val(reduce((u128(v) + MOD) * T128)) { }
    u64 get() const {
        u64 res = reduce(val);
        return res >= MOD ? res - MOD : res;
    }
    
    // mod getter and setter
    static u64 get_mod() { return MOD; }
    static void set_mod(u64 mod) {
        assert(mod < (1LL << 62));
        assert((mod & 1));
        MOD = mod;
        T128 = -u128(mod) % mod;
        INV_MOD = get_inv_mod();
    }
    static u64 get_inv_mod() {
        u64 res = MOD;
        for (int i = 0; i < 5; ++i) res *= 2 - MOD * res;
        return res;
    }
    static u64 reduce(const u128 &v) {
        return (v + u128(u64(v) * u64(-INV_MOD)) * MOD) >> 64;
    }
    
    // arithmetic operators
    mint operator - () const { return mint() - mint(*this); }
    mint operator + (const mint &r) const { return mint(*this) += r; }
    mint operator - (const mint &r) const { return mint(*this) -= r; }
    mint operator * (const mint &r) const { return mint(*this) *= r; }
    mint operator / (const mint &r) const { return mint(*this) /= r; }
    mint& operator += (const mint &r) {
        if ((val += r.val) >= 2 * MOD) val -= 2 * MOD;
        return *this;
    }
    mint& operator -= (const mint &r) {
        if ((val += 2 * MOD - r.val) >= 2 * MOD) val -= 2 * MOD;
        return *this;
    }
    mint& operator *= (const mint &r) {
        val = reduce(u128(val) * r.val);
        return *this;
    }
    mint& operator /= (const mint &r) {
        *this *= r.inv();
        return *this;
    }
    mint inv() const { return pow(MOD - 2); }
    mint pow(u128 n) const {
        mint res(1), mul(*this);
        while (n > 0) {
            if (n & 1) res *= mul;
            mul *= mul;
            n >>= 1;
        }
        return res;
    }

    // other operators
    bool operator == (const mint &r) const {
        return (val >= MOD ? val - MOD : val) == (r.val >= MOD ? r.val - MOD : r.val);
    }
    bool operator != (const mint &r) const {
        return (val >= MOD ? val - MOD : val) != (r.val >= MOD ? r.val - MOD : r.val);
    }
    friend istream& operator >> (istream &is, mint &x) {
        long long t;
        is >> t;
        x = mint(t);
        return is;
    }
    friend ostream& operator << (ostream &os, const mint &x) {
        return os << x.get();
    }
    friend mint modpow(const mint &r, long long n) {
        return r.pow(n);
    }
    friend mint modinv(const mint &r) {
        return r.inv();
    }
};

typename MontgomeryModInt64::u64
MontgomeryModInt64::MOD, MontgomeryModInt64::INV_MOD, MontgomeryModInt64::T128;

// Miller-Rabin
bool MillerRabin(long long N, vector<long long> A) {
    using mint = MontgomeryModInt64;
    mint::set_mod(N);
    
    long long s = 0, d = N - 1;
    while (d % 2 == 0) {
        ++s;
        d >>= 1;
    }
    for (auto a : A) {
        if (N <= a) return true;
        mint x = mint(a).pow(d);
        if (x != 1) {
            long long t;
            for (t = 0; t < s; ++t) {
                if (x == N - 1) break;
                x *= x;
            }
            if (t == s) return false;
        }
    }
    return true;
}

bool is_prime(long long N) {
    if (N <= 1) return false;
    else if (N == 2) return true;
    else if (N % 2 == 0) return false;
    else if (N < 4759123141LL)
        return MillerRabin(N, {2, 7, 61});
    else
        return MillerRabin(N, {2, 325, 9375, 28178, 450775, 9780504, 1795265022});
}

void YosupoPrimarityTest() {
    cin.tie(0);
    ios::sync_with_stdio(false);
    int N;
    cin >> N;
    for (int i = 0; i < N; ++i) {
        long long A;
        cin >> A;
        if (is_prime(A))
            cout << "Yes" << endl;
        else
            cout << "No" << endl;
    }
}

int main() {
    YosupoPrimarityTest();
}

6. おわりに

ここまで、Miller-Rabin 法の原理と実装を見てきました。Miller-Rabin 法の原理は、それ自体が整数論的にとても面白いと思います。みんなもぜひ素敵な整数論ライフを！