enum_prime_number.cpp

enum_prime_number.cpp

#include <algorithm>
#include <vector>
#include <deque>
#include <iostream>
#include <utility>
#include <limits>
#include <array>
#include <functional>
#include <cstdint>
#include <cmath>
#include <cstdio>

namespace project_60{
    template<class UInt = std::uint32_t, class DoubleUInt = std::uint64_t, class DoubleInt = std::int64_t, DoubleUInt BitNum = sizeof(UInt) * 8>
    struct quint{
        template<DoubleUInt Count>
        struct make_bitmask{
            static const DoubleUInt value = (static_cast<DoubleUInt>(1) << Count) | make_bitmask<Count - static_cast<DoubleUInt>(1)>::value;
        };

        template<>
        struct make_bitmask<0>{
            static const DoubleUInt value = static_cast<DoubleUInt>(1);
        };

        static const DoubleUInt base_mask = make_bitmask<BitNum - 1>::value;
        UInt data[4];

        quint(){
            for(int i = 0; i < 4; ++i){ data[i] = 0; }
        }

        quint(const quint &other){
            for(int i = 0; i < 4; ++i){ data[i] = other.data[i]; }
        }

        quint(quint &&other){
            for(int i = 0; i < 4; ++i){ data[i] = other.data[i]; }
        }

        ~quint() = default;

        quint(int x){
            data[0] = static_cast<UInt>(x);
            data[1] = data[2] = data[3] = 0;
        }

        quint(UInt x){
            data[0] = x;
            data[1] = data[2] = data[3] = 0;
        }

        quint &operator =(const quint &other){
            for(int i = 0; i < 4; ++i){ data[i] = other.data[i]; }
            return *this;
        }

        bool operator <(const quint &other) const{
            for(int i = 3; i > -1; --i){
                if(data[i] < other.data[i]){ return true; }
            }
            return false;
        }

        bool operator >(const quint &other) const{
            return other < *this;
        }

        bool operator <=(const quint &other) const{
            for(int i = 3; i > -1; --i){
                if(data[i] < other.data[i]){ return true;  }
                if(data[i] > other.data[i]){ return false; }
            }
            return true;
        }

        bool operator >=(const quint &other) const{
            return other <= *this;
        }

        bool operator ==(const quint &other) const{
            for(int i = 0; i < 4; ++i){
                if(data[i] != other.data[i]){ return false; }
            }
            return true;
        }

        bool operator !=(const quint &other) const{
            return !(*this == other);
        }

        quint add(const quint &other) const{
            quint result;
            DoubleUInt c = 0;
            for(int i = 0; i < 4; ++i){
                DoubleUInt v = static_cast<DoubleUInt>(data[i]) + static_cast<DoubleUInt>(other.data[i]) + c;
                result.data[i] = static_cast<UInt>(v & base_mask);
                c = v >> BitNum;
            }
            return result;
        }

        quint sub(const quint &other) const{
            quint result;
            for(int i = 0; i < 4; ++i){ result.data[i] = ~other.data[i]; }
            DoubleUInt c = 1;
            for(int i = 0; i < 4; ++i){
                DoubleUInt v = static_cast<DoubleUInt>(data[i]) + static_cast<DoubleUInt>(result.data[i]) + c;
                result.data[i] = static_cast<UInt>(v & base_mask);
                c = v >> BitNum;
            }
            return result;
        }

        quint mul(const quint &other) const{
            UInt extended_data[8] = { 0 };
            for(int i = 0; i < 4; ++i){
                for(int j = 0; j < 4; ++j){
                    DoubleUInt c = static_cast<DoubleUInt>(data[i]) * static_cast<DoubleUInt>(other.data[j]);
                    for(int k = 0; i + j + k < 4; ++k){
                        DoubleUInt v = static_cast<DoubleUInt>(extended_data[i + j + k]) + c;
                        extended_data[i + j + k] = static_cast<UInt>(v & base_mask);
                        c = v >> BitNum;
                    }
                }
            }
            quint result;
            for(int i = 0; i < 4; ++i){
                result.data[i] = extended_data[i];
            }
            return result;
        }

        quint div(const quint &other) const{
            quint r = *this, result;
            int n = degree(), m = other.degree(), i = n - m;
            UInt u = other.data[m - 1];
            while(i >= 0){
                int rd = r.degree();
                if(rd == m + i){
                    result.data[i] = r.data[rd - 1] / u;
                    r = r.sub(other.mul(result.data[i]));
                }else{
                    result.data[i] = 0;
                }
                --i;
            }
            return result;
        }

        quint rem(const quint &other) const{
            return sub(other.mul(div(other)));
        }

        int degree() const{
            int i;
            for(i = 4; i > 0; --i){
                if(data[i - 1] > 0){ break; }
            }
            return i;
        }

