# “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 <cassert>
#include <functional>
#include <iostream>
#include <memory>
#include <queue>
#include <utility>

template<typename Int>
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& operator++() { value += step; return *this; }
iotarator operator++(int) { auto ret = *this; ++*this; return ret; }

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

template<class Int>
class sieverator {
struct erased_iterator {
virtual Int dereference() = 0;
virtual void increment() = 0;
};
template<class It>
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<erased_iterator> m_ptr;

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

Int operator*() const {
return m_current;
}

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

sieverator& 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<(const pair& rhs) const {
return next_crossed_off_value < rhs.next_crossed_off_value ? true
: next_crossed_off_value > rhs.next_crossed_off_value ? false
: prime_increment < rhs.prime_increment;
}
bool operator>(const pair& rhs) const { return rhs < *this; }
};
template<class T> using min_heap = std::priority_queue<T, std::vector<T>, std::greater<>>;
min_heap<pair> m_pq;

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<__int128_t> iota(2);
sieverator<__int128_t> sieve(iota);
for (int p : sieve) {
std::cout << p << 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 <vector>
#include <iostream>
#include <locale>

unsigned long primes = 0;

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

for(int i = 2; i<number; i++) {
if(!sieve[i]) {
++primes;
for (int temp = 2*i; temp<number; temp += i)
sieve[temp] = true;
}
}
std::cout.imbue(std::locale(""));
std::cout << "found: " << primes << " 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<bool>. 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).

# Simplify the interface

I'll start by looking at how the code gets used in main(). We have to construct an iterator in order to create a sieverator; we'll most likely want to use an iotarator. I think that it would be more user-friendly if we could ask for a sieverator with a default iotarator. That can be achieved with a pair factory functions, to allow us to construct from either an iterator or from a (defaultable) integer value:

#include <iterator>
#include <type_traits>

template <typename T, typename V = typename std::iterator_traits<T>::value_type>
sieverator<V>
make_sieve(T value)
{
return sieverator<V>{value};
}

template <typename T = unsigned long long>
typename std::enable_if<std::is_integral<T>::value, sieverator<T>>::type
make_sieve(T value = 2)
{
return make_sieve(iotarator<T>{value});
}


For the iterator_traits to work, we must derive iotarator from public std::iterator<std::input_iterator_tag, Int>. With these changes, main() becomes

int main()
{
auto sieve = make_sieve<__uint128_t>();
for (int p : sieve) {
std::cout << p << std::endl;
}
}


# Consider whether the type erasure is worth the benefit

We only use one kind of input to the sieverator in this program, so we don't save on instantiations by erasing the type.

If we don't erase the type, the sieverator becomes much simpler:

template<class Iter, class Int = typename Iter::value_type>
class sieverator {
Iter m_ptr;

explicit sieverator() : m_ptr{} {}  // used by .end()
public:
template<class It>
explicit sieverator(It it) :
m_ptr(std::move(it))
{}
sieverator begin() { return std::move(*this); }
sieverator end() const { return sieverator{}; }
bool operator==(const sieverator&) const { return false; }
bool operator!=(const sieverator&) const { return true; }

Int operator*() const {
return *m_ptr;
}

sieverator& operator++() {
cross_off_multiples_of_prime(*m_ptr);
return *this;
}

sieverator& operator++(int) = delete;

//... unchanged after here


We have to change our helper functions to match the new type signature:

template <typename I, typename V = typename std::iterator_traits<I>::value_type>
sieverator<I,V>
make_sieve(I iter)
{
return sieverator<I,V>{iter};
}

template <typename T = unsigned long long, typename I = typename std::enable_if<std::is_integral<T>::value, iotarator<T>>::type>
sieverator<I, T>
make_sieve(T value = 2)
{
return make_sieve(I{value});
}


What happens when we use char as the value type? It's permitted, but we'll quickly end up with nonsense output. With a sufficiently powerful machine, we could (conceptually) overflow an unsigned 128-bit number, and it may be valuable to think about how to handle it (even if we'll never get there within our lifetimes!)

# Use std::tie to implement comparisons

The logic for comparing pair objects is a fairly standard "compare first element; if equal, compare next, and so on", which is implemented by std::tuple. So we can replace the operator< with

        return std::tie(next_crossed_off_value, prime_increment)
< std::tie(rhs.next_crossed_off_value, rhs.prime_increment);


Alternatively, just use std::pair, and live with the less descriptive member names first and second.

# Non-const begin() method

This makes me uncomfortable:

sieverator begin() { return std::move(*this); }


A begin() method is usually const, but this goes so far as to leave this in a moved-from state. I can see why it does this, but I would like to see a big fat warning comment!

# One-off type definition

There's a definition of min_heap that's used only once:

template<class T> using min_heap = std::priority_queue<T, std::vector<T>, std::greater<>>;
min_heap<pair> m_pq;


Whilst this serves to give a name to the type, it does seem like overkill here. We could make it a non-template:

using min_heap = std::priority_queue<pair, std::vector<pair>, std::greater<>>;
min_heap m_pq;


or just use a comment:

// m_pq is a min-heap
std::priority_queue<pair, std::vector<pair>, std::greater<>> m_pq;


# Naming

I don't like pair as a type name. It's uninformative, but more importantly, looks too much like std::pair.

# An observation

With GCC 6.3.0, this program feeds the optimizer well - to generate one million primes (as in the question) takes over 50 seconds with -O0; that reduces to 4 seconds at -O1 and under 3½ seconds at -O3!

• Re your first two sections: I did want to preserve the ability to use a sieverator with any iterator type. But you're certainly right that I should have done sieverator<It> instead of sieverator<decltype(*It{})>. That would have eliminated the type erasure and also the repetition of __int128_t! – Quuxplusone May 17 '17 at 14:55
• I should have been clearer - the helper functions come as a pair, so you can use any iterator you want; it just defaults the common case. Perhaps my edit makes this clearer? – Toby Speight May 17 '17 at 14:57
• Re "think about overflow": I thought about it. ;) I couldn't figure out anything really clever to do about it, especially since someone might very plausibly want to use a bignum as the input. – Quuxplusone May 17 '17 at 15:00
• Locally, I removed the type erasure and gained a bit more performance (maybe 5%). I'd be curious whether GCC does even better with it after the removal. :) (Removing the type erasure is pretty mechanical. I didn't modify the code in the question, per general CodeReview policy.) – Quuxplusone May 17 '17 at 21:03
• I should have been clear - my timings were on my modified code, not the original. I'd already started playing with it before I realised that -O3 would speed my edit-compile-test loop! – Toby Speight May 17 '17 at 21:13

Following your interest and comparison with the classic algorithm reminded by Jerry Coffin, I'd like to address your concern about the fact that in this form, it can't go to unplanned heights.

It is actually quite easy to modify the code to make that possible, which, considering the difference of efficiency we have, render your complicated implementation that uses a priority queue mostly useless. If you read the paper carefully, you will easily find that there is no problem in the classic imperative implementation that mutates a buffer of bools.

So here we go:

#include <cstdio>
#include <cstdint>
#include <vector>
#include <chrono>

class SieveOfEratosthenes
{
static constexpr int chunksize = 100'000;
public:
SieveOfEratosthenes() : m_compo(chunksize, false), m_next(2) {
for (std::int64_t i = 2; i < chunksize; i++)
if (!m_compo[i])
for (std::int64_t c = i * i; c < chunksize; c += i)
m_compo[c] = true;
}
std::int64_t next() {
const auto result = m_next++;
do {
for (; m_next < m_compo.size() && m_compo[m_next]; ++m_next);
} while (m_next >= m_compo.size() && extend());
return result;
}
private:
bool extend() {
const auto start = m_compo.size();
const auto end = start + chunksize;
m_compo.resize(end);
for (std::int64_t i = 2; i * i < end; i++)
if (!m_compo[i])
for (std::int64_t c = start - (start % i); c < end; c += i)
m_compo[c] = true;
return true;
}

std::vector<bool> m_compo;
std::int64_t m_next;
};

int main()
{

std::int64_t count = 0;
SieveOfEratosthenes sieve;
while (sieve.next() < 100'000'000)
count++;

std::printf("%lli\n", count);
std::printf("%d\n", (int)std::chrono::duration_cast<std::chrono::milliseconds>(t2 - t1).count());
}


I did not output the values found: I think this is not a good idea because this makes us going from the subject of discussing the sieve of Eratosthenes algorithm to the subject of benchmarking it at the same time as benchmarking IO of arbitrary systems. We should concentrate on the algorithm itself.

I expect this algo to be very roughly ten times faster than your "optimized" one, both for small and big values (~10'000'000). If a very small number of primes is needed, arguably the speed does not matter. If a very large number of primes is needed, some serious algorithmic optimizations will be needed anyway. I have not analyzed the complexity figures but I expect that they do not matter in practice except maybe for absurdly large numbers of primes, which like I explained consider out of scope for such toy algos.

Your code is nearly unreadable. There is so much clutter in the classes. Your sieverator has a public then a private then a public section, with a struct sprinkled in, there is a typedef with an immediate use... Just separate the code, it is not really gaining you too much here is it?

I am only 90% sure, but code like this

explicit iotarator(Int v)


denies construction by reference or move semantics. Although as you deal with basic types this will most likely not be a bottleneck

• move semantics are disabled when the members don't support it or copy semantics are deleted or not supported. – Incomputable May 16 '17 at 7:22
• Using references for fundamental types like int or int128 does usually not make sense – JVApen Apr 27 '18 at 6:28

Jerry Coffin's answer finally piqued my interest enough to look into it.

My first observation is that you can't meaningfully compare a program that does std::cout in its inner loop to a program that doesn't produce any output. :)

