# Parallel sieve of Eratosthenes, version 2

This question is a revision of Parallel sieve of Eratosthenes. The goal is to implement a sieve of Eratosthenes with parallel strikes out from the boolean array. I tried to fix the data races and all the threading-related errors as well as to add some of the ideas from the previous answers. Now, the implementation works as follows:

• Compute the prime numbers $p$ such as $p <= \sqrt{n}$ thanks to a sequential sieve of Eratosthenes.
• Initialize a std::vector<std::atomic<bool>> with true for indices between $0$ and $n$ (inclusive).
• Compute the number of threads that should be used to concurrently strike out values from the vector. This number depends on the maximum number of concurrent threads allowed by the implementation and on the number of precomputed prime numbers.
• Spawn threads that will strike out the multiples of a given set of prime numbers from the boolean vector.
• Join the threads.
• Find the actual remaining prime numbers.

To differentiate between the sequential and parallel versions of the sieve, I had the function take an execution policy parameter first, inspired from the C++ parallelism TS N3960. The execution policy classes can be trivially implemented as such (I know... we shouldn't add anything new to std::):

namespace std
{
namespace parallel
{
struct sequential_execution_policy {};
struct parallel_execution_policy {};
struct vector_execution_policy {};

constexpr sequential_execution_policy   seq = sequential_execution_policy();
constexpr parallel_execution_policy     par = parallel_execution_policy();
constexpr vector_execution_policy       vec = vector_execution_policy();
}}


Here is the sequential version of sieve_eratosthenes. In the end, I did not make strike_out_multiples a lambda since I use it multiple times.

template<typename Integer, typename T>
void strike_out_multiples(Integer n, std::vector<T>& vec)
{
for (Integer i = n*2u ; i < vec.size() ; i += n)
{
vec[i] = false;
}
}

template<typename Integer>
auto sieve_eratosthenes(std::parallel::sequential_execution_policy, Integer n)
-> std::vector<Integer>
{
if (n < 2u)
{
return {};
}

std::vector<char> is_prime(n+1u, true);

// Strike out the multiples of 2 so that
// the following loop can be faster
strike_out_multiples(2u, is_prime);

// Strike out the multiples of the prime
// number between 3 and end
auto end = static_cast<Integer>(std::sqrt(n));
for (Integer n = 3u ; n <= end ; n += 2u)
{
if (is_prime[n])
{
strike_out_multiples(n, is_prime);
}
}

std::vector<Integer> res = { 2u };
for (Integer i = 3u ; i < is_prime.size() ; i += 2u)
{
if (is_prime[i])
{
res.push_back(i);
}
}
return res;
}


And now, here is the new parallel version of sieve_eratosthenes:

template<typename Integer>
auto sieve_eratosthenes(std::parallel::parallel_execution_policy, Integer n)
-> std::vector<Integer>
{
if (n < 2u)
{
return {};
}

// Only the prime numbers <= sqrt(n) are
// needed to find the other ones
auto end = static_cast<Integer>(std::sqrt(n));
// Find the primes numbers <= sqrt(n) thanks
// to a sequential sieve of Eratosthenes
const auto primes = sieve_eratosthenes(std::parallel::seq, end);

std::vector<std::atomic<bool>> is_prime(n+1u);
for (auto i = 0u ; i < n+1u ; ++i)
{
is_prime[i].store(true, std::memory_order_relaxed);
}

// Computes the number of primes numbers that will
// be handled by each thread. This number depends on
// the maximum number of concurrent threads allowed
// by the implementation and on the total number of
// elements in primes
static_cast<std::size_t>(std::ceil(
static_cast<float>(primes.size()) /
));

for (std::size_t first = 0u;
first < primes.size();
{
auto last = std::min(first+nb_primes_per_thread, primes.size());
// Spawn a thread to strike out the multiples
// of the prime numbers corresponding to the
// elements of primes between first and last
[&primes, &is_prime](Integer begin, Integer end)
{
for (std::size_t i = begin ; i < end ; ++i)
{
auto prime = primes[i];
for (Integer n = prime*2u ; n < is_prime.size() ; n += prime)
{
is_prime[n].store(false, std::memory_order_relaxed);
}
}
},
first, last);
}

for (auto& thr: threads)
{
thr.join();
}

std::vector<Integer> res = { 2u };
for (Integer i = 3u ; i < is_prime.size() ; i += 2u)
{
{
res.push_back(i);
}
}
return res;
}


And of course, an example main:

int main()
{
auto primes = sieve_eratosthenes(std::parallel::par, 10000u);

for (auto prime: primes)
{
std::cout << prime << " ";
}
}


Is this code correct (I think that I have resolved the threading-related issues), and if so, how can I make it even better? I can think of two ideas from the previous question that I did not implement:

• Using a thread pool (still have to learn how to create one).
• Using some heuristic to decide when using multiple threads may be better than using only one thread. That heuristic would probably rely on the number of prime numbers handled by each thread.

Note: the code does not work with clang++ 3.5: it seems that libc++ std::vector<std::atomic<T>>(std::size_t count) somehow tries to copy some std::atomic<T> instances while it should not copy anything.

• Note: The constructor you mention is the fill constructor. vector(size_type n, const value_type& val = value_type()). The second argument (if not supplied) is defaulted; but it is also "copied" into every cell on the vector. It may work better if you use resize() as that value initialize the extra elements. – Martin York Apr 13 '14 at 20:13
• @LokiAstari This constructor has been split into two constructors when we passed to C++11, and the one with only count guarantees that no copies are made. – Morwenn Apr 13 '14 at 20:24
• Looks like it should work. But you have all the threads accessing all the portions of memory. So behind the scenes there will be lots of page locking to that memory is correct. You current algorithm splits the read memory prime across threads. But each thread must also access all write memory is_prime. I woud reverse that. There is no contention if all threads use all parts of read memory prime then divide the write memory is_prime across the threads. So each thread is writing to its own personal chunk of write memory. – Martin York Apr 13 '14 at 20:25
• @LokiAstari I will have to see what I can do about it. Considering how the algorithm works, I have troubles seeing how I can split the write memory between threads. I understand what you mean, but implmenting it does not seem easy. – Morwenn Apr 13 '14 at 20:30
• Just a quick note, std::thread::hardware_concurrency returns 0 if it cannot compute a value, which would cause a SIGFPE in nb_primes_per_thread. – Yuushi Apr 14 '14 at 3:17

I would reverse the splitting. I would split the writable memory across the threads so each thread wrote to a unique section of memory (it would be good if this were divided by page size as it may help performance). We will call the section of memory that a thread writes to a page.

Then each thread will loop over all the primes and remove them from the page associated with the thread. Using the original code as a starting point it would look like this:

std::size_t pageSizePerThread =
static_cast<std::size_t>(std::ceil(
static_cast<float>(n) /        // Size of write memory.

for (std::size_t first = 0u;
first < n;
{
auto last = std::min(first+pageSizePerThread, n);
// Spawn a thread to strike out the multiples
// of the prime numbers corresponding to the
// elements of primes between first and last
[&primes, &is_prime](Integer begin, Integer end)
{
for (auto prime: primes)
{
Integer first = begin/prime;
first += (first < begin) ? prime : 0;
for (Integer n = first ; n < end ; n += prime)
{
is_prime[n].store(false, std::memory_order_relaxed);
}
}
},
first, last);
}


Note: It does not matter if read memory is_prime is shared across all the threads as each thread can have its own copy in the local cache and not be affected by other threads (as it is read only).

• (begin-1)/prime + 1 would remove a branch. And you're missing a multiply by prime. – Ben Voigt Apr 13 '14 at 21:21
• @BenVoigt: Can you expand on that. I am sure you are correct but cant quite place it in the code. – Martin York Apr 13 '14 at 21:23
• Integer first = std::max((begin-1)/prime + 1, prime)*prime; – Ben Voigt Apr 13 '14 at 21:24
• I presume that in your original code, you meant for first = begin/prime*prime ? – Ben Voigt Apr 13 '14 at 21:29
• @BenVoigt: Yes. But integer rounding means that was not enough. Hence the extra conditional. I assume your code works. But have not thought it through. – Martin York Apr 13 '14 at 21:36

Equally-sized sequential subsets of primes don't evenly divide the task. The smallest numbers lead to a lot more multiples.

Test-before-set might significantly reduce cache contention. Even better, don't start at prime*2u, start at prime*prime.

• The first remark is really valuable. While it now seems obvious, it didn't strike me at all when I wrote the code. – Morwenn Apr 13 '14 at 21:51

It has already been mentioned in an answer to the previous question, but there are really many places where the code does useless iterations on multiples of 2. Actually, these useless operations can be removed almost everywhere:

• In the sequential version of the sieve, strike_out_multiples(2u, is_prime); can be removed since the multiples of 2 are not even considered by the algorithm in the last part of the function where the actual prime numbers are added to the vector. Therefore, we don't need strike_out_multiples as a separate function anymore since it's only used once.
• Similarly, first can be declared as std::size_t first = 1u in the parallel version of the algorithm: 0u corresponds to the index of the vector where the prime number 2u is stored.
• There are loops that begin with 2u*prime at several places in the code. It has been proposed in another answer that these loops should begin at prime*prime instead, which is always odd since we ignore 2u. Consequently, the increment of these loops can be changed to n += 2u*prime (instead of n += prime) to ignore the even values.
• The initialization of is_prime in the parallel version of the sieve can begin at the index 3u (the previous ones can be ignored) and have an increment of 2u to skip the even values.

This highlights the fact that the vector is_prime contains twice as many values as it could contain: the values below 3u and all the even values are not used. There ought to be a way to avoid uselessly storing that many values.