Parallel sieve of Eratosthenes

Question

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 <thread> #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<std::thread> threads; 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); threads.emplace_back(apply_prime<Integer>, 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
– Loki Astari, Commented Apr 12, 2014 at 20:03
@LokiAstari No I didn't time it. Actually, I didn't care at all about speed, I just wanted to write a multithread program and get it reviewed for the sole purpose of learning. — Morwenn
– Morwenn, Commented Apr 12, 2014 at 20:06
My computer isn't even a multicore to start with. It would probably have troubles beating the cache-friendly sequential approach with threads. — Morwenn
– Morwenn, Commented Apr 12, 2014 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
– vzn, Commented Apr 13, 2014 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.
– Niklas B., Commented Apr 13, 2014 at 5:12

Nobody moving away from SE · Accepted Answer · 2014-04-12 19:49:19Z

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:

Thread safety

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
– Morwenn, Commented Apr 12, 2014 at 19:55
@BenVoigt Fortunately, std::vector<std::atomic<bool>> is fine. — Morwenn
– Morwenn, Commented Apr 13, 2014 at 20:46

Edward · Accepted Answer · 2014-04-13 01:10:48Z

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
– Edward, Commented Apr 12, 2014 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
– Morwenn, Commented Apr 12, 2014 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
– Morwenn, Commented Apr 12, 2014 at 22:00
You might want to read this paper for more ideas on how to parallel-ize this algorithm. — Edward
– Edward, Commented Apr 12, 2014 at 22:09
That said, even if the results are now consistent, I am pretty sure that useless threads are still spawned. However, they now are only useless and not harmful. — Morwenn
– Morwenn, Commented Apr 12, 2014 at 22:27

Loki Astari · Accepted Answer · 2014-04-12 20:44:38Z

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).

I usually prefer n+1 where n is the number of processors. Possibly n+2 in IO heavy cases. — Emily L.
– Emily L., Commented Jul 15, 2015 at 14:26

gnasher729 · Accepted Answer · 2014-04-13 08:46:56Z

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
– Morwenn, Commented Apr 13, 2014 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
– Ben Voigt, Commented Apr 13, 2014 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
– Morwenn, Commented Apr 13, 2014 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
– gnasher729, Commented Apr 13, 2014 at 23:24
@gnasher729 I did even more research since yesterday, and you're totally right. — Morwenn
– Morwenn, Commented Apr 14, 2014 at 11:41

Morwenn · Accepted Answer · 2014-04-13 10:30:15Z

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 moving away from SE
– Nobody moving away from SE, Commented Apr 13, 2014 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
– Morwenn, Commented Apr 13, 2014 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
– Edward, Commented Apr 13, 2014 at 15:45
@Edward But that also means, that I could use i += 2u*prime for everywhere prime beside 2. — Morwenn
– Morwenn, Commented Apr 13, 2014 at 16:13

Stack Exchange Network

Parallel sieve of Eratosthenes

5 Answers 5

Naming

Efficiency

Thread safety

You must log in to answer this question.

Linked

Hot Network Questions

Parallel sieve of Eratosthenes

5 Answers 5

Naming

Efficiency

Thread safety

You must log in to answer this question.

Linked

Related

Hot Network Questions