# Parallel sieve of Eratosthenes

For the first time, I tried to use threads by myself in order to implement a parallel sieve of Eratosthenes. The trick is that each time a prime is found, a thread is spawned to eliminate all the multiples of this prime number from the boolean vector (the one that tells whether a number is a prime or not). Here is my code:

#include <cmath>
#include <functional>
#include <vector>

// Invalidates the multiples of the given integer
// in a given boolen vector
template<typename Integer>
void apply_prime(Integer prime, std::vector<bool>& vec)
{
for (Integer i = prime*2u ; i < vec.size() ; i += prime)
{
vec[i] = false;
}
}

template<typename Integer>
auto sieve_eratosthenes(Integer n)
-> std::vector<Integer>
{
std::vector<bool> is_prime(n, true);
std::vector<Integer> res;

auto end = static_cast<Integer>(std::sqrt(n));
for (Integer i = 2u ; i <= end ; ++i)
{
// When a prime is found,
// * add it to the res vector
// * spawn a thread to invalidate multiples
if (is_prime[i])
{
res.push_back(i);
i, std::ref(is_prime));
}
}

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

// Add the remaining primes to the res vector
for (Integer i = end+1u ; i < is_prime.size() ; ++i)
{
if (is_prime[i])
{
res.push_back(i);
}
}
return res;
}


The primes are added in two steps to the res vector: every prime $p$ such as $p < \sqrt{n}$ is added when the prime is found, before the corresponding thread is thrown. The other primes are added at the end the of the function. Here is an example main:

