Solovay–Strassenの素数判定法

Solovay–Strassenの素数判定法をC++で書いて見ました。

Solovay–Strassenの素数判定法とは？

与えられた数が素数かどうかを判定する乱択アルゴリズムです。RSA暗号の実現可能性を示すために重要な役割を担った歴史的価値の高いアルゴリズムです[1]。

Solovay–Strassenの素数判定法は、オイラーの基準

a^(p-1)/2 = (a / p)  (mod p)  (1)

を利用します。ただし、pは素数です。(a / p)はルジャンドルの記号で、その値は

0　 a = 0 (mod p)の場合
1     aが法pにおける平方剰余の場合
-1　aが法pにおける平方剰余ではない場合

のように定義されます。

与えられたnが式(1)を満たせば素数、満たさなければ合成数と判定します。
左辺は繰り返し二乗法を使えば高速（nを多倍長整数とするとnの桁数に比例する時間）に計算できます。右辺は平方剰余の相互法則[2]を使えば高速に計算できます。

フェルマーテストのときと同様に、この判定法では疑陽性が発生します。
ただしaをランダムにt回選んで判定を行うことで、疑陽性の発生確率は (1/2)^t 以下になることが知られています。またこの判定法では、フェルマーテストのときのカーマイケル数のように、疑陽性の発生確率が大きくなるような特別な数はありません。

C++で書いたスパゲッティなソース

以下ソースコードです。

#include <iostream>
#include <boost/multiprecision/cpp_int.hpp>

using namespace std;

typedef boost::multiprecision::cpp_int bigint;

bigint mpow(bigint x, bigint p, bigint mod) {
    bigint ret = 1;
    while (p) {
        if (p & 1)
            ret = ret * x % mod;
        x = x * x % mod;
        p >>= 1;
    }
    return ret;
}

bigint jacobi(bigint a, bigint n) {
    if (a == 0) return 0;
    if (a == 1) return 1;

    if ((a & 1) == 0) {
        if (n % 8 == 1 || n % 8 == 7)
            return jacobi(a>>1, n);
        else
            return -jacobi(a>>1, n);
    }

    if (a >= n)
        return jacobi(a % n, n);
    
    if (a % 4 == 3 && n % 4 == 3)
        return -jacobi(n, a);
    else
        return jacobi(n, a);
}

bool SolovayStrassen(bigint n, int t = 30) {
    if (n <= 2) 
        return n == 2;

    for (int i = 0; i < t; i++) {
        bigint a = rand() % n;
        if (mpow(a, (n-1)/2, n) != (jacobi(a, n) + n) % n)
            return false;
    }

    return true;
}

int main(int argc, char **argv) {
    srand((unsigned)time(NULL));

    bigint n;
    cin >> n;
    cout << SolovayStrassen(n) << endl;

    return 0;
}

参考URL

[1] http://en.wikipedia.org/wiki/Solovay%E2%80%93Strassen_primality_test
[2] http://math.fau.edu/richman/jacobi.htm