Parallel factorial algorithm using std::thread

Question

This code calculates the factorial of a number on multiple threads. My issue: it is only a little bit faster than the sequential version of it (and I think I know why, I just can't find a way to solve this).

I use boost::multiprecision::cpp_int so the limits of default integers are not a problem, the size of integers is only limited by memory.

Only showing the relevant parts:

// ... other includes ...
#include <boost/multiprecision/cpp_int.hpp>

#define THREAD_COUNT 4
std::atomic<int> thread_num(1);         // global variable

// stuff...

void threaded_factorial(unsigned long long int num, boost::multiprecision::cpp_int& bigInt)
{
    int threadid = thread_num++;     // thread_num is atomic, so this is safe
    boost::multiprecision::cpp_int N = 1;
    for (unsigned long long int i = threadid; i <= num; i = i + THREAD_COUNT)
    {
        N *=(i);

    }

    std::lock_guard<std::mutex> lock(mu);      // race condition --> mutex needed
    bigInt *= N;
}

// more stuff ...

And the call of the function:

// ...

boost::multiprecision::cpp_int result = 1;
std::thread workers[THREAD_COUNT];

for (int i = 0; i < THREAD_COUNT; ++i)
{
    workers[i] = std::thread(threaded_factorial, num, std::ref(result));
}
for (int i = 0; i < THREAD_COUNT; ++i)
{
    workers[i].join();
}

// ...

The results seem correct, but as I said, this is not much faster than sequential code.

For example. The calculation of the factorial of 325253 took

• 67586 ms on 4 threads
• 76226 ms on a single thread

That is some really poor performance.

The reason, I think is that the for cycle in the threaded_factorial function roughly takes the same amount of time for each thread to complete, so when the std::mutex mu is locked, (THREAD_COUNT-1) threads have to wait for the one which locked the mutex.

This way, most of the work (by far the largest multiplications) is happening in a sequential manner, so the algorithm is really slow.

How can I work around this issue and make this work efficiently?

Answer 1

Firstly, some issues with API usage. Using raw std::threads is generally not the way you want to go with this sort of thing: prefer to use std::async. This also means you don't need to pass in by reference for updates, as it can return a value. The other big win from this is that you don't need to lock and perform an update in the thread that is running the calculation; this can be done independently.

Firstly, let's modify the threaded_factorial function:

constexpr static auto threads = 4U;
constexpr static auto test = 325253U;

namespace mp = boost::multiprecision;

mp::cpp_int thread_fact(unsigned num, int start)
{
    mp::cpp_int n = start;
    for (auto i = start + threads; i <= num; i += threads) {
        n *= i;
    }
    return n;
}

To call this, we setup some arrays for the std::futures that will be returned, as well as an array of partial results.

std::array<std::future<mp::cpp_int>, threads> futures;
std::array<mp::cpp_int, threads> results;
for (auto i = 1; i <= threads; ++i) {
    futures[i - 1] = std::async(std::launch::async, thread_fact, test, i);
}
for (auto i = 0; i < threads; ++i) {
    results[i] = futures[i].get();
}

Now, the step where you combine these is actually pretty expensive. Multiplying two numbers that are of this magnitude will be time consuming; let's launch the multiplications in separate threads as well:

std::future<mp::cpp_int> x = std::async([&results]() -> mp::cpp_int { return results[0] * results[1]; });
std::future<mp::cpp_int> y = std::async([&results]() -> mp::cpp_int { return results[2] * results[3]; });
auto x_val = x.get();
auto y_val = y.get();
auto z = x_val * y_val;

Making these changes, this runs in a bit under 10 seconds for me. In fact, from the profile graph, most of that time is spent doing the combining.

Others have pointed out that further algorithmic improvements are possible if you need more speed.

Answer 2

I would change this to a map/reduce problem.

1) Have a set of `N` mappers.
   Each mapper calculates the value for a range.
   Then saves the value for use by the reducer.

2) Reducer waits for all mappers to finish
   Then calculates a result based on the value generated
   by the mappers.

Using this technique you don't need any data locks. You just need a way to know when all the mappers have finished working.

// Sort of pseudo code.
int fact(int n) {

   int partitions = calculateNumberOfPartitions(n);
   int worker     = calculateNumberOfWorkers(n);

   int valuesPerPart = n+1 / partitions;
   if (partitions * valuesPerPart <= n) {
       ++valuesPerPart;
   }

   std::vector<boost::multiprecision::cpp_int>  data(partitions);
   boost::multiprecision::cpp_int               result;

   std::vector<std::function<void()>  jobs;

   // Calculate all factorial for all the partitions.
   for(int loop=0;loop < partitions; loop++) {
       jobs.push_back([&data, loop, n, valuesPerPart](){
             int low  = loop * valuesPerPart;
             int high = low  + valuesPerPart;
             high = high > n ? n+1 : high;

             boost::multiprecision::cpp_int  part = 1;
             for(int val = low; val < high; ++val) {
                 part *= val;
             }
             data[loop] = part;
       });
   }
   // The first (n-1) workers will finish
   // When they do force them to just wait for the last guy.
   std::vector<std::condition_variable>  wait(worker-1);
   for(int loop=0;loop < (worker-1); ++loop) {
       jobs.push_back([&wait, loop](){
           wait[loop].wait();
       });
   }
   // When the last worker finishes.
   // Let him do the reduce job.
   jobs.push_back([&data, &result](){
       for(auto& val: data) {
           result *= val;
       }
   });

   runJobsInParallel(jobs);

   // Now you can release the other workers you put to sleep.
}

