Pollard Rho 整数分解

Pollard Rho 寻找合数的非平凡因子，再配合素性测试递归分解。在模最小素因子 $p$ 的序列表现近似随机的启发式假设下，通常用 $O(\sqrt p)$ 次迭代找到碰撞；由 $p\le\sqrt n$，得到常说的 $O(n^{1/4})$ 量级。这里统计的是模运算次数，不是关于输入位数的多项式时间保证，也不是任意迭代函数的严格期望界。

Miller–Rabin 素性测试

对奇数 $n>2$，写成 $n-1=d2^s$，其中 $d$ 为奇数。若 $n$ 为素数，则对任意非零底数 $a\bmod n$，必有 $a^d\equiv1\pmod n$，或序列

$$a^d,a^{2d},\ldots,a^{2^{s-1}d}\pmod n$$

中出现 $-1$。否则该底数证明 $n$ 是合数。这来自费马小定理与素数模下平方根性质：$x^2\equiv1$ 只有 $x\equiv\pm1$ 两种解。

单个底数通过测试不能证明任意大整数为素数。但对 $0\le n<2^{64}$，使用固定底数

$$\{2,325,9375,28178,450775,9780504,1795265022\}$$

足以得到确定结果。不能把仅在较小范围有效的一组小素数底数用于全部 64 位整数。固定底数的范围可参见 Miller–Rabin 实现说明。

通过碰撞寻找因子

迭代 $f(x)=(x^2+c)\bmod n$。若两个序列值模 $p$ 相等但模 $n$ 不同，$\gcd(|x-y|,n)$ 就会暴露非平凡因子。

用快慢指针分别每次推进一步和两步，并计算差值的 GCD：结果为 1 时继续，为 $1<g<n$ 时成功，为 $n$ 时本轮失败，重新选择起点和常数。不能把 $n$ 本身送回递归分解，否则递归规模没有缩小。

C++ 实现

以下代码使用 C++17 及 GCC/Clang 的 __uint128_t 扩展，支持 $1\le n<2^{64}$。输入与因子都使用无符号 64 位整数；128 位中间值避免乘法溢出。get_prime_factors(1) 返回空列表，0 不在分解接口的定义域。

#include <algorithm>
#include <cstdint>
#include <numeric>
#include <random>
#include <stdexcept>
#include <vector>

using u64 = std::uint64_t;
using u128 = __uint128_t;

struct PollardRho {
    std::mt19937_64 rng;

    explicit PollardRho(u64 seed = std::random_device{}()) : rng(seed) {}

    static u64 mul_mod(u64 a, u64 b, u64 m) {
        return static_cast<u64>(static_cast<u128>(a) * b % m);
    }

    static u64 power(u64 a, u64 e, u64 m) {
        u64 result = 1;
        while (e) {
            if (e & 1) result = mul_mod(result, a, m);
            a = mul_mod(a, a, m);
            e >>= 1;
        }
        return result;
    }

    static bool miller_rabin(u64 n) {
        if (n < 2) return false;
        for (u64 p : {2ULL, 3ULL, 5ULL, 7ULL, 11ULL, 13ULL}) {
            if (n % p == 0) return n == p;
        }
        u64 d = n - 1;
        int s = 0;
        while ((d & 1) == 0) {
            d >>= 1;
            ++s;
        }
        for (u64 a : {2ULL, 325ULL, 9375ULL, 28178ULL, 450775ULL,
                      9780504ULL, 1795265022ULL}) {
            a %= n;
            if (a == 0) continue;
            u64 x = power(a, d, n);
            if (x == 1 || x == n - 1) continue;
            bool passed = false;
            for (int r = 1; r < s; ++r) {
                x = mul_mod(x, x, n);
                if (x == n - 1) {
                    passed = true;
                    break;
                }
            }
            if (!passed) return false;
        }
        return true;
    }

    u64 get_factor(u64 n) {
        if (n < 2) throw std::invalid_argument("n must be at least 2");
        if (n % 2 == 0) return 2;
        if (n % 3 == 0) return 3;
        if (miller_rabin(n)) return n;
        std::uniform_int_distribution<u64> pick(1, n - 1);
        for (;;) {
            const u64 c = pick(rng);
            u64 x = pick(rng), y = x;
            auto next = [&](u64 v) {
                return static_cast<u64>((static_cast<u128>(mul_mod(v, v, n)) + c) % n);
            };
            // 限制单轮长度；失败或迟迟未找到因子时重新随机选择。
            for (int step = 0; step < 1000000; ++step) {
                x = next(x);
                y = next(next(y));
                u64 diff = x > y ? x - y : y - x;
                u64 g = std::gcd(diff, n);
                if (g == n) break;
                if (g > 1) return g;
            }
        }
    }

    void factorize(u64 n, std::vector<u64>& factors) {
        if (n == 1) return;
        if (miller_rabin(n)) {
            factors.push_back(n);
            return;
        }
        u64 d = get_factor(n);
        factorize(d, factors);
        factorize(n / d, factors);
    }

    std::vector<u64> get_prime_factors(u64 n) {
        if (n == 0) throw std::invalid_argument("zero has no prime factorization");
        std::vector<u64> factors;
        factorize(n, factors);
        std::sort(factors.begin(), factors.end());
        return factors;
    }
};

每次返回的因子都由 GCD 验证，但找到因子所需的重试次数没有固定上限。批量累乘差值后再求 GCD 可以减少开销，不过乘积的 GCD 为 $n$ 时需要逐项回退或重新开始。上述实现采用逐步 GCD，保留较简单的失败处理。