Miller-Rabin 与 Pollard-Rho

naoliaok_lovely · 2023-12-12 23:18:11 · 算法·理论

Miller-Rabin 素性测试

就像他的名字，这个算法是由 Miller 和 Rabin 二人所提出，用于判断一个数是不是质数。

在此之前，我们使用的都是复杂度为 O\left(\dfrac{\sqrt n}{\ln n}\right) 的方法，但是它最多也只能跑 10^{15} 的数量级，再大的数就很吃力了。而在平时的生活中，很可能会遇到更大的数的分解。这时，我们以前的分解方法不再使用，需要更快的方法，而 miller-rabin 就可以做到这一点。

事实上，素数判定为 P 问题，目前已知的确定性素数判断方法为 APRCL，复杂度 O((\log n)^{c\log\log\log n})。很显然，他也不是多项式复杂度的，但是可以用于约 1000 位的素数判断。
既然如此，我们的 miller-rabin 算法便选择了“随机”的策略。一般在工业中还有一种非确定性算法叫 ECCP，不过他是绝对正确的，“非确定”指的是因为随机导致的程序运行效率不同，为 O(\log^{5+\epsilon}n)。

费马测试

在介绍 miller-rabin 算法之前，先让我们看看费马测试。我们都知道，对于质数 p，那么有费马小定理 a^{p-1}\equiv1\pmod p 存在，其中 \gcd(a,p)=1。那么，如果 a 满足 \gcd(a,p)=1，存在 a^{p-1}\equiv1\pmod p 成立，是否可以说明 p 为质数呢？再换句话说，费马小定理的逆命题是否成立呢？

这是一个很不错的想法，但很遗憾他是错误的……不过幸运的是，在大多数情况下，这个方法都是正确的。
于是我们可以尝试进行改进——既然找一个底数不行，那多找几个呢？可惜的是有一类数，他们是合数，却也满足上面费马小定理的性质！这些数被称为 Carmichael 数，不过分布十分稀疏，在 10^{18} 以内大约有 1.4\times10^6 个。（不过总不可能打表吧啊喂）

性质：Carmichael 数至少有 3 个不同的质因子。

二次探测定理

如果 p 为奇素数，满足 a^2\equiv1\pmod p，则 a\equiv\pm1\pmod p。

考虑一个问题：什么样的合数可以通过二次探测？
（易证）
结论：当且仅当 n=p^k 时。

改进

Miller-Rabin 就是将二者结合了一下。事实上， Carmichael 数有一个性质：其一定不是 p^k 。所以便产生了如下的方式：

将 p 分解为 2^km，即找出当中的偶数质因子。
对于底数 a，计算 a^m\bmod p，然后不断平方。由二次探测定理得，设当前的余数为 r，且 r^2\equiv1\pmod p，则 r\equiv\pm1\pmod p。如果出现矛盾则不是质数。
如果最终的余数不为 1 则不是质数。

多次随机底数 a 判定即可。

时间复杂度：
在 OI 中，miller-rabin 常常用来计算 long long 范围内的数，所以乘法和取模均可看做 O(1)，所以复杂度为 O(T\log n)。但是在工业中，需要用到高精度，乘法和取模是 O(\log^2n) 的，所以总复杂度为 O(T\log^3n)，这时需要使用 FFT 将其优化至 O(T\log^2n)。

优化

上面的随机看着就不靠谱。事实上，对于 int 范围内的数，我们只需要验证 \{2,7,61\} 即可；对于 long long 范围内的数，我们只需要验证前 12 个质数即可。

bool miller_rabin(LL n)
{
    if(n <= 1 || n > 2 && ~n & 1) return false;
    LL m = n - 1, k = 0;
    while(~m & 1) m >>= 1, k++;
    vector<int> s = {2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37};
    for(int a : s)
    {
        if(a >= n) break;
        i128 x = ksm(a, m, n), y;
        for(int i = 1; i <= k; i++, x = y)
        {
            y = x * x % n;
            if(y == 1 && x != 1 && x != n - 1) return false;
        }
        if(x != 1) return false;
    }
    return true;
}

Rollard-Rho

P4718 【模板】Pollard-Rho
可以在 O(\sqrt[4]n) 的复杂度以内分解因数。

质因数分解本身为 NP 问题、反 NP 问题、BQP 问题。

生日悖论

先引入生日悖论，不过相信很多人都知道，即：

在一个班的同学中，存在两个人生日相同的概率是多少？

相信大家都很熟悉，这里就不证明了，总之一句话就是，看起来概率很小，实则概率接近 100%。也正因如此，被称为了“悖论”。

随机与 \gcd

同样考虑我们以前的分解方法。最优的就是先枚举 \sqrt n 以内的所有质因数，然后再组合，复杂度 O(\frac{\sqrt n}{\ln n})。这似乎已经走到头了，再想优化复杂度就只能随机了。

但是，随机一个数简直太荒谬了，可以说运气不好的话跑几分钟甚至几十分钟绰绰有余。对于我们随机到的数字 a，当 a|n 的时候，我们才去递归处理 a 和 n/a，这概率确实很小。
我们可以采用求 \gcd 的方式增大这个概率：如果 \gcd(a,n)>1，我们是不是就找到了一个因数呢？考虑最劣的情况：n=t_1\times t_2，其中 t_1,t_2 都接近于 \sqrt n。此时，t1,t2 的倍数约为 O(\sqrt n)，那么总复杂度也就降低到了 O(\sqrt n)。可惜还得加上 \gcd 的复杂度，所以复杂度为 O(\sqrt n\log n)。