int main()
{
auto primes = sieve_eratosthenes(1000u);

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


I was pretty sure that I would get some problems due to parallelism, but for some reason, it seems to work. I got the expected results in the right order. Just to be sure, I would like to know whether my program is or correct or whether it has some threading issues that I couldn't see.

Note: I used many of the ideas from the answer to improve the code and wrote a follow-up question.

-
Have you timed this? It would surprise me if a parallel sieve is faster because of the memory contention. The speed up you receive from cached memory over real memory seem more likely to give you a speed up (as you can't do as much local caching with threads as different threads may be on different cores and you need to keep pushing things backwards and forwards across caches and memory). –  Loki Astari Apr 12 '14 at 20:03
My computer isn't even a multicore to start with. It would probably have troubles beating the cache-friendly sequential approach with threads. –  Morwenn Apr 12 '14 at 20:12
+1; parallelism in prime searching is a very interesting advanced idea & wonder if there is more scientific analysis of this somewhere...refs anyone? anyway note that the science of prime detection in general is highly advanced and theoretical and sieve of eratosthenes while respectable as a programming exercise is regarded by experts as basically a "toy" algorithm for the problem... my understanding GNFS is the leading/typical algorithm, wonder if it has been parallelized by anyone? –  vzn Apr 13 '14 at 0:41
@LokiAstari Actually if you precompute the primes up to sqrt(maxp), you can then partition the sieving space evenly between the processors and get no contention at all. The same trick works to make the sequential sieve algorithm much more cache-efficient. –  Niklas B. Apr 13 '14 at 5:12
@NiklasB.: Its interesting that you think that one sentence written by you is the same as a whole page of detailed algorithms by somebody else? :-) –  Loki Astari Apr 13 '14 at 21:08

As you are correctly assuming there are multiple threading related issues with your code, but lets tease you and start with the usual suspects.

# Naming

The name apply_prime is misleading and inexpressive. Neither does the function really require a prime nor does apply do it justice. You should name it something along the lines: strike_out_multiples or something similar.

# Efficiency

Your loop can be twice as fast by doing steps of two instead of one.

And now for the main act:

I am no expert here but I can spot at least two problems:

Ordering:

Nothing prevents the system from "favoring" your main thread and letting all others run after the main thread is hitting the join loop. This results in wrong results as the mainthread plowed through the completely true vector before anyone could tell it that for example 4 is not prime.

Parallel writes:

This might be a bit of a corner case and not a problem anymore but at least in the old standard std::vector<bool> was a specialization that used only one bit per value. While this saves you some space it comes with the cost of actually accessing 8 bits when only working on one. This means that two threads might well be working on different bits but in the same byte. Consider this:

Thread 1 writes false to bit 4 while thread 2 writes false to bit 6. Both are located in byte 0. Now both start out with reading the initial value of the byte, say (true, true, true, false, true, true, true, true) and stores them! Now one of the two threads will write its value later than the other and overwrites the false of the other one -> lost update.

Maybe this can even happen with byte sized booleans as some architectures do not allow for only byte-wise access but only word wise (which results in the same problem, only with wider sizes).

-
I wanted to make apply_prime a lambda, but thought multi-lines lambda were not allowed. Actually, they are allowed and a lambda is fine. –  Morwenn Apr 12 '14 at 19:55
std::vector<bool> is still a problem –  Ben Voigt Apr 13 '14 at 15:11
@BenVoigt Fortunately, std::vector<std::atomic<bool>> is fine. –  Morwenn Apr 13 '14 at 20:46

Generally it's nice, well-structured code, but it relies on a promise that might not be kept.

Specifically, simultaneous access to different elements in std::vector<bool> is not guaranteed to be thread-safe because storage bytes may be shared by multiple bits in the vector.

Consider an alternative way to slice things. Each thread could be responsible for its own section of the boolean array. As primes become known, they could be dispatched to each of the threads for simultaneous elimination from the corresponding section.

It might also be nice to have a tuning parameter in which the the size of the subarray is balanced against the cost of spawning another thread.

Edit: I modified main as follows:

int main()
{
auto primes = sieve_eratosthenes(20000000u);

long long s=0;
for (auto prime: primes)
s += prime;
std::cout << s << '\t' << primes.size() << std::endl;
}


I then ran the code four times and got this result:

12273796368896  1270814
12273126258541  1270843
12273106282821  1270780
12272824476679  1270794


So either the number of prime numbers is actually changing from iteration to iteration (which mathematicians generally consider to be an unlikely possibility!) or the thread contention issue is manifesting itself. The correct numbers are:

12272577818052  1270607


Edit 2: I ran across this paper which nicely explains a number of possible approaches to parallelization of the Sieve of Eratosthenes algorithm.

-
I also just attempted to have it calculate all of the primes up to 200000000 and it aborted on my machine after it ran out of threads, throwing an instance of std::system_error with what() = "Resource temporarily unavailable." –  Edward Apr 12 '14 at 21:32
It probably boils down to this problem: threads are spawn on supposed prime numbers. That means that depending on the order of the threads, 4 may be considered prime, be added to res and spawn a thread. While having threads spawned for non-prime numbers is not a problem (4 assigns false to places where 2 would also have assigned false), the push_back to res is probably what adds the wrong elements to the table. –  Morwenn Apr 12 '14 at 21:59
Could you try to remove the push_back and have the last loop go from 2u to is_prime.size() and check again? Also, replace the vector<bool> by a vector<unsigned>, it should be better :) –  Morwenn Apr 12 '14 at 22:00
Because I'm parsimonious, I made it a vector<char> but yes, it works consistently now. –  Edward Apr 12 '14 at 22:07
You might want to read this paper for more ideas on how to parallel-ize this algorithm. –  Edward Apr 12 '14 at 22:09

Ignoring the memory problems with cache invalidation (which will slow the code down).

Creating a thread is relatively expensive (as you have to allocate a stack and maintain it). So rather than creating and destroying threads it is better to maintain a thread pool and reuse the threads.

The number of threads to put in the pool should be slightly larger than the total number of cores you have (as a rule of thumb * 1.x (where x is in the range 2=>5) is the number of threads you should put in the pool).

-
It seems that I totally have to learn about thread pools :) –  Morwenn Apr 12 '14 at 21:53