Adding std::cout << i << std::endl; to Jerry's inner loop (inside the if) slows it down from 0.197s to 3.520s, which is at least of the same general order of magnitude as my original program.

Jerry-with-no-I/O  0.197s
Jerry-with-cout    3.520s
Arthur-with-cout  11.995s


The next step is to eliminate the overhead of std::cout. Replacing std::cout << i << std::endl with printf("%d\n", i) reduces both of our running times by about 3 seconds:

Jerry-with-printf  0.528s
Arthur-with-printf 8.549s


The next step is to implement one's own priority_queue, implementing not just push and pop but also replace_top, since that operation is the bottleneck for this (and many other) algorithms.

Jerry                    0.528s
Arthur-with-replace_top  3.433s


Finally, I protest that we're comparing arbitrarily sized apples to fixed-size oranges here. Computing the 10-millionth prime (namely, 179'424'673) can be done with the "genuine sieve" program out of the box, but not with Jerry's — the latter requires that you recompile it with number > 179'424'673. And then to compute the 20-millionth prime, you have to recompile again. And so on.

Jerry-with-small-vector-computing-1Mth-prime            0.528s
Jerry-with-big-vector-computing-1Mth-prime              2.916s
Jerry-with-even-bigger-vector-computing-1Mth-prime      6.397s
Arthur-unmodified-computing-1Mth-prime                  3.433s


Jerry's definitely always runs (much!) faster at the "bleeding edge", but by the time you've altered it to compute the 20-millionth prime correctly, that same program is now running slower to compute the 1-millionth prime. (If you knew, up front, which prime you were going to be stopping at, you could choose a safe bound for the array in advance. But you can't do that if we're just piping the infinite stream of primes to head -n.)