Baillie-PSW素数判定テストをC++で実装する方法

Baillie-PSW素数判定テストは、Robert Baillie、Carl Pomerance、John Selfridge、Samuel Wagstaffという4人の数学者にちなんで名付けられた素数判定法です。このテストは、与えられた数が合成数であるか、それともおそらく素数であるかを判定します。

Baillie-PSWテストの特徴

Baillie-PSWテストは、ミラー・ラビン素数判定法と強いルーカス疑似素数判定を組み合わせた手法として知られています。現在までに、このテストが誤って素数と判定してしまう合成数（反例）は一つも発見されていません。そのため、暗号鍵の生成など、大きな数の素数性を確認したい場面で活用される、信頼性の高い確率的素数判定法です。
ここでは、その中核となるミラー・ラビン型の判定処理をC++で実装した例を紹介します。

アルゴリズム

MillerTest() の流れ

開始
　　ブール型の関数 MillerTest を宣言する。
　　整数型の MT_dt と MT_num を宣言し、引数として受け取る。
　　整数型の MT_a と MT_x を宣言する。
　　　　MT_a = 2 + rand( ) % (MT_num - 4) で初期化する。
　　　　MT_x = pow(MT_a, MT_dt, MT_num) で初期化する。
　　if (MT_x == 1 || MT_x == MT_num - 1) ならば
　　　　true を返す。
　　while (MT_dt != MT_num - 1) の間、繰り返す
　　　　MT_x = (MT_x * MT_x) % MT_num を計算する。
　　　　MT_dt *= 2 を実行する。
　　　　if (MT_x == 1) ならば
　　　　　　false を返す。
　　　　if (MT_x == MT_num - 1) ならば
　　　　　　true を返す。
　　false を返す。
終了

C++での実装例

#include <iostream>
#include<stdlib.h>
using namespace std;
// 繰り返し二乗法により冪乗の剰余を高速に計算する関数
int pow(int pow_a, unsigned int pow_b, int pow_c) {
    int result = 1;
    pow_a = pow_a % pow_c;
    while (pow_b > 0) {
        if (pow_b & 1)
        result = (result * pow_a) % pow_c;
        pow_b = pow_b >> 1;
        pow_a = (pow_a * pow_a) % pow_c;
    }
    return result;
}
// ミラー・ラビン判定の1回分のテスト
bool MiillerTest(int MT_dt, int MT_num) {
    // 2 ～ MT_num - 2 の範囲からランダムな底を選ぶ
    int MT_a = 2 + rand( ) % (MT_num - 4);
    int MT_x = pow(MT_a, MT_dt, MT_num);
    if (MT_x == 1 || MT_x == MT_num - 1)
        return true;
    while (MT_dt != MT_num - 1) {
        MT_x = (MT_x * MT_x) % MT_num;
        MT_dt *= 2;
        if (MT_x == 1)
            return false;
        if (MT_x == MT_num - 1)
            return true;
    }
    return false;
}
// 素数判定を行う関数（P_K 回テストを繰り返す）
bool prime(int P_N, int P_K) {
    if (P_N <= 1 || P_N == 4)
        return false;
    if (P_N <= 3)
        return true;
    int P_D = P_N - 1;
    while (P_D % 2 == 0)
        P_D /= 2;
    for (int i = 0; i < P_K; i++)
        if (MiillerTest(P_D, P_N) == false)
            return false;
        return true;
}
int main() {
    int iter = 50; // テストの繰り返し回数
    long num1;
    long num2;
    cout<< "Enter the first number: ";
    cin>>num1;
    cout<<endl;
    if (prime(num1, iter))
        cout<<num1<<" is a prime number
"<<endl;
    else
        cout<<num1<<" is a composite number
"<<endl;
    cout<<"Enter another number: ";
    cin>>num2;
    cout<<endl;
    if (prime(num2, iter))
        cout<<num2<<" is a prime number
"<<endl;
    else
        cout<<num2<<" is a composite number
"<<endl;
    return 0;
}

実行結果

Enter the first number: 23
23 is a prime number
Enter another number: 45
45 is a composite number

プログラムのポイント

• pow() 関数：繰り返し二乗法を用いて a^d mod n を効率的に計算します。通常の冪乗計算と比べて大幅に高速です。
• MiillerTest() 関数：ランダムに選んだ底 a を使って、n が合成数である証拠を探します。条件を満たさなければ、その時点で「合成数」と断定できます。
• prime() 関数：n − 1 を 2 で割れるだけ割って d を求め、その d を使ってテストを P_K 回（ここでは50回）繰り返します。すべてのテストを通過すれば「おそらく素数」と判定します。
• 反復回数 iter = 50：テストの回数を増やすほど誤判定の確率は指数的に低下します。50回程度繰り返せば、実用上ほぼ確実な判定が可能です。