primes.h

素数や素因数分解に関するライブラリ。VC10β2で先取りしたC++0xの機能を用いている。
破壊的な変更がされる場合は別にエントリーを立てる。

//
// primes.h
// by M.Inamori
// written 3/25/10
// revised 4/16/10

#pragma once
#include "itertools.h"

namespace primes {
    std::vector<int>    primes(1, 2);
    
    class genPrimes {
        int k;
        static const int    N = 20000;
        
    public:
        genPrimes(int k = 0) : k(k) { init_primes(); }
        int next();
        bool exists_next() { return true; }
        
    private:
        void init_primes() {
            if((int)primes.size() > 2)
                return;
            
            // ふるいには範囲の上端の平方根までの素数が必要
            for(int n = 3; n <= int(sqrt((double)N)) * 2 + 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();
        }
    }
    
    bool is_prime(int n) {
        using namespace itertools;
        if(n < 2) return false;
        return all([n] (int p) { return n % p != 0; },
                    takewhile([n] (int p) { return p * p <= n; },
                    genPrimes()));
    }
    
    typedef std::tuple<int,int> Fact;
    typedef std::vector<Fact>   Facts;
    
    std::ostream& operator<<(std::ostream& os, const Fact& f) {
        int     p = itertools::fst(f);
        int     e = itertools::snd(f);
        os << p;
        if(e > 1)
            os << "^" << e;
        return os;
    }
    
    std::ostream& operator<<(std::ostream& os, const Facts& f) {
        if(f.size() > 0) {
            for(int k = 0; k < (int)f.size() - 1; k++) {
                os << f[k] << " * ";
            }
            os << f.back();
        }
        else {
            os << 1;
        }
        return os;
    }
    
    class factorize {
        int     n;
        Fact    f;
        factorize   *fs;
        bool    first;
        
    public:
        factorize(int n, genPrimes *pg = nullptr) : n(n), first(true) {
            if(pg == nullptr)
                pg = new genPrimes();
            
            int k = 0;
            int p = pg->next();
            while(p * p <= n) {
                k++;
                int e = 0;
                while(n % p == 0) {
                    n /= p;
                    e++;
                }
                if(e > 0) {
                    f = Fact(p, e);
                    if(n > 1)
                        fs = new factorize(n, new genPrimes(k));
                    else
                        fs = nullptr;
                    delete pg;
                    return;
                }
                p = pg->next();
            }
            f = Fact(n, 1);
            fs = nullptr;
            delete pg;
        }
        ~factorize() {
            if(fs != nullptr)
                delete fs;
        }
        
        Fact next() {
            if(first) {
                first = false;
                return f;
            }
            else {
                if(fs == nullptr)
                    throw("error : no more elements in factorize.");
                return fs->next();
            }
        }
        bool exists_next() {
            if(first)
                return true;
            else if(fs == nullptr)
                return false;
            else
                return fs->exists_next();
        }
    };
    
    class divisors {
        Facts::iterator begin;
        Facts::iterator end;
        divisors    *g;
        int         p, e;
        int         e_max;
        std::tuple<Facts,Facts>     f;
        
    public:
        divisors() { }
        divisors(Facts::iterator b, Facts::iterator e) : begin(b), end(e) {
            if(begin == end) {
                g = nullptr;
            }
            else {
                g = new divisors(begin + 1, end);
                p = std::get<0>(*begin);
                e_max = std::get<1>(*begin);
            }
            this->e = 0;
        }
        
        std::tuple<Facts,Facts> next() {
            using namespace std;
            
            if(g == nullptr) {
                if(e == 0) {
                    e++;
                    return make_tuple(Facts(),Facts());
                }
                else {
                    throw("error : no more elements in divisors.");
                }
            }
            else {
                if(e == 0)
                    f = g->next();
                else if(e == 1)
                    get<0>(f).insert(get<0>(f).begin(), Fact(p, 1));
                else
                    get<1>(get<0>(f).front()) = e;
                
                if(e == 0)
                    get<1>(f).insert(get<1>(f).begin(), Fact(p, e_max));
                else if(e == e_max)
                    get<1>(f).erase(get<1>(f).begin());
                else
                    get<1>(get<1>(f).front()) = e_max - e;
                
                if(e == e_max)
                    e = 0;
                else
                    e++;
                
                return f;
            }
        }
        bool exists_next() {
            if(g == nullptr)
                return e == 0;
            else if(e == 0)
                return g->exists_next();
            else
                return true;
        }
    };
    
    int value_f(const Facts& f) {
        int     v = 1;
        for(int i = 0; i < (int)f.size(); i++) {
            int p = std::get<0>(f[i]);
            int e_max = std::get<1>(f[i]);
            for(int e = 0; e < e_max; e++)
                v *= p;
        }
        return v;
    }
    
    Facts mul_f(const Facts& f1, const Facts& f2) {
        Facts   f;
        int     k1 = 0;
        int     k2 = 0;
        while(true) {
            if(k1 < (int)f1.size()) {
                int p1 = itertools::fst(f1[k1]);
                int e1 = itertools::snd(f1[k1]);
                if(k2 < (int)f2.size()) {
                    int p2 = itertools::fst(f2[k2]);
                    int e2 = itertools::snd(f2[k2]);
                    if(p1 == p2) {
                        f.push_back(Fact(p1, e1 + e2));
                        k1++;
                        k2++;
                    }
                    else if(p1 < p2) {
                        f.push_back(f1[k1]);
                        k1++;
                    }
                    else {
                        f.push_back(f2[k2]);
                        k2++;
                    }
                }
                else {
                    f.push_back(f1[k1]);
                    k1++;
                }
            }
            else {
                if(k2 < (int)f2.size()) {
                    f.push_back(f2[k2]);
                    k2++;
                }
                else {
                    break;
                }
            }
        }
        return f;
    }
    
    // 割り算の結果は自然数となる仮定で
    Facts div_f(const Facts& f1, const Facts& f2) {
        Facts   f;
        int     k2 = 0;
        for(int k1 = 0; k1 < (int)f1.size(); k1++) {
            int p1 = itertools::fst(f1[k1]);
            int e1 = itertools::snd(f1[k1]);
            if(k2 < (int)f2.size()) {
                int p2 = itertools::fst(f2[k2]);
                int e2 = itertools::snd(f2[k2]);
                if(p1 == p2) {
                    if(e1 > e2)
                        f.push_back(Fact(p1, e1 - e2));
                    k2++;
                }
                else {
                    f.push_back(f1[k1]);
                }
            }
            else {
                f.push_back(f1[k1]);
            }
        }
        return f;
    }
}