Fast factorial computation

#define _CRT_SECURE_NO_WARNINGS
#pragma GCC Optimize("Ofast")
#include <string.h>
#include <iostream>
#include <vector>
#include <set>
#include <chrono>
#include <string>
#include <time.h>
#include <unordered_set>
#include <iomanip>
#include <cmath>
#include <map>
#include <queue>
#include <numeric>
#include <unordered_map>
#include <algorithm>
#include <fstream>
#include <sstream>
#include <cstddef>
#include <cstdio>
#include <iostream>
#include <memory>
#include <stdexcept>
#include <string>
#include <cmath>
#include <complex>
#include <array>
using namespace std;
typedef long long ll;
const int maxn = 1e2 + 14, lg = 15;

const int BASE = 3e3;
int64_t count_ops = 0;
const double PI = acos(-1);
typedef complex<double> base;

double SIN[1000000];
#define MAXN 10000000

void fft(vector<base>& a, bool invert) {
    int n = (int)a.size();

    for (int i = 1, j = 0; i < n; ++i) {
        int bit = n >> 1;
        for (; j >= bit; bit >>= 1)
            j -= bit;
        j += bit;
        if (i < j)
            swap(a[i], a[j]);
    }

    for (int len = 2; len <= n; len <<= 1) {
        double ang = 2 * PI / len * (invert ? -1 : 1);
        base wlen(cos(ang), sin(ang));
        for (int i = 0; i < n; i += len) {
            base w(1);
            for (int j = 0; j < len / 2; ++j) {
                base u = a[i + j], v = a[i + j + len / 2] * w;
                a[i + j] = u + v;
                a[i + j + len / 2] = u - v;
                w *= wlen;
            }
        }
    }
    if (invert)
        for (int i = 0; i < n; ++i)
            a[i] /= n;
}

vector<int> multiply(const vector<int>& a, const vector<int>& b) {
    vector<base> fa(a.begin(), a.end()), fb(b.begin(), b.end());
    size_t n = 1;
    while (n < max(a.size(), b.size()))  n <<= 1,
        ++count_ops;
    n <<= 1;
    fa.resize(n), fb.resize(n);

    fft(fa, false), fft(fb, false);
    for (size_t i = 0; i < n; ++i)
        fa[i] *= fb[i],
        ++count_ops;
    fft(fa, true);
    vector<long long> res(n);
    for (size_t i = 0; i < n; ++i) {
        res[i] = (long long)(fa[i].real() + 0.5),
            ++count_ops;
        if (res[i] < 0 || res[i] > 1e15) {
            cout << -1;
            res[i] = res[i];
        }
    }
    int carry = 0;
    for (size_t i = 0; i < n; ++i) {
        res[i] += carry;
        carry = res[i] / BASE;
        ++count_ops;
        res[i] %= BASE;
    }
    vector<int> r2(res.begin(), res.end());
    while (r2.size() && r2.back() == 0) r2.pop_back();
    if (!r2.size()) r2.push_back(0);
    return r2;
}

#define typ vector<int>

const int N = 1000000;
int lp[N + 1];
vector<int> pr;

typ prod(vector<typ>& a, int l, int r) {
    if (l > r) return { 1 };
    if (l == r) return { a[l] };
    if (r - l == 1) return multiply(a[l], a[r]);
    return multiply(prod(a, l, (l + r) >> 1), prod(a, ((l + r) >> 1) + 1, r));
}

typ primeSwing(int n) {
    vector<typ> mult;
    for (int i = 0; pr[i] <= n; ++i) {
        int q = n;
        typ prime = { pr[i] };
        typ val = { 1 };
        while (q) {
            q /= pr[i];
            if (q & 1) {
                val = multiply(val, prime);
            }
        }
        if (prime.size() > 1 || prime[0] != 1) {
            mult.push_back(val);
        }
    }
    return prod(mult, 0, mult.size() - 1);
}

typ factorial(int n) {
    if (n < 2) {
        return { 1 };
    }
    typ f = factorial(n / 2);
    return multiply(multiply(f, f), primeSwing(n));
}

void sieve() {
    for (int i = 2; i <= N; ++i) {
        if (lp[i] == 0) {
            lp[i] = i;
            pr.push_back(i);
        }
        for (int j = 0; j < (int)pr.size() && pr[j] <= lp[i] && i * pr[j] <= N; ++j)
            lp[i * pr[j]] = pr[j];
    }
}

signed main() {
    cin.tie(0);
    cout.tie(0);
    ios::sync_with_stdio(false);
    sieve();
    std::chrono::steady_clock::time_point begin = std::chrono::steady_clock::now();
    auto it = factorial(100000);

    cout << it.size() * 3;
    std::chrono::steady_clock::time_point end = std::chrono::steady_clock::now();
    std::cout << "Time difference = " << std::chrono::duration_cast<std::chrono::milliseconds>(end - begin).count() << "[µs]" << std::endl;
    cout << count_ops;

}