Are numbers prime multithreading

Question

I'm trying to learn some multithreading and this seems like a good example, want to apply this to frustum culling later. The goal is to find prime numbers in randomly generated vector of numbers. Single threaded solution is simple, multithreaded way is to divide vector of random numbers in two parts, solve each one in the same time and merge results.

Here is the time it took to solve:

mt:

solved in: 17185.4 ms

solved in: 15832.2 ms

solved in: 16012.5 ms

solved in: 16354.9 ms

st:

solved in: 26942 ms

solved in: 30577.9 ms

solved in: 26403.4 ms

solved in: 29625.7 ms

It obviously works, but is the code correct?

This is from main.cpp:

 #include "primeNumbers.h"
 int main(){
    PrimeNumbers numbers(100000, 0, 5000000);
    PrimeNumberSolver* solver = new PrimeNumberSingleThread;
    //PrimeNumberSolver* solver = new PrimeNumberMultiThread;
    numbers.Solve(solver);
    delete solver;
 }

This is the rest of the code:

//primeNumbers.h
#include <mutex>
#include <thread>
#include <vector>
#include <iostream>
class PrimeNumberSolver abstract {
public:
    bool isPrime(int n) {
        for (int i = 2; i < n/2; ++i) {
            if (n % i == 0)
                return false;
        }

        return true;
    }

    virtual void solve(const std::vector<int>& source, std::vector<int>& result) = 0;

};

class PrimeNumberSingleThread final : public PrimeNumberSolver {
    friend class PrimeNumbers;
private:
    virtual void solve(const std::vector<int>& source, std::vector<int>& result) override{
        for (unsigned i = 0; i < source.size(); ++i) {
            if (isPrime(source[i]))
                result.push_back(source[i]);
        }
    }
};

class PrimeNumberMultiThread final : public PrimeNumberSolver {
    friend class PrimeNumbers;

private:
    virtual void solve(const std::vector<int>& source, std::vector<int>& result) override {
        std::vector<int> r1;
        std::vector<int> r2;

        //split source
        int half = source.size() / 2;
        int end = source.size();

        std::mutex m;

        //run two threads 
        std::thread t1([&]() {
            for (int i = 0; i < half; ++i) {
                //std::lock_guard<std::mutex> lock(m);  //do i need to lock mutex before accessing source
                if (isPrime(source[i]))
                    r1.push_back(source[i]);
            }
        });
        std::thread t2([&]() {
            for (int i = half; i < end; ++i) {
                //std::lock_guard<std::mutex> lock(m);  //do i need to lock mutex before accessing source
                if (isPrime(source[i]))
                    r2.push_back(source[i]);
            }
        });

        //join threads
        if(t1.joinable()) t1.join();
        if(t2.joinable()) t2.join();

        //merge results
        result.insert(result.end(), r1.begin(), r1.end());
        result.insert(result.end(), r2.begin(), r2.end());
    }
};

class PrimeNumbers {
public:
    PrimeNumbers(int size, unsigned min, unsigned max) {
        std::mt19937 rng;
        for (int i = 0; i < size; ++i) {
            rng.seed(std::random_device()());
            std::uniform_int_distribution<std::mt19937::result_type> dist6(min, max);
            randoms.push_back(dist6(rng));
        }
        std::cout << "Numbers are generated" << std::endl;
    }

    void Solve(PrimeNumberSolver* solver) {
        double time = glfwGetTime();
        solver->solve(randoms, primes);
        std::cout << "solved in: " << (glfwGetTime() - time) * 1000.0 << " ms" << std::endl;
    }

    void printPrimes() {
        for (unsigned i = 0; i < primes.size(); ++i)
            std::cout << primes[i] << std::endl;
    }

    std::vector<int> randoms;
    std::vector<int> primes;
};

Answer 1

This is a poor way to test primality:

bool isPrime(int n) {
    for (int i = 2; i < n/2; ++i) {
        if (n % i == 0)
            return false;
    }
    return true;
}

Even sticking with trial division, we can stop at √n, which is generally much less than n/2. It's also well worth treating 2 as a special case, to reduce the iterations by half (i.e. i += 2 instead of ++i).

Why are we using signed integers, given that we don't need any negative values?

You say the "output is tested", but don't show any of the test code. It's useful to see how the test works, to identify tests that are missing or erroneously passing. I mention this because it's notoriously difficult to create reliable tests for multithreaded code.

In PrimeNumberMultiThread::solve, we have opportunity to reduce duplication, and we don't need to hard-code the number of threads to two. Instead, we can use as many threads as we have cores available, using a for loop to start them and another for loop to join them all.

The mutex m and the commented-out lock_guard should just be removed: these aren't necessary, because only one thread accesses each vector between the worker threads starting and finishing (which points are memory barriers).

Instead of result.insert(), the general case might benefit from std::move algorithm - in any case, it's well worth a reserve before adding the results.

With OpenMP, much of the work can be automated for you - you'll just need a custom reduction to combine vectors - I think that's available from OpenMP 4.0 onwards.

An alternative to using vectors is to use an array of flags to indicate which values are to be included in the result (the flags could even be in the input vector: overwrite the non-primes with zero to discard them). Then there's no reallocation within either thread, and you can use OpenMP of older versions because we no longer need a custom reduction.

Here's one version that uses the OpenMP approach I've just described:

#include <algorithm>
#include <iostream>
#include <iterator>
#include <random>
#include <vector>

bool isPrime(unsigned n) {
        if (n < 2) return false;
        if (n < 4) return true;
        if (n%2 == 0) return false;
        for (auto i = 3ul;  i*i < n;  i += 2) {
            if (n % i == 0)
                return false;
        }

        return true;
}

std::vector<unsigned> solve(std::vector<int>& source)
{
    auto const len = source.size();
    std::size_t result_size = 0;

#pragma omp parallel for reduction(+:result_size)
    for (std::size_t i = 0;  i < len;  ++i) {
        if (isPrime(source[i])) {
            ++result_size;
        } else {
            source[i] = 0;
        }
    }

    std::vector<unsigned> result;
    result.reserve(result_size);

    std::copy_if(source.begin(), source.end(),
                 std::back_inserter(result),
                 [](unsigned i){ return i; });

    return result;
}


int main(){
    const std::size_t size = 1000000u;
    const unsigned min = 0u;
    const unsigned max = 5000000u;

    std::vector<int> randoms;
    randoms.reserve(size);

    {
        std::uniform_int_distribution<std::mt19937::result_type> dist(min, max);
        std::mt19937 rng;
        rng.seed(0);            // for a reproducible test
        for (std::size_t i = 0;  i < size;  ++i) {
            randoms.push_back(dist(rng));
        }
    }

    auto primes = solve(randoms);
    std::clog << "Found " << primes.size() << " primes "
              << " from " << randoms.size() << " inputs\n";
}

Compile with g++ -fopenmp, or the equivalent for your compiler.

To test a single-threaded version, we just have to remove (manually, or with #ifdef) the #pragma omp line.

---

Note that since C++17, we're able to specify an execution policy for std::copy_if, so we could reduce that code even further, especially if we're not fussy about the order of the results:

#include <execution>

std::vector<unsigned> solve(std::vector<int>& source)
{
    std::vector<unsigned> result;

    std::copy_if(std::execution::par_unseq
                 source.begin(), source.end(),
                 std::back_inserter(result),
                 [](unsigned i){ return isPrime(i); });

    return result;
}

I haven't been able to test this, as I'm still working with an old standard library that doesn't have <execution>.