Actually, the threading is just totally wrong. There are race conditions all over the place. Due to the kind of problem, it doesn't have many visible effects, but it is still wrong.

The out loop in sieve_eratosthenes loops over indices from 2 to end and checks whether array elements are marked as "prime". In that loop it starts threads which change array elements from "prime" to "non-prime". There is no guarantee how much progress these threads have been making. So when the outer loop checks if 4, 6, 8, 9, 10, 12 etc. are primes, there is no guarantee that they actually have been marked as non-prime. Worst case, the outer loop starts a thread for each i from 2 to end before any of these threads are running. In the case of a prime sieve that doesn't change the correctness, but in general this program doesn't do what it is supposed to do.

The threads set boolean values to false. They tend to set the same values to false multiple times. For example the threads for 2 and 3 both set each multiple of 6 to false. Two threads writing to the same variable is known as a "race condition" and is more than just bad. In this case, it doesn't have an effect because both threads set the same boolean value to false. In general, it will cause serious bugs.

Most people looking seriously at primes will use large numbers and use bit arrays to safe space. When you do that, this threading problem will kill you, because if one thread tries to clear one bit in a memory word while another thread tries to clear another bit, you can be quite confident that one of these operations will be lost.

-
I considered the race condition when writing the code (multiple threads wirting at the same place), but since the threads are known to perform exactly the same operation at the same place, I considered that it was not a problem (I hope that they just both write false and don't trigger undefined behaviour). All the other concerns are real problems though. –  Morwenn Apr 13 '14 at 8:59
@Morwenn: There's a write->read race as well as redundant writes, the read occurs in if (is_prime[i]) –  Ben Voigt Apr 13 '14 at 15:11
I read a bit more about data races today and I indeed have an undefined behaviour. I don't think that redundant writes are a problem though: they would probably not be much slower than a read then a write if the the type was atomic. –  Morwenn Apr 13 '14 at 15:14
Undefined behaviour is undefined behaviour. It means you absolutely cannot predict what happens. It may work on your current compiler. It may work on your current compiler 99% of the time, which is worse than not working at all. With multithreading, the approach "well, it's not officially working, but it's good enough for me" is an absolute recipe for desaster. Anything less than "it's guaranteed to work" will come back and bit you. –  gnasher729 Apr 13 '14 at 23:24
@gnasher729 I did even more research since yesterday, and you're totally right. –  Morwenn Apr 14 '14 at 11:41

None of the other answers mentioned that there is actually a problem when the value passed to the sieve is a prime number. For example, sieve_eratosthenes(7u) returns a vector containing 2 3 5. This is due to the vector is_prime being too short by one element: it considers the elements between $0$ and $n-1$ while it should consider the elements between $0$ and $n$. It should have been declared as:

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


Also, as mentioned in one of the comments, I can replace apply_prime by a lambda. Taking is_prime by reference in the lambda capture also allows to drop std::ref and the corresponding #include <functional>. Here is the modified threads.emplace_back:

threads.emplace_back([&is_prime](Integer prime)
{
for (Integer i = prime*2u ; i < is_prime.size() ; i += prime)
{
is_prime[i] = false;
}
}, i);

-
To be fair: You did not specify whether the input number should be tested or not, which indicates a lack of specification (preconditions). –  Nobody Apr 13 '14 at 13:43
@Nobody That's true, but I have the feeling that mathematical integer bounds generally tend to be inclusive. Also the pseudocode not Wikipedia says "not exceeding n". It seems to be inclusive. –  Morwenn Apr 13 '14 at 13:49
You can also replace i=prime*2u with i=prime*prime because all composite numbers less than $p^2$ will be eliminated as multiples of smaller primes. –  Edward Apr 13 '14 at 15:45
@Edward But that also means, that I could use i += 2u*prime for everywhere prime beside 2. –  Morwenn Apr 13 '14 at 16:13
I actually do way too many operations on even numbers that could be avoided... –  Morwenn Apr 13 '14 at 19:47