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.

Answer 1

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

Answer 2

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.

Answer 3

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

Answer 4

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.

Answer 5

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