Project Euler 214

プロジェクトオイラー
http://projecteuler.net/index.php

Q214.
オイラーのφ関数を次々に適用していくと、
φ(5) = 4, φ(4) = 2, φ(2) = 1
などと最後には1に辿りつく鎖が出来る。この鎖の長さを上の場合4とする。4000万より小さい素数で鎖の長さが25のものの総和を求めよ。

とにかくメモリとの戦いだった。なぜなのかよくわからないが、Pythonを実行するとやたらとメモリを食う。そしてほとんど動かなくなる。メモリを1GBくらい搭載していれば問題はなかったのだろうが。しかたがなくC++で組んだら、200MBも使わずにあっさり答えが出た。
ロジックとしては、最初はよく使う素因数分解を生成する方法を用いたが、途中でφをそのまま生成するほうが断然速いことに気がついた。しかし、やたらとメモリを食う。それで、エラトステネスのふるいのように、素数pの倍数を(p-1)/p倍するという方法に切り替えた。これだと4000万までの素数を用意しなければならないが。
素数だとφ(n)はn-1だから、2で割れるのでその半分だけ&phiの計算をすればいいとも考えたが、あまり美しくないのでやめた。

#include 

#include 

#include 
using namespace std;
vector primes;
bool is_prime(int n) {

    if(n < 2)

        return false;

    

    typedef vector::const_iterator CIT;

    for(CIT q = primes.begin(); q != primes.end(); ++q) {

        int     p = *q;

        if(p * p > n)

            return true;

        else if(n % p == 0)

            return false;

    }

    return true;

}
void make_primes(int n) {

    int     m = (int)sqrt( (n + 0.5));

    for(int p = 3; p <= m; p += 2) {

        if(is_prime(p))

            primes.push_back(p);

    }

    

    int     size = (n + 15) / 16;

    int     *a = new int[size];

    if(!a) {

        cout << size << endl;

        exit(1);

    }

    for(int i = 0; i < size; i++)

        a[i] = 0;

    

    typedef vector::const_iterator CIT;

    for(CIT q = primes.begin(); q != primes.end(); ++q) {

        int     p = *q;

        if(p * p > n)

            break;

        if(p == 2)

            continue;

        for(int k = 3 * p; k < n; k += 2 * p)

            a[k>>4] |= 1 << (k & 15);

    }

    

    for(int k = (m + 1) / 2 * 2 + 1; k <= n; k += 2) {

        if((a[k>>4] & (1 << (k & 15))) == 0)

            primes.push_back(k);

    }

    cout << (int)primes.size() << endl;

}
int chain_length(int n, int *phi) {

    if(n <= 1)

        return 1;

    

    for(int k = 2; ; k++) {

        n = phi[n];

        if(n == 1)

            return k;

    }

}
int main() {

    const int   N = (int)4e7;

    const int   M = 25;

    primes.push_back(2);

    make_primes(N);

    

    int     *phi = new int[N+1];

    for(int n = 0; n <= N; n++)

        phi[n] = n;

    

    typedef vector::const_iterator CIT;

    for(CIT q = primes.begin(); q != primes.end(); ++q) {

        int     p = *q;

        for(int n = p; n <= N; n += p) {

            phi[n] /= p;

            phi[n] *= p - 1;

        }

    }

    

    long long   s = 0;

    typedef vector::const_iterator CIT;

    for(CIT q = primes.begin(); q != primes.end(); ++q) {

        int     p = *q;

        if(chain_length(p, phi) == M)

            s += p;

    }

    cout << s << endl;

    delete [] phi;

}