Parallel factorial algorithm using std::thread
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 &lt;= 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&lt;std::future&lt;mp::cpp_int&gt;, threads&gt; futures;
std::array&lt;mp::cpp_int, threads&gt; results;
for (auto i = 1; i &lt;= threads; ++i) {
    futures[i - 1] = std::async(std::launch::async, thread_fact, test, i);
}
for (auto i = 0; i &lt; 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&lt;mp::cpp_int&gt; x = std::async([&amp;results]() -&gt; mp::cpp_int { return results[0] * results[1]; });
std::future&lt;mp::cpp_int&gt; y = std::async([&amp;results]() -&gt; 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.