フェルマーテストで素数判定

　ちょっと気晴らしにフェルマーテストを書いてみました。

この前boostの多倍長ライブラリを入れたので、C++で書いてみました。テストの式の微妙な違いで大きな差が出たのでちょっとメモしておきます。

a^p = a (mod p)    (1)

のようなテストをすると、大量のカーマイケル数が陽性になりました。ちなみに(1)は正しいフェルマーテストのチェック式ではなくて、何も見ずにプログラムを書いてみたときに使った式です。

で疑陽性多すぎるなと思って、よくよく調べてみると正しい式は(2)でした。

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

上式でテストすると、疑陽性は少なくなりました。[1, 1000000)まででテストしましたが、疑陽性は出ませんでした（各数に対するテスト反復回数=100）。(1)と(2)って同じように見えるけど、(2)の方が厳しい条件ですね。

例えばカーマイケル数561は、
237⁵⁶¹ = 237 (mod 561) のようになります。
しかし、 237⁵⁶⁰ = 375 (mod 561)
なので、(1)のようなテストをすると疑陽性が発生し、(2)のようなテストをすると正しく陰性になるというわけです。

んー、ちゃんと定義は見た方がいいな。ということで、カーマイケル数の定義を確認してみました。

wikipediaを見てみると、
合成数 n がカーマイケル数であるとは、自身と互いに素である任意の自然数 a に対し、
a^n-1 = 1 (mod n)    (3)
を満たすことをいう。

とあります。

カーマイケル数は互いに素でないaに対しては(3)を満たすとは限らないということか。。

確かに、例で見た237と561はどちらも3で割れるので、互いに素ではないです。

んーということは、カーマイケル数nに対するフェルマーテストを行う際に、nと互いに素ではないaがうまく選ばれれば疑陽性の発生頻度が少なくなる可能性があると考えてもよさそうです。

Fermat witnessの濃度を考えると実験で設定した反復回数=100ってちょっと多すぎやろと思いましたが、カーマイケル数に対するテストでカーマイケル数と互いに素ではない数がたくさん選ばれるようにするという意味で100という数はなかなかいいチョイスだったのかもしれません。。

今回のテストで使った範囲[0, 1000000)内のカーマイケル数のtotient functionを計算してみると、
Φ(561) = 320
Φ(1105) = 768

Φ(1729) = 1296

Φ(2465) = 1792
......
Φ(488881) = 466560
......

のように。

488881の場合は、coprimeじゃない数は5%くらいの確率でしか出ない。んー、やはり100という設定値は絶妙なチョイスだったようだ。

と、意外と発見があったフェルマーテストの実装でした。
コードを載せておきます。

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

using namespace std;

typedef boost::multiprecision::cpp_int bigint;

bigint pow(bigint x, bigint p, bigint mod) {
    bigint ret = 1;

    while (p > 0) {
        if (p & 1)
            ret = ret * x % mod;
        x = x * x % mod;
        p >>= 1;
    }

    return ret;
}

bool fermatTest(bigint p, int k = 100) {
    if (p < 2)
        return false;
    
    srand((unsigned)time(NULL));
    for (int i = 0; i < k; i++) {
        bigint a = rand() % (p-1) + 1;
        if (pow(a, p-1, p) != 1)  // this condition is better than pow(a, p, p) != a
            return false;
    }

    return true;
}

#ifndef TEST
int main() {
    bigint n;
    cin >> n;

    cout << fermatTest(n) << endl;
    
    return 0;
}
#else

bigint primalityTest(bigint x) {
    for (bigint i = 2; i * i <= x; i++) {
        if (x % i == 0)
            return false;
    }

    return x > 1;
}

int main() {
    bool res = true;
    for (int i = 1; i < 1000000; i++) {
        if (fermatTest(i) != primalityTest(i)) {
            res = false;
            cout << i << endl;
        }
    }
    if (res)
        cout << "passed." << endl;

    return 0;
}
#endif