        UInt lc() const{ return data[degree() - 1]; }
    };

    template<class UInt, class DoubleUInt, class DoubleInt, DoubleUInt BitNum>
    quint<UInt, DoubleUInt, DoubleInt, BitNum> operator +(
        const quint<UInt, DoubleUInt, DoubleInt, BitNum> &lhs,
        const quint<UInt, DoubleUInt, DoubleInt, BitNum> &rhs
    ){ return lhs.add(rhs); }

    template<class UInt, class DoubleUInt, class DoubleInt, DoubleUInt BitNum>
    quint<UInt, DoubleUInt, DoubleInt, BitNum> operator -(
        const quint<UInt, DoubleUInt, DoubleInt, BitNum> &lhs,
        const quint<UInt, DoubleUInt, DoubleInt, BitNum> &rhs
    ){ return lhs.sub(rhs); }

    template<class UInt, class DoubleUInt, class DoubleInt, DoubleUInt BitNum>
    quint<UInt, DoubleUInt, DoubleInt, BitNum> operator *(
        const quint<UInt, DoubleUInt, DoubleInt, BitNum> &lhs,
        const quint<UInt, DoubleUInt, DoubleInt, BitNum> &rhs
    ){ return lhs.mul(rhs); }

    template<class UInt, class DoubleUInt, class DoubleInt, DoubleUInt BitNum>
    quint<UInt, DoubleUInt, DoubleInt, BitNum> operator /(
        const quint<UInt, DoubleUInt, DoubleInt, BitNum> &lhs,
        const quint<UInt, DoubleUInt, DoubleInt, BitNum> &rhs
    ){ return lhs.div(rhs); }

    template<class UInt, class DoubleUInt, class DoubleInt, DoubleUInt BitNum>
    quint<UInt, DoubleUInt, DoubleInt, BitNum> operator %(
        const quint<UInt, DoubleUInt, DoubleInt, BitNum> &lhs,
        const quint<UInt, DoubleUInt, DoubleInt, BitNum> &rhs
    ){ return lhs.rem(rhs); }

#define PROJECT_60_QUINT_OPERATOR_OVERLOAD_FOR_BUILT_IN_TYPE(built_in_type, symbol, fn) \
    template<class UInt, class DoubleUInt, class DoubleInt, DoubleUInt BitNum> \
    quint<UInt, DoubleUInt, DoubleInt, BitNum> operator symbol( \
        const quint<UInt, DoubleUInt, DoubleInt, BitNum> &lhs, \
        const built_in_type &rhs \
    ){ return lhs.fn(quint<UInt, DoubleUInt, DoubleInt, BitNum>(static_cast<unsigned>(rhs))); } \
    template<class UInt, class DoubleUInt, class DoubleInt, DoubleUInt BitNum> \
    quint<UInt, DoubleUInt, DoubleInt, BitNum> operator symbol( \
        const built_in_type &lhs, \
        const quint<UInt, DoubleUInt, DoubleInt, BitNum> &rhs \
    ){ return quint<UInt, DoubleUInt, DoubleInt, BitNum>(static_cast<unsigned>(lhs)).fn(rhs); }

    PROJECT_60_QUINT_OPERATOR_OVERLOAD_FOR_BUILT_IN_TYPE(int, +, add);
    PROJECT_60_QUINT_OPERATOR_OVERLOAD_FOR_BUILT_IN_TYPE(int, -, sub);
    PROJECT_60_QUINT_OPERATOR_OVERLOAD_FOR_BUILT_IN_TYPE(int, *, mul);
    PROJECT_60_QUINT_OPERATOR_OVERLOAD_FOR_BUILT_IN_TYPE(int, /, div);
    PROJECT_60_QUINT_OPERATOR_OVERLOAD_FOR_BUILT_IN_TYPE(int, %, rem);

    PROJECT_60_QUINT_OPERATOR_OVERLOAD_FOR_BUILT_IN_TYPE(unsigned, +, add);
    PROJECT_60_QUINT_OPERATOR_OVERLOAD_FOR_BUILT_IN_TYPE(unsigned, -, sub);
    PROJECT_60_QUINT_OPERATOR_OVERLOAD_FOR_BUILT_IN_TYPE(unsigned, *, mul);
    PROJECT_60_QUINT_OPERATOR_OVERLOAD_FOR_BUILT_IN_TYPE(unsigned, /, div);
    PROJECT_60_QUINT_OPERATOR_OVERLOAD_FOR_BUILT_IN_TYPE(unsigned, %, rem);

