フェルマーテストは乱数の性能にあんまり影響されない

試しにフェルマーテストが乱数の質で性能に違いがあるかやってみました。
1千万までの合成数について、何回の判定で合成数と認識されるかというのを、100回やってみました。

1回 - 433428013
2回 - 97607
3回 - 7565
4回 - 2207
5回 - 1160

中略

213回 - 1
223回 - 2
231回 - 1
320回 - 1
平均1.001判定 分散 0.0130 偏差 0.1141

Javaの標準Randomと、メルセンヌ・ツイスターで試してみましたが、だいたいおんなじ感じでした。実行時間も、ちょっとメルセンヌ・ツイスターが遅いくらい。
結局、フェルマーテストでは、ほとんどの合成数が一回の判定で認識できてしまっているので、あまり乱数の質の影響を受けないのではないかと思います。また、200回を超えても判定できてないものがいくつかあるので、フェルマーテストでは判定回数が数百回程度では信頼性があまり高いとはいえないです。

実際には、数回の判定をすりぬけるのはカーマイケル数で、そのカーマイケル数も多くは19までの素数の倍数になっているので、19までの素数の倍数になっていないかチェックすることと、19までの素数を約数に持たないカーマイケル数を配列に持たせて判定することで、intの範囲では実用的になると思います。
まあ、intの範囲なら、√21億≒46500までの素数を求めておいて、その全部で割れるか試してみればいいだけの話ですが。

ということで、ソースはこれ。Core Duo L2500でだいたい15分くらいかかりました。
手軽に試すには、eratosの引数を100000くらい、iのループを10回くらいにするといいと思います。

import java.util.*;

public class FermetEfficient {
    public static void main(String[] args){
        r.nextLong();
        boolean[] comp = eratos(10000000);
        Map<Integer, Integer> counter = new TreeMap<Integer, Integer>();
        for(int i = 0; i < 100; ++i){
            for(int n = 9; n < comp.length; n += 2){
                if(!comp[n]) continue;//素数はとばす
                int count = fermetCheckCount(n);
                if(!counter.containsKey(count)) counter.put(count, 0);
                counter.put(count, counter.get(count) + 1);
            }
        }
        int sum = 0;//判定回数
        int total = 0;//試行回数
        for(int i : counter.keySet()){
            int c = counter.get(i);
            System.out.printf("%d回 - %d%n", i, c);
            total += c;
            sum += i * c;
        }
        double avg = (double)sum / total;//平均
        //分散を求める
        double disp = 0;
        for(int i : counter.keySet()){
            int c = counter.get(i);
            disp += (avg - i) * (avg - i) * c;
        }
        System.out.printf("平均%.3f判定 分散 %.4f  偏差 %.4f%n", 
            avg, disp / total, Math.sqrt(disp / total));
    }
    static Random r = new Random();

    /** エラトステネスの篩により合成数ならtrue/素数ならfalseの配列 */
    private static boolean[] eratos(int n){
        boolean[] ret = new boolean[n + 1];
        ret[0] = true;
        ret[1] = true;
        for(int i = 4; i <= n; i += 2){
            ret[i] = true;
        }
        for(int i = 3; i <= n; i += 2){
            if(ret[i]) continue;
            for(int j = i * 2; j <= n; j += i){
                ret[j] = true;
            }
        }
        return ret;
    }

    /** フェルマーテストで合成数と判定されるまでの回数 */
    private static int fermetCheckCount(int n){
        for(int i = 1; ; ++i){
            if(n % 2 == 0) return i;//偶数なら合成数
            int a = r.nextInt(n - 2) + 2;//2以上n未満の値
            int g = (int) gcd(n, a);
            if(g != 1) return i;//最大公約数がみつかったなら合成数
            if(modPow(a, n - 1, n) != 1) return i;
        }
    }

    /** (a ^ b) % m を計算 */
    static int modPow(int a, int b, int m){
        if(a == 0) return 0;

        int result = 1;
        a %= m;
        while(b != 0){
            if(b % 2 == 1){
                result = (int)(((long)result * a) % m);
            }
            a = (int)(((long)a * a) % m);
            b = b >> 1;
        }

        return result;
    }

    /** 最大公約数を求める(a > b) */
    private static long gcd(long a, long b){
        if(b == 0) return a;
        long m = a % b;

        return gcd(b, m);
    }
}