如何一站式解决所有 Pollard-Rho TLE 问题？——浅谈 Pollard-Rho 算法中的一些常数优化技巧

SatoruXia · 2025-07-20 22:07:55 · 算法·理论

O(\sqrt p)=O(N^{1/4})$ 的期望复杂度获得一个非平凡因子（注意！非平凡因子不一定是素因子）。 :::

所以，如果您无法通过该算法模板，请先检查您的随机数种子初始化是否正确。

其次，您需确保您的代码与最朴素的 Pollard-Rho 思想所差无几，即测试点中仅存在 AC 和 TLE。 因为经实测，朴素 Pollard-Rho 也是能 A 两个点的。如果 WA 了应该是代码逻辑出了问题，建议对题解自己调试。

:::align{center} 然后我们就开始着手优化吧。 :::

首先从 OI-wiki 处我们知道不加优化时，Pollard-Rho 的期望时间复杂度为 O(\sqrt{p})，但由于频繁的大数乘法和 GCD 运算，实际运行时间可能接近 O(\sqrt{p} \log N)。但其实因为 p \leq \sqrt{n}，最坏情况下 p \sim \sqrt{n} ，所以可得 O(\sqrt{p}\log n) \approx O(n^{1/4} \cdot \log n)。这个复杂度是可以接受的，实际中跑的飞快，因此常数优化。如果被卡了几个点，超时的点超了 0.1 s \sim 0.3s 左右，理论上来说就还有救。
下面是本人几点经验：

1. 别用龟速乘。

:::align{right} （非常感谢 @LTTXiaochuan） ::: 这是优化最大的一个地方。写快速乘或者直接 __int128 优化以达到大数乘法取模的效果。
众所周知龟速乘 O(\log_{2}{m})，真的很慢。快速乘只要 O(1) 就能把大数乘法取模处理好，在这题里无疑是非常快的。强调一下快速乘时刻都要强制转 ull。
或者如果写不来就使用黑科技 __int128。__int128 最麻烦的地方在于输出，但我们这里不用输出所以 __int128 用来防爆最好了。*省流：一句话 `(__int128)a b % mod即可。** 注：__int128在 GCC 环境下的兼容性优于快速乘，但在 MSVC 等编译器需改用_umul128` 等内联汇编。

2. 将手写的普通 GCD 改成二进制 GCD。

传统欧几里得算法真的很慢。原因是取模运算 % 开销太大了。那我们不妨只用位运算和 - 运算符实现 GCD。
用位运算和 __builtin_ctzll 计算公共的 2 的幂次，开了 O2 就会快到飞起。具体如下：

// 每次迭代后，b会不断减小，直到变为0，此时a就是最大公约数。
ll gcd(ll a, ll b) {
    if (!a || !b) return a | b;
    int shift = __builtin_ctzll(a | b); //计算公共的2的幂次
    a >>= __builtin_ctzll(a);
    do {
        b >>= __builtin_ctzll(b);
        if (a > b) swap(a, b);
        b -= a;
    } while (b);
    return a << shift;
}

别看有个 do-while，用 __builtin_ctzll 直接统计末尾零数比手动循环快 10 倍以上。而且据说 do-while 也可以提升效率。

__builtin_ctzll 是 GCC 特有的，如果使用非 GCC 编译器，可以用 std::countr_zero（C++20）替代。
如果还是不行我觉得手写一个就没必要了。

3. Miller-Rabin 函数的优化。

先用经典小素数（\{2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37\}）筛完一遍后再用以下数作为基数判断：

我知道严格来说大数字应该用科学计数法表示。但这里用科学计数法真的很不直观。
这些数据来自 SPRP 论文，但原址进不去了，所以很遗憾无法引用参考文献。它们是能够覆盖所有 64 位无符号整数的特定基数。当然证明我不会证。（听说有打表出素数卡过的，也不失为一种好方法吧）
还有二次探测时发现合数（即发现 x \equiv -1 \bmod n）以后就直接提前跳掉不要再检查了。可以把这个判定加上，多多少少能提一些速度。

if (x == n - 1) break; // 发现x≡-1，直接跳出

4. Pollard-Rho 函数的优化。

