動機

現代のコンピューターで $2^{64} = 18446744073709551616$ を超える様な巨大な数に対して素数かどうかを判定する際に試し割りやエラトステネスの篩を実行する事は現実的に不可能なので Miller-Rabin の確率的素数判定の実装を紹介する.
多倍精度整数の基本的な実装は以下より.
多倍精度整数の基礎 ( 四則演算編 )
多倍精度整数の基礎 ( Karatsuba乗算編 )
多倍精度整数の基礎 ( Karatsuba除算編 )

実装

// Miller Rabin's Primality Test.
static bool miller_rabins_test(const integer &n, std::size_t s){
    static std::mt19937 g(0xABABABABu);
    integer a, m = n - 2;
    for(std::size_t i = 0; i < s; ++i){
        if(witness(m.random(g) + 1, n)){
            return false;
        }
    }
    return true;
}

static bool witness(const integer &a, const integer &n){
    integer b = n - 1, d = 1, x, t;
    for(int i = static_cast<int>(b.bit_num()) - 1; i >= 0; --i){
        x = d;
        d = (d * d) % n;
        if(d == 1 && x != 1 && x != b){
            return true;
        }
        t = b >> static_cast<UInt>(i);
        if(t > 0 && t.data[0] % 2 == 1){
            d = (d * a) % n;
        }
    }
    return d != 1;
}

// Random Number
template<class Generator>
integer random(Generator &g) const{
    std::size_t s = g() % data.size() + 1;
    integer x;
    x.data.resize(s);
    for(std::size_t i = 0; i < s - 1; ++i){
        x.data[i] = g() & base_mask;
    }
    if(s == data.size()){
        x.data.back() = g() % data.back();
    }
    x.normalize_data_size();
    return x;
}

実行してみる

実装した miller_rabins_test 関数の第一引数には判定したい数を, 第二引数には何度判定を繰り返すかを指定する.

int main(){
    multi_precision::integer<> a("1111111111111111111111"), b("11111111111111111111111");
    std::cout << multi_precision::integer<>::miller_rabins_test(a, 64) << std::endl;
    std::cout << multi_precision::integer<>::miller_rabins_test(b, 64) << std::endl;

    return 0;
}

is_prime[1111111111111111111111] = false
is_prime[11111111111111111111111] = true

全体像

marionette-of-u / integer.hpp

参考

アルゴリズムイントロダクション第 1 版, 33.8節.