カーマイケル数を求めてみる

やっぱ、どの数がカーマイケル数か気になるので、1千万までのカーマイケル数を求めてみました。全部で105個です。

561
1105	1729	2465	2821	6601	8911
10585	15841	29341	41041	46657	52633	62745	63973	75361
101101	115921	126217	162401	172081	188461	252601	278545	294409
314821	334153	340561	399001	410041	449065	488881	512461	530881
552721	656601	658801	670033	748657	825265	838201	852841	997633
1024651	1033669	1050985	1082809	1152271	1193221	1461241	1569457	1615681
1773289	1857241	1909001	2100901	2113921	2433601	2455921	2508013	2531845
2628073	2704801	3057601	3146221	3224065	3581761	3664585	3828001	4335241
4463641	4767841	4903921	4909177	5031181	5049001	5148001	5310721	5444489
5481451	5632705	5968873	6049681	6054985	6189121	6313681	6733693	6840001
6868261	7207201	7519441	7995169	8134561	8341201	8355841	8719309	8719921
8830801	8927101	9439201	9494101	9582145	9585541	9613297	9890881

ソースはこんな感じ。10分くらいかかりました。

public class CarmichaelNumber {
    public static void main(String[] args){
        boolean[] comp = eratos(10000000);
        for(int i = 2; i < comp.length; ++i){
            if(!comp[i]) continue;
            if(carmichelCheck(i)) System.out.println(i);
        }
    }

    /** エラトステネスの篩により合成数なら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 boolean  carmichelCheck(int n){
        for(int a = 2; a < n ; ++a){
            int g = (int) gcd(n, a);
            if(g != 1) continue;
            if(modPow(a, n - 1, n) != 1) return false;
        }
        return true;
    }

    /** (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);
    }
}