随机函数

Pollard-Rho 最精妙的一点就在于它的随机函数。首先，我们随机 d,x_1，那么序列 x 由 x_i=(x_{i-1}^2+d)\bmod n 给出。即，我们创造了一个伪随机序列。

举个例子，设 n=50,d=2,x_1=1，则序列为 1,3,11,23,31,11,23,31\dots。显然，无论如何最后都会进入一个循环，即下面的形状：（摘自 OI-wiki）
这也是其名字 \rm rho\to\rho 的由来。

生成函数的性质：若 p|n,x_i\equiv x_j\pmod p，则 x_{i+k}\equiv x_{j+k}\pmod p。
所以我们只需判断下标之差不同的数对。

于是，对于下标的差不同的两项，我们取他们的差，在执行 \gcd 的判断即可。Pollard-Rho 给出了这样的构造方式并说明了这样构造后 \gcd>1 的期望随机次数不大。最后给出的结论为，总随机次数约 O(\sqrt[4]n)。

floyd 判环

先让我们根据上面的理论，总结一下实现步骤：

如果 n 不为合数，直接返回。
定义上述的递推函数（可以用下面的 Lambda 表达式写法，看起来清新一点）。
随机 d,x_1，然后进行循环，每次求出下标差不同的两项的差。如果 \gcd>1，那么递归分解。

可惜这个方法并不完美，存在两个很大的缺陷。其一是如何寻找下标差不同的数对，这个先放着，重点在第二个。仔细思考上面的 \rho 形，无论从哪一个数开始最后总会步入循环——那万一运气不好，每一组差都刚好互质呢？
这个情况显然是会出现的，并且很容易构造。考虑修正，一种比较直观的方法是，一旦发现出现了环，那么就重新随机一个 \rho 形。但是怎么判环呢？值域达到了 10^{18} 级别，如果直接用 unordered_map 那又要带一个大常数，手写哈希表当然也可以，但是空间又不好把握。

这里引入 floyd 判环的思想，可以同时规避以上两个问题。简单点说就是在“龟兔赛跑”。我们设置两个变量 x,y，其中 x 是“乌龟”，y 是兔子。我们规定，x 每次跳一步，y 每次跳两步，将上面“相邻两项的差”更改为“x,y 的差”。一旦 x=y，就说明 y 一定比 x 多跑了至少 1 圈，此时说明出现了环，重新随机即可。更改一下上面的步骤：

随机 d,x，令 y=(x^2+d)\bmod n。重复执行，如果 \gcd(|x-y|,n)>1 则递归分解，否则 x 跳一步，y 跳两步。
如果 x=y，重新随机 d,x。

优化

仔细思考，如果对于每次我们即将判断的数都取一遍 \gcd，那实际复杂度应该为 \sqrt[4]n\log n，这个 \log 在平时的运用可能无所谓，但是多了一个 50 的常数对于一些卡常题（比如模板）肯定是过不去的。事实上，我们没有必要每次都求一个 \gcd，可以设定一个阈值，每次取多个数相乘，用乘积取 \gcd，即减少计算 \gcd 的次数。
在题解中有一个倍增优化的东西，是在对这个阈值进行倍增处理，即可优化掉这个 \log。不过我一般喜欢分块优化，直接将这个阈值设置为常数（一般设为 32 或者 64），代码会相对简单一些。

#include<bits/stdc++.h>
using namespace std;
#define LL long long 
#define i128 __int128
mt19937_64 rd(time(0));

LL n;

i128 ksm(i128 x, i128 y, LL mod)
{
    i128 res = 1;
    x %= mod;
    while(y)
    {
        if(y & 1) res = res * x % mod;
        y >>= 1;
        x = x * x % mod;
    }
    return res;
}

bool miller_rabin(LL n)
{
    if(n <= 1 || n > 2 && ~n & 1) return false;
    LL m = n - 1, k = 0;
    while(~m & 1) m >>= 1, k++;
    vector<int> s = {2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37};
    for(int a : s)
    {
        if(a >= n) break;
        i128 x = ksm(a, m, n), y;
        for(int i = 1; i <= k; i++, x = y)
        {
            y = x * x % n;
            if(y == 1 && x != 1 && x != n - 1) return false;
        }
        if(x != 1) return false;
    }
    return true;
}

vector<LL> p;
void pollard_rho(LL n)
{
    if(n <= 1) return;
    if(miller_rabin(n))
    {
        p.push_back(n);
        return;
    }
    LL d = rd() % (n - 1) + 1;
    auto f = [&](LL x)
    {
        return ((i128)x * x + d) % n; 
    };
    LL x = rd() % (n - 1) + 1, y = f(x), t = 1;
    while(t == 1)
    {
        for(int i = 1; i <= 32; i++)
        {
            while(x == y)
                d = rd() % (n - 1) + 1, x = rd() % (n - 1) + 1, y = f(x);
            if((i128)t * abs(x - y) % n == 0) break;
            t = (i128)t * abs(x - y) % n;
            x = f(x), y = f(f(y));
        }
        t = __gcd(n, t);
    }
    pollard_rho(t), pollard_rho(n / t);
}

int main()
{
    int T;
    cin >> T;
    while(T--)
    {
        cin >> n;
        if(miller_rabin(n))
        {
            puts("Prime");
            continue;
        }

        p.clear();
        pollard_rho(n);
        LL ans = 0;
        for(LL i : p) ans = max(ans, i);
        cout << ans << endl;
    }
    return 0;
}