ソロベイ・シュトラッセン素数性テストをC++で実装して、指定した数が素数かどうかを判定する方法

ソロベイ・シュトラッassen素数性テストとは

ソロベイ・シュトラッassen（Solovay–Strassen）素数性テストは、ある数が合成数か、それともおそらく素数であるかを判定するための確率的アルゴリズムです。オイラーの規準とヤコビ記号を利用しており、ランダムな試行を複数回繰り返すことで、誤判定の確率を極めて低く抑えられるのが特徴です。

オイラーの規準によれば、p が素数である場合、任意の整数 a に対して次の関係が成り立ちます。

a(p−1)/2 ≡ (a/p) (mod p)

ここで (a/p) はヤコビ記号です。この等式が成立しない数を見つけられれば、その数は確実に合成数だと判断できます。

プログラムの構成要素

本プログラムは、次の3つの関数で構成されています。

1. べき乗剰余を高速に計算する modulo 関数

二進法（バイナリ計算）を用いて a^b mod m を効率的に求めます。指数を半分ずつ減らしながら処理するため、非常に大きな指数でも高速に計算できます。

2. ヤコビ記号を計算する Jacobian 関数

与えられた整数 a と奇数 n に対するヤコビ記号 (a/n) を反復的に計算します。負数の処理、2で割れる間の処理、ユークリッドの互除法に似た入れ替え操作などを組み合わせて求めます。

3. 判定本体の Solovoystrassen 関数

指定された回数だけ以下の手順を繰り返します。

• p < 2 の場合は合成数として扱う
• p が 2 以外の偶数の場合も合成数として扱う
• 1 以上 p−1 以下の乱数 a を生成する
• jacob = (p + Jacobian(a, p)) % p を計算する
• mod = a^(p−1)/2 mod p を計算する
• jacob が 0、または mod ≠ jacob であれば「合成数」と判定して終了する

すべての試行に合格すれば、「おそらく素数」であると判定します。各試行で合成数が素数と誤認される確率は最大でも 1/2 なので、50回回繰り返せば誤判定の確率は (1/2)⁵⁰ ≒ 8.9×10⁻¹⁶ 以下になります。

C++サンプルコード

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

// 二進法によるべき乗剰余計算
long modulo(long m_base, long m_exp, long m_mod) {
    long a = 1;
    long b = m_base;
    while (m_exp > 0) {
        if (m_exp % 2 == 1)
            a = (a * b) % m_mod;
        b = (b * b) % m_mod;
        m_exp = m_exp / 2;
    }
    return a % m_mod;
}

// ヤコビ記号 (a/n) を計算
int Jacobian(long CJ_a, long CJ_n) {
    if (!CJ_a)
        return 0; // (0/n) = 0
    int answer = 1;
    if (CJ_a < 0) {
        CJ_a = -CJ_a; // (a/n) = (-a/n)*(-1/n)
        if (CJ_n % 4 == 3)
            answer = -answer; // n ≡ 3 (mod 4) のとき (-1/n) = -1
    }
    if (CJ_a == 1)
        return answer; // (1/n) = 1
    while (CJ_a) {
        if (CJ_a < 0) {
            CJ_a = -CJ_a;
            if (CJ_n % 4 == 3)
                answer = -answer;
        }
        while (CJ_a % 2 == 0) {
            CJ_a = CJ_a / 2;
            if (CJ_n % 8 == 3 || CJ_n % 8 == 5)
                answer = -answer;
        }
        swap(CJ_a, CJ_n);
        if (CJ_a % 4 == 3 && CJ_n % 4 == 3)
            answer = -answer;
        CJ_a = CJ_a % CJ_n;
        if (CJ_a > CJ_n / 2)
            CJ_a = CJ_a - CJ_n;
    }
    if (CJ_n == 1)
        return answer;
    return 0;
}

// ソロベイ・シュトラッassen素数性テストを実行
bool Solovoystrassen(long SS_p, int itr) {
    if (SS_p < 2)
        return false;
    if (SS_p != 2 && SS_p % 2 == 0)
        return false;
    for (int i = 0; i < itr; i++) {
        // ランダムな数 a を生成
        long a = rand() % (SS_p - 1) + 1;
        long jacob = (SS_p + Jacobian(a, SS_p)) % SS_p;
        long mod = modulo(a, (SS_p - 1) / 2, SS_p);
        if (!jacob || mod != jacob)
            return false;
    }
    return true;
}

// メイン関数
int main() {
    int iter = 50; // 試行回数
    long num1;
    long num2;
    cout<< "最初の数値を入力してください: ";
    cin>>num1;
    cout<<endl;
    if (Solovoystrassen(num1, iter))
        cout<<num1<<" は素数です\n"<<endl;
    else
        cout<<num1<<" は合成数です\n"<<endl;
    cout<<"別の数値を入力してください: ";
    cin>>num2;
    cout<<endl;
    if (Solovoystrassen(num2, iter))
        cout<<num2<<" は素数です\n"<<endl;
    else
        cout<<num2<<" は合成数です\n"<<endl;
    return 0;
}

実行結果

最初の数値を入力してください: 24
24 は合成数です

別の数値を入力してください: 23
23 は素数です

まとめ

このように、ソロベイ・シュトラッassen素数性テストを使えば、ヤコビ記号とべき乗剰余の計算を組み合わせるだけで、大きな数に対しても高速に素数判定を行えます。試行回数を増やせば判定の信頼性が向上するため、暗号技術などで大きな素数を扱う場面で役立つ手法です。ただし本アルゴリズムは確率的な判定である点には注意が必要です。