Project Euler 9

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

いつものようにピタゴラス数の式を使うと、

a + b + c = 2lm(m + n)

これが1000になる。約数を出すオブジェクトを作ればよいが、これがなかなか大変だった。

itertools.h

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

using namespace std;
using namespace itertools;

typedef tuple<int,int>  Fact;

class factorize : public cGenerator<Fact> {
    int     n;
    Fact    f;
    factorize   *fs;
    bool    first;
    
public:
    factorize(int n) : n(n) {
        first = true;
        
        for(int p = 2; p * p <= n; p++) {
            int e = 0;
            while(n % p == 0) {
                n /= p;
                e++;
            }
            if(e > 0) {
                f = Fact(p, e);
                if(n > 1)
                    fs = new factorize(n);
                else
                    fs = NULL;
                return;
            }
        }
        f = Fact(n, 1);
        fs = NULL;
    }
    
    Fact next() {
        try {
            if(first) {
                first = false;
                return f;
            }
            else {
                if(fs == NULL)
                    throw(1);
                return fs->next();
            }
        }
        catch(int e) {
            delete fs;
            throw(1);
        }
    }
    
    factorize *copy() const { return new factorize(n); }
};

typedef vector<Fact>        Facts;
ostream& operator<<(ostream& os, const Facts& f);

class divisors : public cGenerator<tuple<Facts,Facts>> {
    Facts::iterator begin;
    Facts::iterator end;
    divisors    *g;
    int         p, e;
    int         e_max;
    tuple<Facts,Facts>      f;
    
public:
    divisors(Facts::iterator b, Facts::iterator e) : begin(b), end(e) {
        if(begin == end) {
            g = NULL;
        }
        else {
            g = new divisors(begin + 1, end);
            p = fst(*begin);
            e_max = snd(*begin);
        }
        this->e = 0;
    }
    
    tuple<Facts,Facts> next() {
        if(g == NULL) {
            if(e == 0) {
                e++;
                return tuple<Facts,Facts>(Facts(),Facts());
            }
            else {
                throw(1);
            }
        }
        else {
            if(e == 0)
                f = g->next();
            else if(e == 1)
                get<0>(f).insert(fst(f).begin(), Fact(p, 1));
            else
                get<1>(fst(f).front()) = e;
            
            if(e == 0)
                get<1>(f).insert(snd(f).begin(), Fact(p, e_max));
            else if(e == e_max)
                get<1>(f).erase(snd(f).begin());
            else
                get<1>(snd(f).front()) = e_max - e;
            
            if(e == e_max)
                e = 0;
            else
                e++;
            
            return f;
        }
    }
    
    divisors *copy() const { return new divisors(begin, end); }
};

int value_f(const Facts& f) {
    int     v = 1;
    for(int i = 0; i < (int)f.size(); i++) {
        int p = get<0>(f[i]);
        int e_max = get<1>(f[i]);
        for(int e = 0; e < e_max; e++)
            v *= p;
    }
    return v;
}

typedef tuple<int,int,int>  triplet;

class divs : public cGenerator<triplet> {
    divisors    *g1;
    divisors    *g2;
    int         n1;
    
public:
    divs(int n) {
        Facts   f = list(new factorize(n));
        g1 = new divisors(f.begin(), f.end());
        next_f2();
    }
    
    divs(divisors *g1, divisors *g2, int n1) : g1(g1), g2(g2), n1(n1) { }
    
    tuple<int,int,int> next() {
        try {
            tuple<Facts,Facts>  ff = g2->next();
            return triplet(n1, value_f(fst(ff)), value_f(snd(ff)));
        }
        catch(int e) {
            delete g2;
            next_f2();
            tuple<Facts,Facts>  ff = g2->next();
            return triplet(n1, value_f(fst(ff)), value_f(snd(ff)));
        }
    }
    
    divs *copy() const { return new divs(g1, g2, n1); }
    
private:
    void next_f2() {
        tuple<Facts,Facts>  ff = g1->next();
        n1 = value_f(fst(ff));
        Facts   f2 = snd(ff);
        g2 = new divisors(f2.begin(), f2.end());
    }
};

template<class T>
T gcd(T a, T b) {
    if(a % b == 0)
        return b;
    else
        return gcd(b, a % b);
}

bool is_valid(triplet d) {
    int l = get<0>(d);
    int m = get<1>(d);
    int n = get<2>(d) - m;
    return 0 < n && n < m && (m + n) % 2 == 1 && gcd(m, n) == 1;
}

int abc(triplet d) {
    int l = get<0>(d);
    int m = get<1>(d);
    int n = get<2>(d) - m;
    int a = l * (m * m - n * n);
    int b = 2 * l * m * n;
    int c = l * (m * m + n * n);
    return a * b * c;
}

const int   N = 1000;

int main() {
    cout << head(new map<int,triplet>(abc,
                new filter<triplet>(is_valid, new divs(N / 2)))) << endl;
}