Answer 3

Forget the threads for a while. I've tried it in Java and it took 21 seconds (no threading). All it needed was

com.google.common.math.BigIntegerMath.factorial(325253)

This tells us that some optimizations should be done, before you try to parallelize. The reason for the slowness is that you're working with increasingly bigger numbers step by step. Serially, you'd compute

 damnBigNumber *= 325251;
 damnBigNumber *= 325252;
 damnBigNumber *= 325253;

so there are many very costly operations there. Your threading should help here (no idea why it didn't), but reordering the multiplications would help surely a lot.

Imagine to compute two numbers a and b of similar size and finally multiplying them. This means that you have many operations with much shorter number and only the final multiplication needs to deal with the full length.

Another optimization is dealing with odd numbers only. Simply strip out all the factors of two and add them in the final step (it's left shift, a rather cheap operation).

---

With these optimizations, you'll get a different algorithm and its parallelization would differ a lot. As the time needed to compute a product is proportional to the product of operand sizes, the final multiplication time will probably dominate. So parallelizing this multiplication might offer the biggest gain. Maybe something like

using boost::multiprecision;
cpp_int mul(cpp_int a, cpp_int b) {
    // half of the number of bits of a
    // there's most probably a better method for this
    int split = int(log(a) / log(2) / 2);

    cpp_int al = lower_half(a, split);
    cpp_int ah = upper_half(a, split);
    cpp_int bl = lower_half(b, split);
    cpp_int bh = upper_half(b, split);

    // this part is to be parallelized
    cpp_int albl = al * bl;
    cpp_int albh = al * bh;
    cpp_int ahbl = ah * bl;
    cpp_int ahbh = ah * bh;

    return shiftLeft(ahbh, 2*split) + shiftLeft(albh+ahbl, split) + albl;
}

Answer 4

The OP's solution is breaks down at the end where it does NTHREADS multiplications of very large numbers sequentially (due to the mutex). This gets worse the more threads you use. It would be preferable to do them in parallel, but this would require more code with your solution. If you think of the whole problem as "reduce", it is easy to make it parallel all the way down. The trick is to do it without recursion.

Update for C++17: Instead of dealing with threads and locks yourself, with std::execution::par, you can use common std algorithms to solve the problem with a lot less code. The std algo implementations will (hopefully) be very optimized, possibly for the specifics of the platform.

Here's a "one line" solution that runs in 0.3s

using mp = boost::multiprecision::cpp_int;
constexpr static auto test = 325253U;

std::reduce(std::execution::par, boost::counting_iterator<long long>(1), boost::counting_iterator<long long>(test+1), mp(1), [](const mp a, const mp b) { return a*b;})};

because std algorithms use iterators, you might think you need to allocate and fill a container to pass in. The trick is to use a "fake" iterator such as boost::counting_iterator.

Here is the complete code:

#include <iostream>
#include <boost/multiprecision/cpp_int.hpp>
#include <boost/iterator/counting_iterator.hpp>
#include <numeric>
#include <execution>
#include <chrono>

using mp = boost::multiprecision::cpp_int;
constexpr static auto test = 325253U;

using namespace std;

int main(int argc, char * argv[]) {
    auto eval = [](auto fun) {
        const auto t1 = std::chrono::high_resolution_clock::now();
        const auto [name, result] = fun();
        const auto t2 = std::chrono::high_resolution_clock::now();
        const std::chrono::duration<double, std::milli> ms = t2 - t1;
        std::cout << std::fixed << std::setprecision(1) << name << " result "
                  << " took " << ms.count() << " ms\n";
    return result;
    };

    cout << "Running factorial(" << test << ")\n";

    const auto r1 = eval([]{ return std::pair{"std::reduce (par, mp)",
    std::reduce(std::execution::par,  boost::counting_iterator<long long>(1), boost::counting_iterator<long long>(test+1), mp(1), [](const mp a, const mp b) { return a*b;})}; });

    const auto r2 = eval([]{ return std::pair{"std::reduce (seq, mp)",
    std::reduce(std::execution::seq,  boost::counting_iterator<long long>(1), boost::counting_iterator<long long>(test+1), mp(1), [](const mp a, const mp b) { return a*b;})}; });

    cout << ( (r1==r2) ? "OK\n" : "FAIL\n" );
    
    return 0;
}

$ ./a.out
Running factorial(325253)
std::reduce (par, mp) result took 297.0 ms
std::reduce (seq, mp) result took 12763.6 ms
OK