"The Genuine Sieve of Eratosthenes" in C++14 The following C++ code implements the "Genuine Sieve of Eratosthenes" algorithm as described in Melissa O'Neill's classic paper.  On my MacBook it computes the first 1,000,000 primes in about 11 seconds.

$ ./a.out | head -1000000 | tail -1
15485863
$


Just looking for general code review comments here. At least two parts of sieverator smell really bad to me, but I'm not sure of the proper way to fix them while still preserving the general "STLishness" of this code. For example, I really want to keep using the ranged for-loop in main.

#include &lt;cassert&gt;
#include &lt;functional&gt;
#include &lt;iostream&gt;
#include &lt;memory&gt;
#include &lt;queue&gt;
#include &lt;utility&gt;

template&lt;typename Int&gt;
class iotarator {
    Int value = 0;
    Int step = 1;
public:
    explicit iotarator() = default;
    explicit iotarator(Int v) : value(v) {}
    explicit iotarator(Int v, Int s) : value(v), step(s) {}
    Int operator*() const { return value; }
    iotarator&amp; operator++() { value += step; return *this; }
    iotarator operator++(int) { auto ret = *this; ++*this; return ret; }

    bool operator==(const iotarator&amp; rhs) const {
        return value == rhs.value &amp;&amp; step == rhs.step;
    }
    bool operator!=(const iotarator&amp; rhs) const { return !(*this == rhs); }
};

template&lt;class Int&gt;
class sieverator {
    struct erased_iterator {
        virtual Int dereference() = 0;
        virtual void increment() = 0;
    };
    template&lt;class It&gt;
    class derived_iterator : public erased_iterator {
        It it;
    public:
        derived_iterator(It it) : it(std::move(it)) {}
        Int dereference() override { return *it; }
        void increment() override { ++it; }
    };

    Int m_current;
    std::unique_ptr&lt;erased_iterator&gt; m_ptr;

    explicit sieverator() {}  // used by .end()
public:
    template&lt;class It&gt;
    explicit sieverator(It it) :
        m_current(*it),
        m_ptr(std::make_unique&lt;derived_iterator&lt;It&gt;&gt;(std::move(it)))
    {}
    sieverator begin() { return std::move(*this); }
    sieverator end() const { return sieverator{}; }
    bool operator==(const sieverator&amp;) const { return false; }
    bool operator!=(const sieverator&amp;) const { return true; }

    Int operator*() const {
        return m_current;
    }

    sieverator&amp; operator++() {
        cross_off_multiples_of_prime(m_current);
        do {
            m_ptr-&gt;increment();
            m_current = m_ptr-&gt;dereference();
        } while (is_already_crossed_off(m_current));
        return *this;
    }

    sieverator&amp; operator++(int) = delete;

private:
    struct pair {
        Int next_crossed_off_value;
        Int prime_increment;
        explicit pair(Int a, Int b) : next_crossed_off_value(a), prime_increment(b) {}
        bool operator&lt;(const pair&amp; rhs) const {
            return next_crossed_off_value &lt; rhs.next_crossed_off_value ? true
                 : next_crossed_off_value &gt; rhs.next_crossed_off_value ? false
                 : prime_increment &lt; rhs.prime_increment;
        }
        bool operator&gt;(const pair&amp; rhs) const { return rhs &lt; *this; }
    };
    template&lt;class T&gt; using min_heap = std::priority_queue&lt;T, std::vector&lt;T&gt;, std::greater&lt;&gt;&gt;;
    min_heap&lt;pair&gt; m_pq;

    bool is_already_crossed_off(Int value) {
        if (value != m_pq.top().next_crossed_off_value) {
            return false;
        } else {
            do {
                auto x = m_pq.top();
                m_pq.pop();
                m_pq.emplace(x.next_crossed_off_value + x.prime_increment, x.prime_increment);
            } while (value == m_pq.top().next_crossed_off_value);
            return true;
        }
    }

    void cross_off_multiples_of_prime(Int value) {
        m_pq.emplace(value * value, value);
    }
};

int main()
{
    iotarator&lt;__int128_t&gt; iota(2);
    sieverator&lt;__int128_t&gt; sieve(iota);
    for (int p : sieve) {
        std::cout &lt;&lt; p &lt;&lt; std::endl;
    }
}
At least to me, this seems a little like another instance of the same basic problem pointed out in linked paper: although there are a few parts that I guess are sort of neat, overall it strikes me as a huge amount of code to solve a trivial problem, and (at least based on the numbers you're posting) apparently runs really slowly.

It doesn't attempt to make much in the way of style points, but lets consider sort of a reference implementation--roughly the simplest thing that works to some degree:

#include &lt;vector&gt;
#include &lt;iostream&gt;
#include &lt;locale&gt;

unsigned long primes = 0;

int main() {
    int number = 20'000'000;
    std::vector&lt;bool&gt; sieve(number,false);
    sieve[0] = sieve[1] = true;

    for(int i = 2; i&lt;number; i++) {
        if(!sieve[i]) {
            ++primes;
            for (int temp = 2*i; temp&lt;number; temp += i)
                sieve[temp] = true;
        }
    }
    std::cout.imbue(std::locale(""));
    std::cout &lt;&lt; "found: " &lt;&lt; primes &lt;&lt; " Primes\n";
    return 0;
}


On my notebook (4 or 5 years old, and wasn't tremendously fast when it was new) this finds 1,270,607 primes in 147 milliseconds (so in round numbers, it's running something like 75 or 80 times as fast as the code in the question). Even compared to @Toby's optimized version, it's still something like 20 times faster.

Don't get me wrong: I'm all for using nice algorithms and making code generic where it's reasonable--and I'm pretty sure there are ways to do that in this case without paying much of a penalty. But, at least to me, it seems to me as if the code you presented is paying a pretty high penalty--and given its much greater length and overall complexity, I'm not convinced we're getting any real improvement in return for that speed penalty either.

I'd re-emphasize that this is not really optimized code, or anything like that either--we can do quite a bit better by taking the square root into account, using a segmented sieve, and switching to something like a sieve of Atkins, and using a more specialized bit vector instead of std::vector&lt;bool&gt;. If we did those, we could probably expect to improve speed by around another order of magnitude (and depending on the number of primes we decided to find, maybe even more than that).