Pollard rho（Pollard ρ）因数分解算法

#include <bits/stdc++.h>
using namespace std;

using u64 = uint64_t;
using u128 = __uint128_t;

u64 mul_mod(u64 a, u64 b, u64 mod) {
    return (u128)a * b % mod;
}

u64 pow_mod(u64 a, u64 e, u64 mod) {
    u64 r = 1;
    while (e) {
        if (e & 1) r = mul_mod(r, a, mod);
        a = mul_mod(a, a, mod);
        e >>= 1;
    }
    return r;
}

bool is_prime(u64 n) {
    if (n < 2) return false;

    for (u64 p : {2ULL, 3ULL, 5ULL, 7ULL, 11ULL,
                  13ULL, 17ULL, 19ULL, 23ULL, 29ULL,
                  31ULL, 37ULL}) {
        if (n % p == 0) return n == p;
    }

    u64 d = n - 1;
    int s = 0;

    while ((d & 1) == 0) {
        d >>= 1;
        ++s;
    }

    // 对 uint64_t 范围足够的确定性底数
    for (u64 a : {
        2ULL, 325ULL, 9375ULL, 28178ULL,
        450775ULL, 9780504ULL, 1795265022ULL
    }) {
        if (a % n == 0) continue;

        u64 x = pow_mod(a % n, d, n);

        if (x == 1 || x == n - 1) continue;

        bool witness = true;

        for (int r = 1; r < s; ++r) {
            x = mul_mod(x, x, n);

            if (x == n - 1) {
                witness = false;
                break;
            }
        }

        if (witness) return false;
    }

    return true;
}

mt19937_64 rng(chrono::steady_clock::now().time_since_epoch().count());

u64 pollard_rho(u64 n) {
    if (n % 2 == 0) return 2;
    if (n % 3 == 0) return 3;

    while (true) {
        u64 c = uniform_int_distribution<u64>(1, n - 1)(rng);
        u64 x = uniform_int_distribution<u64>(0, n - 1)(rng);
        u64 y = x;
        u64 d = 1;

        auto f = [&](u64 v) {
            return (mul_mod(v, v, n) + c) % n;
        };

        while (d == 1) {
            x = f(x);
            y = f(f(y));

            u64 diff = x > y ? x - y : y - x;
            d = gcd(diff, n);
        }

        if (d != n) return d;
    }
}

void factor(u64 n, vector<u64>& result) {
    if (n == 1) return;

    if (is_prime(n)) {
        result.push_back(n);
        return;
    }

    u64 d = pollard_rho(n);

    factor(d, result);
    factor(n / d, result);
}

int main() {
    ios::sync_with_stdio(false);
    cin.tie(nullptr);

    u64 n;
    cin >> n;

    vector<u64> factors;
    factor(n, factors);

    sort(factors.begin(), factors.end());

    for (int i = 0; i < (int)factors.size();) {
        int j = i;

        while (j < (int)factors.size() &&
               factors[j] == factors[i]) {
            ++j;
        }

        cout << factors[i] << "^" << (j - i);

        if (j != (int)factors.size()) {
            cout << " * ";
        }

        i = j;
    }

    cout << '\n';
    return 0;
}