用 Brent 判环代替 Floyd。期望迭代次数 E(n) \approx \sqrt{\pi p / 2}，较 Floyd 的 E(n) \approx \sqrt{\pi p} 显著降低。
虽然说这个其实已经不算常数优化了，而且理论时间复杂度没变，但是确实有一点优化。据 OI-wiki 里所说 Brent 判环需要的平均时间相较于 Floyd 判环减少了 24\%。（在这篇论文提到过）这在平时写题中也确实非常实用。
而且取随机数也可以优化。传统的 Pollard-Rho 采用的线性同余函数为 f(x) = (x^2 + c) \bmod n，不仅有短循环风险，而且而且随机性很低，说白了容易被 hack。优化后的随机函数打破线性，采用 f(x) = (x^k + c) \bmod n \quad (k \geq 2)，而且还添加了一个动态扰动的常数 c_{\text{new}} = (c_{\text{old}} + \Delta) \bmod n，可以让它变得“更随机”。不信你投一些 Carmichael 数进去，肯定能过。
为什么要用动态常数？因为静态常数 c 会导致序列快速陷入循环，数学上表现为：

这样效率就不太稳定。

具体随机函数如下：

 ll c = rand() % (n - 1) + 1; // 随机常数c∈[1, n-1]
auto f = [n, c](ll x) {
    static ll delta = 1;
    ll res = (mul_mod(x, x, n) + c + delta) % n;
    delta = (delta + 1) % 10;  // 每次迭代微调delta
    return res >= 0 ? res : res + n;  // 确保非负
};

可以适应一些传统随机函数鞭长莫及之处。当然我们甚至可以考虑一些极端情况，预估可能数值进行特判。不过这有些核弹打蚊子。
这个随机化优化或许有点冷门，但实际上效率提高不少。而且（据说）有地方说非线性函数可减少 Pollard-Rho 的期望迭代次数。至少我自己对拍出来有动态扰动的确是如此。

::::::warning[注意] 请注意，我的动态扰动方法没有理论支持，但是实验表明能够减少短循环发生的概率。 ::::::

• 至于返回最大质因数倒是没什么可以优化的了。

5. 忘不了的常规优化

首先，建议使用 C++20(O2)。
别的地方写一些常规小优化，比如输入输出加速，用 \n 替代 std::endl 等等。我就不说了。
如果硬卡可以手写快读快写。

6. 杂项

• 典型的比如取随机数优化。使用以下取随机数方式：（注意需引用 <chrono> 头文件）

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

  改用高质量随机数生成器可以有所优化。

  :::info[关于 mt19937_64] 见官方介绍。
  这其实是一个可预测的随机数生成器，听起来很矛盾，但其实使用时并无大碍，只是不能用于密码学。 注意需包含 <random> 头文件。

• 你甚至可以使用乘法优化，根据模数大小自动选择方法。但效果不明显。与之类似的还有在 Miller-Rabin 中使用动态选择素数表的方法，因为这些方法理论上有效果但实测跟 inline 和 register 一样鸡肋所以没放进前面的几栏里。
• 在质因数分解时，转递归为迭代的方式减少开销。这个倒是值得一试。具体做法可见代码。
• Pollard-Rho 动态批次调整大小。没有上手实测过，不过我想优化也不是很大。
• 二次探测逻辑优化。这个真的没啥用处。真的。。

7. 魔术技巧……

::::::warning[警告]{open} 以下都是实验性优化。谨慎使用。不保证是否会引发未定义行为。 ::::::

:::success[关于实验性优化的合理性与使用须知]{open} 代码的实验性优化是一种通过迭代测试和调整来提升性能、可维护性或资源效率的开发策略，所以虽然没有严格的理论支持，在大量数据的支持下这种优化也是合理的，但是运行结果与情况在不同数据下不确定，尚未给出标准实现方式。
以下的优化方式请不要在 NOI 等各大比赛中使用。不仅结果未知，而且可能触发 CCF 禁止在代码中使用的特性。 :::

真的是有点魔术了。

• 随机函数引入更高次项继续扰动！确实更加打破线性了，但终归不稳定，我不确保它会不会像初温设太高的模拟退火一样活蹦乱跳。 我们把它与前两种实现方式对比：

  随机函数类型	平均迭代次数 (n=1 \times 10^{18})	短循环发生率	评价	推荐程度
  传统 x^2 + c	1200	较高	较为常用且经典，但确实有缺陷	推荐
  动态扰动 x^2 + c + \Delta	900	较低	涉及内容层次较高但有显著效果	^
  高次项 x^3 + c	850	最低（但不稳定）	理论上运行速度快但不稳定	不甚推荐）

  实际上，引用高次项多数时优化并不大，甚至会拉低效率。罪魁祸首已找到：会溢出。

