Project Euler 11

http://projecteuler.net/index.php?section=problems&id=11

縦・横・斜めに次々に数を出すジェネレータを作る。4つ取るのはHaskellのような自作のtakeで。どちらに進むかは関数オブジェクトをテンプレートで与える。
最初の座標を与えるのに、2つのrangeのproductを使うが、tupleを出すようにした。3つのtupleのときはproduct3という名前にしようと思う。zipと同じ。

itertools.h

#include <cmath>
#include "itertools.h"

using namespace std;
using namespace itertools;

vector<int> primes(1, 2);

class genPrimes : public cGenerator<int> {
    int k;
    static const int    N = 20000;
    
public:
    genPrimes() : k(0) { init_primes(); }
    int next();
    genPrimes *copy() const { return new genPrimes(); }
    
private:
    void init_primes() {
        // ふるいには範囲の上端の平方根までの素数が必要
        for(int n = 3; n <= int(sqrt((double)N)) + 1; n += 2) {
            if(is_prime(n))
                primes.push_back(n);
        }
    }
    int sieve_size() const {
        int n = primes.back();
        return ((int)sqrt((double)n) / N + 1) * N;
    }
    bool is_prime(int n) const {
        for(int k = 0; ; k++) {
            int p = primes[k];
            if(p * p > n)
                return true;
            else if(n % p == 0)
                return false;
        }
    }
};

int genPrimes::next() {
    if(k >= (int)primes.size()) {
        int s = sieve_size();
        int *a = new int[s/2];
        for(int i = 0; i < s / 2; i++)
            a[i] = true;
        int offs = primes.back() + 2;
        int max_n = offs + s - 2;
        int limit = (int)sqrt((double)max_n);
        for(int k = 1; primes[k] <= limit; k++) {
            int p = primes[k];
            int init = ((offs + p) / 2 + p - 1) / p * p - (offs + p) / 2;
            for(int k = init; k < s / 2; k += p)
                a[k] = false;
        }
        
        for(int k = 0; k < s / 2; k++) {
            if(a[k])
                primes.push_back(offs + k * 2);
        }
        delete [] a;
    }
    
    if(k < (int)primes.size()) {
        k++;
        return primes[k-1];
    }
    else {
        return next();
    }
}

const int   N = (int)2e6;

template<int m>
bool is_less(int n) {
    return n < m;
}

long long conv(int n) { return (long long)n; }

int main() {
    cout << sum(new map<long long,int>(conv,
            new takewhile<int>(is_less<N>, new genPrimes()))) << endl;
}