Parallel sieve of Eratosthenes, version 2

Question

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);
    }
    std::vector<std::thread> threads;

    // 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
    std::size_t nb_primes_per_thread =
        static_cast<std::size_t>(std::ceil(
            static_cast<float>(primes.size()) /
            static_cast<float>(std::thread::hardware_concurrency())
        ));

    for (std::size_t first = 0u;
         first < primes.size();
         first += nb_primes_per_thread)
    {
        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
        threads.emplace_back(
            [&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)
    {
        if (is_prime[i].load(std::memory_order_relaxed))
        {
            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.

Answer 1

Expanding on comments:

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.
        static_cast<float>(std::thread::hardware_concurrency())

for (std::size_t first = 0u;
     first < n;
     first += pageSizePerThread)
{
    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
    threads.emplace_back(
        [&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).

Answer 2

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.

Answer 3

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.