Length of longest subsequence having gcd > 1, C++ (gcc)

Language: Layout:

#include <iostream>
#include <cstdint>
#include <set>
#include <cmath>
#include <chrono>
#include <vector>
#include <random>
#include <unordered_map>
#include <algorithm>
#include <ctime>

std::set<int64_t> primeFactors(int64_t);

int main() {
    std::mt19937 prbg{(unsigned int)time(0)};
    auto start = std::chrono::steady_clock::now();
    std::unordered_map<int64_t, int> freq;
    for (int i = 0; i < 1000; ++i) {
        int64_t n = std::uniform_int_distribution<int64_t>{
            2LL, 1'000'000'000'000'000'000}(prbg);
        for (auto p : primeFactors(n)) freq[p]++;
    }
    int result = std::max_element(begin(freq), end(freq), [](auto& a, auto& b) {
                     return a.second < b.second;
                 })->second;
    std::chrono::duration<double, std::milli> elapsed =
        std::chrono::steady_clock::now() - start;
    std::cout << result << "\n" << elapsed.count() << "ms\n";
}

// compute a^d mod n
int64_t modPow(int64_t a, int64_t d, int64_t n) {
    int64_t ret = 1;
    int64_t apow = a;
    for (; d; d >>= 1) {
        if (d % 2) ret = __int128(ret) * apow % n;
        apow = __int128(apow) * apow % n;
    }
    return ret;
}

bool witness(int64_t a, int64_t n, int64_t d, int64_t r) {
    auto x = modPow(a, d, n);
    if (x == 1 || x == n - 1) return true;
    for (auto i = r - 1; i-- > 0;) {
        x = modPow(x, 2, n);
        if (x == 1) return false;
        if (x == n - 1) return true;
    }
    return false;
}

bool isPrimeDetMillerRabin(int64_t n) {
    int64_t d, r; // n-1 = d * 2^r
    for (d = n - 1, r = 0; d % 2 == 0; d >>= 1, r++) {}
    const int64_t arr[] = {2, 325, 9375, 28178, 450775, 9780504, 1795265022};
    for (auto a : arr)
        if (!witness(a, n, d, r)) return false;
    return true;
}

template <class T>
T gcd(T a, T b) {
    return b ? gcd(b, a % b) : a;
}

int64_t pollardsRho(int64_t n) {
    auto g = [n](int64_t x) { return (int64_t)((__int128(x) * x + 1) % n); };
    int64_t x = 2, y = 2, d = 1;
    while (d == 1) {
        x = g(x);
        y = g(g(y));
        d = gcd(x > y ? x - y : y - x, n);
    }
    // if (d == n) throw "Failture";
    return d;
}

std::set<int64_t> primeFactors(int64_t n) {
    std::set<int64_t> result;
    // 2
    if (n % 2 == 0) result.insert(2);
    while (n % 2 == 0) n /= 2;
    // odds
    const int iMax = 31623;
    for (int sqrtn = (int)sqrt(n), i = 3; i < sqrtn && i < iMax; i += 2) {
        if (n % i == 0) result.insert(i);
        while (n % i == 0) n /= i;
        sqrtn = (int)sqrt(n);
    }
    while (n > 1) {
        // MR
        if (isPrimeDetMillerRabin(n)) {
            result.insert(n);
            break;
        } else {
            // Pollard's rho
            auto d = pollardsRho(n);
            result.insert(d);
            n /= d;
        }
    }
    return result;
}

Show compiler warnings [ + ] Compiler args [ + ] Show input
Compilation time: 1.13 sec, absolute running time: 0.73 sec, cpu time: 0.68 sec, memory peak: 3 Mb, absolute service time: 1,96 sec