Длинная арифметика

#include <iostream>
#include <fstream>
#include <iterator>
#include <vector>
#include <list>
#include <time.h>
#include <boost/multiprecision/cpp_int.hpp>
 
using namespace boost::multiprecision;
 
 
std::vector<int> primes(size_t n) {
    ++n;
    if (n < 3) return {};
 
    std::vector<bool> sieve(n);
    std::fill(sieve.begin(), sieve.end(), true);
 
    for (size_t i(2); i * i < n; ++i) {
        if (!sieve[i]) continue;
        for (size_t j(i * i); j < n; j += i) {
            sieve[j] = false;
        }
    }
 
    std::vector<int> primes;
    primes.reserve(n >> 1);
    for (size_t i(2); i < n; ++i) {
        if (sieve[i]) primes.push_back(i);
    }
    primes.shrink_to_fit();
 
    return std::move(primes);
}
 
cpp_int product(const std::vector<cpp_int>& multipliers, int i, int j) {
    if (i > j) return 1;
    if (i == j) return multipliers[i];
    int k = (i + j) >> 1;
    return product(multipliers, i, k) * product(multipliers, k + 1, j);
}
 
cpp_int primeSwing(const int& n, const std::vector<int>& primes) {
    std::vector<cpp_int> multipliers;
    multipliers.reserve(n >> 1);
 
    for (int i{}; i < primes.size() && primes[i] <= n; ++i) {
        cpp_int p = 1;
        int q = n;
        while (q != 0) {
            q = q / primes[i];
            if (q % 2 == 1) {
                p *= primes[i];
            }
        }
        if (p != 1) {
            multipliers.push_back(p);
        }
    }
    return product(multipliers, 0, multipliers.size() - 1);
}
 
cpp_int factorial(int n, const std::vector<int>& primes) {
    if (n < 2) return 1;
    cpp_int f = factorial(n / 2, primes);
    cpp_int ps = primeSwing(n, primes);
    return f * f * ps;
}
 
cpp_int factorial(int n) {
    return factorial(n, primes(n));
}
 
 
int main()
{
    std::ofstream os("output.txt");
 
    clock_t tStart = clock();
    os << factorial(3000);
    printf("Time taken: %.4lfs\n", (double)(clock() - tStart) / CLOCKS_PER_SEC);
 
    return 0;
}
 
 
/*
    Time taken: 0.0090s
*/