コルセルトの判定法を使ったカーマイケル数の検出

　コルセルトの判定法を使って[0, 1000000)内のカーマイケル数を列挙してみました。

コルセルトの判定法とは、以下のような条件です。

正の合成数nがカーマイケル数であるための必要十分条件は以下のとおりである。
・nはsquare-free（平方数で割り切れない）である。
・nのすべての素因数pについて、p-1 | n-1が成り立つ。

与えられた数の素因数分解が出来ればいいので、O(√n)で判定可能です。

　[0, 1000000)までの各整数に対してコルセルトの条件が成り立つかどうかを調べると、全体でO(n√n)かかってしまいます。少し工夫してエラトステネスの篩を応用すると、O(n log log n)でカーマイケル数の列挙が出来ます。

　以下ソースコードです。

#include <iostream>
#include <cassert>
#include <map>

using namespace std;

bool carmichael[1000000];
bool prime[1000000];

#ifdef debug
map<long long, int> factorize(long long n) {
    map<long long, int> ret;
    
    for (long long i = 2; i * i <= n; i++) {
        while (n % i == 0) {
            ++ret[i];
            n /= i;
        }
    }
    
    if (n != 1)
        ret[n] = 1;
    
    return ret;
}

void korseltCriterion(int x) {
    map<long long, int> fact = factorize(x);

    for (const auto &pr : fact) {
        assert(pr.second == 1);
        assert((x-1) % (pr.first-1) == 0);
    }
}
#endif

int main(int argc, char **argv) {
    fill(carmichael, carmichael+1000000, true);
    fill(prime, prime+1000000, true);

    carmichael[0] = carmichael[1] = false;
    prime[0] = prime[1] = false;

    for (int i = 2; i < 1000000; i++) {
        if (!prime[i])
            continue;

        carmichael[i] = false;   // i is a prime number

        long long sq = (long long)i * i;
        for (int j = 2 * i; j < 1000000; j += i) {
            prime[j] = false;
            if (j % sq == 0)
                carmichael[j] = false;  // j is not square-free
            if ((j - 1) % (i - 1))
                carmichael[j] = false;  // (i-1) | (j-1) does not hold
        }
    }
    
    for (int i = 0; i < 1000000; i++) {
        if (carmichael[i]) {
            cout << i << endl;
#ifdef debug
            korseltCriterion(i);
#endif
        }
    }

    return 0;
}

以下、[0, 1000000)のカーマイケル数のリストです。
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