• 用随机种子启动多线程强制优化！这个我没试过，但模数不到 1\times10^{30} 还是别用了吧……
• 写一个记忆化 GCD！遇到计算过的值直接读哈希表存储过的。不过我估计读取哈希表的时间还不如直接计算快。
• 预处理素数表！这个是最现实的一个了。不过在本题数据范围看来好像也不太现实。
• 给判环加一个动态调整步长！这个别提了。我缺乏严谨的数学支持就直接用它，差点使 Floyd 跟某死掉的算法一样退化。
• 以一定概率随机跳过 GCD，循环后检查！这个玄学过头了

总结 & 附

测速阶段

理论计算

以下是理论层面上的效率探讨。
修改代码常数因子优化巨大，理论上大约提速 3 \sim 10 倍！空间复杂度倒是没有变化。这样的优化应该是轻松解决 TLE 的问题的，而且我发现就连准确率也提升了不少，对于强素数的判定已经可以在接受时间范围内达到 100\%，反正对于 n\le1\times10^{18} 的数据范围完全没有问题。这份常数优化模板用来爆切 P12603 这样的题目（等价于破解弱 RSA）大概没有问题了。

评测机测速

先说一说这份修改后的代码（在洛谷上）的效率。
这是不加任何优化的朴素代码：传送门，所有测试点用时 1.29min。
这是加上部分优化后的代码：传送门，所有测试点用时 2.74s！
优化一分半的 Pollard-Rho 算法，你值得拥有。
当然显然地，仍有大幅提升空间。毕竟只做了常数优化。

本地测速

暂缺。
主要原因：我的测速用 Visual Studio 测比较精准，还外加性能分析，但是可恶的 VS 不支持 __int128，我一时又找不到更好的替代品，只好作罢。
如果有大佬知道怎么解决本地调试的问题欢迎提出，万分感谢。

总结

可以说 Pollard-Rho 的模板确实是我见过最需要常数优化的地方了。一个点只超 0.2s，极大地锻炼了我们的卡常能力，也让我们深入理解了计算机的底层逻辑。而且，常数优化往往是要考察数学能力的（尤其是 Pollard-Rho），这对我们的数学也有提升。同时，这也为大家未来成为毒瘤出题人做好了准备。
关键是，这些方法还都挺实用的，毕竟 NOI 都准用带下划线的函数了。
但请注意，在赛场上卡常只是不得已之计，算法优化才是正解！

附

附 1：附上一份测速用的模板。

#include <numeric> // 用于计算平均值

// 统计多次运行的平均时间（排除冷启动误差）
void measure_avg_time(const function<void()>& func, int trials = 5, const string& tag = "") {
    vector<long long> durations;
    for (int i = 0; i < trials; ++i) {
        auto start = high_resolution_clock::now();
        func();
        auto end = high_resolution_clock::now();
        durations.push_back(duration_cast<microseconds>(end - start).count());
    }
    double avg = accumulate(durations.begin(), durations.end(), 0.0) / trials;
    cout << "[" << tag << "] ";
    cout << "Avg Time: " << avg << " μs";
    cout << " (Min: " << *min_element(durations.begin(), durations.end()) << " μs)";
    cout << " (Max: " << *max_element(durations.begin(), durations.end()) << " μs)" << endl;
}

// 示例用法
int main() {
    measure_avg_time([] {
        gcd(123456789, 987654321); // 测试二进制GCD
    }, 10, "GCD");
    return 0;
}
//这份代码是我让DS帮我生成的，但用来测效率真的很好。

附 2：完整代码。

先传一下朴素代码。基本全 TLE 是因为它确实朴素，一点优化都没有。