    template<class UInt, class DoubleUInt, class DoubleInt, DoubleUInt BitNum>
    std::ostream &operator <<(
        std::ostream &os,
        const quint<UInt, DoubleUInt, DoubleInt, BitNum> &rhs
    ){
        UInt sum = 0;
        for(int i = 0; i < 4; ++i){ sum |= rhs.data[i]; }
        if(sum == 0){
            os << '0';
            return os;
        }
        char stack[40] = { 0 };
        int stack_count = 0;
        quint<UInt, DoubleUInt, DoubleInt, BitNum> t = rhs;
        while(
            [&](const quint<UInt, DoubleUInt, DoubleInt, BitNum> &v){
                sum = 0; for(int i = 0; i < 4; ++i){ sum |= v.data[i]; } return sum > 0;
            }(t)
        ){
            char a[4] = { 0 };
            std::sprintf(a, "%d", t.rem(10).data[0]);
            t = t / 10;
            stack[stack_count++] = a[0];
        }
        while(stack_count){
            char a[2] = { stack[stack_count--], '\0' };
            os << a;
        }
        return os;
    }

    //using qint = quint<std::uint32_t, std::uint64_t, std::int64_t, 32>;
    using qint = std::uint_fast32_t;
    struct end_of_elf{};

    std::vector<qint> prime_table, composite_prime_table;

    void elf(std::size_t limit_num){
        std::vector<qint> &prime_table(prime_table);
        prime_table.reserve(limit_num);
        composite_prime_table.reserve(limit_num / 2);

        std::size_t skip_idx = 0;
        qint line_limit = 1, n_line_limit, next;

        auto search = [&](
            const qint &x,
            std::size_t idx_range_first,
            std::size_t idx_range_last,
            std::function<void(const qint&)> push_back
        ){
            for(std::size_t i = idx_range_first; i < idx_range_last; ++i){
                const qint &y = prime_table[i];
                if(y * y > x){
                    push_back(x);
                    if(prime_table.size() >= limit_num){
                        throw end_of_elf();
                    }
                    break;
                }
                if((x % y) == 0){ break; }
            }
        };

        std::function<void(std::size_t, std::size_t, qint)> composite_search;
        composite_search = [&](
            std::size_t idx_range_first,
            std::size_t prime_idx,
            qint x
        ){
            {
                qint y = x * prime_table[prime_idx];
                {
                    qint z = y + n_line_limit;
                    if(z > next){ return; }
                    search(
                        z,
                        idx_range_first,
                        prime_table.size(),
                        [&](const qint &p){
                            auto iter = std::lower_bound(composite_prime_table.begin(), composite_prime_table.end(), p);
                            if(iter == composite_prime_table.end() || *iter != p){
                                composite_prime_table.insert(iter, p);
                            }
                        }
                    );
                }
                composite_search(idx_range_first, prime_idx, y);
            }
            if(prime_idx + 1 < prime_table.size()){
                composite_search(idx_range_first, prime_idx + 1, x);
                composite_search(idx_range_first, prime_idx + 1, prime_table[prime_idx + 1]);
            }
        };

        // bootstrap prime numbers for line-limit, p_1 = 2, p_2 = 3.
        prime_table.push_back(qint(2));
        prime_table.push_back(qint(3));

        // bootstrap prime number, for next line-limit, p_3 = 5.
        prime_table.push_back(qint(5));

        for(; skip_idx < prime_table.size() - 1; ++skip_idx){
            line_limit = line_limit * prime_table[skip_idx];
        }

        for(; ; ){
            next = line_limit * prime_table[skip_idx];
            std::size_t s = prime_table.size(), composite_counter = 0;

            auto push_back = [&](const qint &x){
                while(composite_counter < composite_prime_table.size() && composite_prime_table[composite_counter] < x){
                    search(
                        composite_prime_table[composite_counter++],
                        skip_idx,
                        prime_table.size(),
                        [&](const qint &y){ prime_table.push_back(y); }
                    );
                }
                if(composite_counter < composite_prime_table.size() && composite_prime_table[composite_counter] == x){
                    ++composite_counter;
                }
                prime_table.push_back(x);
            };

            for(n_line_limit = line_limit; n_line_limit < next; n_line_limit = n_line_limit + line_limit){
                search(n_line_limit + 1, skip_idx, prime_table.size(), push_back);

                for(std::size_t a = skip_idx; a < s; ++a){
                    composite_search(a, a, prime_table[a]);
                    search(n_line_limit + prime_table[a], skip_idx, prime_table.size(), push_back);
                }
            }

            composite_prime_table.clear();
            ++skip_idx;
            line_limit = next;
            next = line_limit * prime_table[skip_idx];
        }

        return;
    }
}

#include <boost/timer.hpp>

int main(){
    boost::timer t;
    try{
        project_60::elf(512);
    }catch(project_60::end_of_elf){
        // end...
    }
    
    std::cout << project_60::prime_table.back() << std::endl;
    std::cout << t.elapsed() << std::endl;
    
    char a;
    std::cin >> a;
    return 0;
}