# Parallel factorial algorithm using std::thread

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>

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;

for (int i = 0; i < THREAD_COUNT; ++i)
{
}
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?

• The cost of spinning up a thread is non trivial. Usually you spin up a bunch of threads then give them lots of small pieces of work (only spin them down when things are finished). You are spinning your threads up on a task bases (thats a bad idea). – Martin York Jun 24 '15 at 20:52
• Before I tried boost::multiprecision::cpp_int, I coded this with std::atomic<int>. That way I got about 2.6 times faster execution times with 4+ threads, the only problem was that the results were not correct as the default integer overflowed. Is there an atomic implementation of arbitrary precision integers? Or is it even possible to implement such integers? (might be a stupid question, as I have no idea how atomic types work, I only know that they are thread-safe) – krispet krispet Jun 24 '15 at 20:59
• You'll be much better off just using a better algorithm. See luschny.de/math/factorial/FastFactorialFunctions.htm. The "If you do not attach great importance to high performance" section gives a naïve implementation which should still significantly outperform this. – Veedrac Jun 24 '15 at 21:24
• Since factorials are known. You should not be recalculating them. The result is fact(int n){return factData[n];} can't get much faster than that. This is not a really good problem to practice parallelism on in my opinion. – Martin York Jun 24 '15 at 21:44

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

• Exactly the kind of answer I was looking for! Thanks. Also, if std::async is preferred in these kind of problems, when should I use std::thread? – krispet krispet Jun 25 '15 at 13:28
• @krispetkrispet Generally when you want to launch something and forget about it - like running a thread that will listen for network traffic etc. I suppose a rule of thumb is if you except it to return a result, you should probably use std::async. – Yuushi Jun 25 '15 at 21:54

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

• @Veedrac: I am not sure it is better either. I just offer it as an alternative worth investigating. I would bet there is a lot of stalling happening as the threads go for the lock (I could be totally wrong). But this avoids processor stalls. – Martin York Jun 24 '15 at 22:18
• @Veedrac: Then its worth timing to find out. – Martin York Jun 24 '15 at 22:23
• @Veedrac: The optimism of youth. But its not the time for the lock that's important. Its the time for the processor stalls which can be orders of magnitude greater than the lock. – Martin York Jun 24 '15 at 22:27
• @Veedrac: In this version no processes is stalled waiting to do work. Which they do on the original (as they can potentially all stall waiting for access to the shared resource). They basically stop when there is no work (so yes they are way under utilized during the reduce phase). But at no cost to the speed. But this is something that is easily resolved by timming. There is no point in arguing about it. Just time it. I have no idea if it is faster or not. – Martin York Jun 24 '15 at 22:56
• Ack, I was being silly. Nvm.I thought you were talking about stalls to do the reduction. – Veedrac Jun 25 '15 at 4:26

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