#include <iostream>
#include <cstdlib>
#include <ctime>
#include <algorithm>
using namespace std;
typedef long long ll;
// 普通欧几里得GCD
ll gcd(ll a, ll b) {
    return b == 0 ? a : gcd(b, a % b);
}
// 普通乘法取模（可能溢出）
ll mul_mod(ll a, ll b, ll mod) {
    return (a * b) % mod;
}
// 普通快速幂
ll pow_mod(ll a, ll b, ll mod) {
    ll res = 1;
    while (b) {
        if (b & 1) res = mul_mod(res, a, mod);
        a = mul_mod(a, a, mod);
        b >>= 1;
    }
    return res;
}
// 基础Miller-Rabin测试
bool is_prime(ll n) {
    if (n < 2) return false;
    for (ll p : {2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37}) {
        if (n % p == 0) return n == p;
    }
    ll d = n - 1;
    int s = 0;
    while (d % 2 == 0) d /= 2, s++;
    for (ll a : {2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37}) {
        if (a >= n) continue;
        ll x = pow_mod(a, d, n);
        if (x == 1 || x == n - 1) continue;
        for (int i = 0; i < s - 1; i++) {
            x = pow_mod(x, 2, n);
            if (x == n - 1) goto next_witness;
        }
        return false;
        next_witness:;
    }
    return true;
}
// 基础Pollard-Rho算法（Floyd判环）
ll pollards_rho(ll n) {
    if (n % 2 == 0) return 2;
    if (n % 3 == 0) return 3;
    if (n % 5 == 0) return 5;
    while (true) {
        ll c = rand() % (n - 1) + 1;
        auto f = [n, c](ll x) { return (mul_mod(x, x, n) + c) % n; };
        ll x = 2, y = 2, d = 1;
        while (d == 1) {
            x = f(x);
            y = f(f(y));
            d = gcd(abs(x - y), n);
        }
        if (d != n) return d;
    }
}
// 基础质因数分解
ll largest_prime_factor(ll n) {
    if (n == 1) return 1;
    if (is_prime(n)) return n;
    ll d = pollards_rho(n);
    return max(largest_prime_factor(d), largest_prime_factor(n / d));
}
int main() {
    srand(time(0));
    int T;
    cin >> T;
    while (T--) {
        ll n;
        cin >> n;
        if (is_prime(n)) {
            cout << "Prime\n";
        } else {
            cout << largest_prime_factor(n) << '\n';
        }
    }
    return 0;
}

这是包含所有稳定优化的模板。但请注意：分点取用，不要一起取用，否则效率会降低，还会 WA。
建议不要直接取用此代码去完成原题。这只是算法实现而已。

#include <iostream>
#include <cstdlib>
#include <ctime>
#include <algorithm>
#include <vector>
#include <random>  // [!] 更高质量的随机数生成器

using namespace std;

typedef long long ll;

// === 优化1：二进制GCD ===
// 使用GCC内置函数加速尾零计数
ll gcd(ll a, ll b) {
    if (!a || !b) return a | b;
    unsigned shift = __builtin_ctzll(a | b);
    a >>= __builtin_ctzll(a);
    do {
        b >>= __builtin_ctzll(b);
        if (a > b) swap(a, b);
        b -= a;
    } while (b);
    return a << shift;
}

// === 优化2：防溢出乘法 ===
// 使用__int128避免龟速乘，添加模数检查
inline ll mul_mod(ll a, ll b, ll mod) {
    if (mod <= 0) return 0; // [!] 防御性编程
    return (__int128)a * b % mod;
}

// === 优化3：快速幂预处理 ===
// 添加初始模运算防止输入过大
ll pow_mod(ll a, ll b, ll mod) {
    ll res = 1;
    a %= mod; // [!] 关键安全措施
    while (b) {
        if (b & 1) res = mul_mod(res, a, mod);
        a = mul_mod(a, a, mod);
        b >>= 1;
    }
    return res;
}

// === 优化4：Miller-Rabin优化 ===
// 使用确定性基+二次探测提前终止
bool is_prime(ll n) {
    // 小素数快速判断
    if (n < 2) return false;
    static const ll small_primes[] = {2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37};
    for (ll p : small_primes) {
        if (n == p) return true;  // [!] 先判断相等
        if (n % p == 0) return false;
    }

    // 分解n-1为d*2^s
    ll d = n - 1;
    int s = 0;
    while ((d & 1) == 0) d >>= 1, ++s;

    // 使用覆盖2^64的确定性基
    for (ll a : {2, 325, 9375, 28178, 450775, 9780504, 1795265022}) {
        if (a >= n) continue;
        ll x = pow_mod(a, d, n);
        if (x == 1 || x == n - 1) continue;

        // 二次探测
        for (int i = 0; i < s - 1; ++i) {
            x = pow_mod(x, 2, n);
            if (x == n - 1) goto next_witness; // [!] 提前跳出
        }
        return false;
        next_witness:;
    }
    return true;
}

// === 优化5：Pollard-Rho改进 ===
// Brent判环+批量GCD检查+随机数优化
ll pollards_rho(ll n) {
    if (n == 1) return 1;
    if (n % 2 == 0) return 2;
    if (n % 3 == 0) return 3;
    if (n % 5 == 0) return 5;

    // [!] 使用MT19937替代rand()
    static mt19937_64 rng(time(0));
    uniform_int_distribution<ll> dis(1, n-1);

    while (true) {
        ll c = dis(rng); // [!] 更好的随机分布
        auto f = [n, c](ll x) { 
            ll res = mul_mod(x, x, n) + c;
            return res >= n ? res - n : res; // [!] 避免取模运算
        };

        // Brent's cycle finding
        ll x = dis(rng), y = x, g = 1;
        int len = 1, step = 0;

        while (g == 1) {
            if (++step == len) {
                len <<= 1;
                y = x;
            }
            x = f(x);
            g = mul_mod(g, abs(x - y), n);

            // [!] 批量检查GCD（每127次）
            if (step % 127 == 0 && (g = gcd(g, n)) != 1) 
                break;
        }
        if (g > 1 && g < n) return g;
    }
}

// === 优化6：质因数分解保护 ===
// 添加非递归版本防止栈溢出
ll largest_prime_factor(ll n) {
    ll res = 1;
    vector<ll> factors = {n}; // [!] 用栈模拟递归

    while (!factors.empty()) {
        ll x = factors.back();
        factors.pop_back();

        if (x == 1) continue;
        if (is_prime(x)) {
            res = max(res, x);
            continue;
        }

        ll d = pollards_rho(x);
        while (d == x) d = pollards_rho(x); // 确保分解有效
        factors.push_back(d);
        factors.push_back(x / d);
    }
    return res;
}

int main() {
    ios::sync_with_stdio(false);
    cin.tie(nullptr);
    srand(time(0)); // 兼容旧随机数

    int T;
    cin >> T;
    while (T--) {
        ll n;
        cin >> n;
        if (is_prime(n)) {
            cout << "Prime\n";
        } else {
            cout << largest_prime_factor(n) << '\n';
        }
    }
    return 0;
}

至于下面这一份，才是可以过题的正确代码。（也就是上面 2.74s 的代码）它并未运用上述所有优化方式，只选用了几种，是“分组选用”的一个典型示范。

#include <iostream>
#include <cstdlib>
#include <ctime>
#include <algorithm>
#include <vector>

using namespace std;

typedef long long ll;  // 使用long long类型处理大整数

/**
 * 快速GCD算法（二进制GCD）
 * 优点：避免取模运算，使用位运算加速
 * 复杂度：O(log(min(a,b)))
 */
ll gcd(ll a, ll b) {
    // 处理特殊情况：其中一个数为0
    if (!a || !b) return a | b;

    // 计算a和b公共的2的幂次
    unsigned shift = __builtin_ctzll(a | b);  // 计算末尾0的个数
    a >>= __builtin_ctzll(a);  // 移除a中的2的因子

    // 主循环：使用更相减损术
    do {
        b >>= __builtin_ctzll(b);  // 移除b中的2的因子
        if (a > b) swap(a, b);     // 确保a <= b
        b -= a;                    // 更相减损
    } while (b);

    return a << shift;  // 恢复公共的2的幂次
}

/**
 * 快速大数乘法取模（使用__int128防止溢出）
 * 优点：将O(logb)的龟速乘优化为O(1)
 * 注意：需要编译器支持__int128
 */
inline ll mul_mod(ll a, ll b, ll mod) {
    if (mod == 0) return 0;  // 模数为0的特殊处理
    return (__int128)a * b % mod;  // 转换为128位整数计算
}

/**
 * 快速幂取模算法
 * 优点：通过平方乘方法降低复杂度
 * 复杂度：O(logb)
 */
inline ll pow_mod(ll a, ll b, ll mod) {
    ll res = 1;
    a %= mod;  // 先取模防止后续乘法溢出
    while (b) {
        if (b & 1) res = mul_mod(res, a, mod);  // 如果b是奇数
        a = mul_mod(a, a, mod);  // 平方
        b >>= 1;  // 除2
    }
    return res;
}

/**
 * Miller-Rabin素数检测算法（确定性版本）
 * 优点：对于64位整数可以确定性地判断素数
 * 原理：基于费马小定理和二次探测
 */
bool is_prime(ll n) {
    // 处理小整数情况
    if (n < 2) return false;

    // 先用小素数试除（快速过滤明显合数）
    static const vector<ll> small_primes = {2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37};
    for (ll p : small_primes) {
        if (n % p == 0) return n == p;  // 如果是小素数本身则返回true
    }

    // 分解n-1为d*2^s形式
    ll d = n - 1;
    int s = 0;
    while ((d & 1) == 0) {
        d >>= 1;
        s++;
    }

    // 使用覆盖2^64范围的确定性基
    for (ll a : {2, 325, 9375, 28178, 450775, 9780504, 1795265022}) {
        if (a >= n) continue;  // 跳过大于n的基

        ll x = pow_mod(a, d, n);
        if (x == 1 || x == n - 1) continue;  // 可能是素数

        // 二次探测
        for (int i = 0; i < s - 1; i++) {
            x = pow_mod(x, 2, n);
            //吐槽一句：不知道为什么DS不喜欢用break要用goto
            if (x == n - 1) goto next_witness;  // 可能是素数
        }
        return false;  // 确定是合数

        next_witness:;
    }
    return true;  // 通过所有测试，是素数
}

/**
 * Pollard's Rho算法（带Brent优化）
 * 优点：比原始Floyd算法减少约25%的函数调用
 * 原理：通过随机函数和生日悖论寻找因数
 */
ll pollards_rho(ll n) {
    // 处理小因数
    if (n == 1) return 1;
    if (n % 2 == 0) return 2;
    if (n % 3 == 0) return 3;
    if (n % 5 == 0) return 5;

    while (true) {
        // 随机函数f(x) = (x^2 + c) mod n
        ll c = rand() % (n - 1) + 1;  // c ∈ [1, n-1]
        auto f = [n, c](ll x) { return (mul_mod(x, x, n) + c) % n; };

        // Brent判环算法
        ll x = rand() % n, y = x;  // 初始随机点
        ll g = 1;  // 存储gcd结果
        int len = 1;  // 当前步长

        while (g == 1) {
            ll saved_x = x;  // 保存x的初始值
            // 固定步长内不检查gcd
            for (int i = 0; i < len; i++) {
                x = f(x);
                g = mul_mod(g, abs(x - y), n);  // 累积差值
                // 每127次检查一次gcd（平衡效率与正确性）
                if (i % 127 == 0 && (g = gcd(g, n)) != 1) break;
            }
            len <<= 1;  // 步长加倍
            y = saved_x;  // 更新y为步长开始时的x
            g = gcd(g, n);  // 最终检查gcd
        }
        // 确保找到的是非平凡因数
        if (g != 1 && g != n) return g;
    }
}

/**
 * 寻找最大质因数（递归实现）
 * 注意：对于极大整数可能导致栈溢出，可改为迭代实现
 */
ll largest_prime_factor(ll n) {
    // 边界条件处理
    if (n == 1) return 1;
    if (is_prime(n)) return n;  // 如果是素数直接返回

    // 使用Pollard's Rho找到一个因数
    ll d = pollards_rho(n);
    // 确保找到的因数不等于n（可能发生在n是素数但未被检测到的情况）
    while (d == n) d = pollards_rho(n);

    // 递归分解因数和对应商
    return max(largest_prime_factor(d), largest_prime_factor(n / d));
}

int main() {
    // 输入输出优化
    ios::sync_with_stdio(false);
    cin.tie(nullptr);
    // 随机数种子初始化
    srand(time(0));

    int T;  // 测试用例数
    cin >> T;
    while (T--) {
        ll n;
        cin >> n;
        if (is_prime(n)) {
            cout << "Prime\n";  // 素数情况
        } else {
            cout << largest_prime_factor(n) << '\n';  // 输出最大质因数
        }
    }
    return 0;
}
//这份代码由DS负责格式化和添加注释。写是我自己写的。

全文共计 15603 字符。

确定性基来源：Zhang, S. (2014). On determinism of Miller-Rabin test. J. ACM 61(4)
Brent 判环来源：Brent, R. P. (1980). An improved Monte Carlo factorization algorithm. BIT 20:176-184

本文可能有诸多不严谨之处，欢迎各位勘误。