Miller-Rabin 与 Pollard-Rho 算法学习笔记

· · 题解

Miller-Rabin 与 Pollard-Rho 算法学习笔记

也许有更好的阅读体验

前言

Miller-Rabin 算法用于判断一个数 p 是否是质数，若选定 w 个数进行判断，那么正确率约是 1-\frac{1}{4^w} ，时间复杂度为 O(\log p+w\log p)。（我的实现）

Pollard-Rho 算法可以在期望 O(n^{\frac{1}{4}}) 的时间复杂度内找到合数 n 的某一个非平凡的（即既不是 1，也不是它本身的）因子。

下文中用 \mathbb{P} 来表示质数集合。

Miller-Rabin 算法

前置知识

费马小定理：若 p\in\mathbb{P},\gcd(a,p)=1，则 a^{p-1}\equiv1\pmod{p}。

二次探测定理：若 p\in\mathbb{P},x^2\equiv 1\pmod{p}，则 x\equiv\pm1\pmod{p}。

注意：费马小定理的逆命题并不成立！

算法流程

首先，将 p-1 表示成 t2^k 的形式。那么，若 p\in\mathbb{P}，则 p^{t2^k}=p^{p-1}\equiv1\pmod{p}。

然后我们选择 w 个数 q_1,q_2,\cdots,q_w 进行判断。假如当前判断到了 q_i，那么用快速幂计算出 a=q_{i}^{t}\bmod{p}。然后让 a 自乘 k 次，就可以得到 p-1。

自乘的时候我们判断，如果 a\equiv1\pmod{p} ，那么此时 p 有一定概率是质数。于是我们看一看 a 自乘前是否满足二次探测定理即可，如果是，则继续自乘，否则表明 p 一定不是质数。

如果自乘得到的数同余 p 不为 1，那么 p 也不一定是质数。否则 p 很可能是合数。

这时您可能会说：不是说过费马小定理逆定理不成立吗？其实，逆定理的反例（卡迈克尔数）是十分稀少的。经过多次判断，合数判成质数的概率十分小（质数不可能判成合数，想一想，为什么）

OI 中可以选择 w=9，q 为前 9 个质数。这样子 10^{18} 范围内一般不会出错。

参考实现

struct {/*Miller Rabin 质数判定算法*/
    vector<int> primes= {2,3,5,7,11,13,17,19,23};
    bool operator()(int p) {
        if (p==1)return 0;
        int t,k;
        for (t=p-1,k=0; !(t&1); k++,t>>=1);
        for (int i : primes) {
            if (p==i) return true;
            int a=fastpow(i,t,p),b=a;
            for (int j=1; j<=k; j++) {
                b=M(((__int128)a)*a,p);
                if (b==1&&a!=1&&a!=p-1) return false;
                a=b;
            }
            if (a!=1) return false;
        }
        return true;
    }
} MillerRabin;

模板题

Pollard-Rho 算法

前置知识

Floyd 判圈算法：该算法可以线性判断一个链表上是否有环。其流程为使用两个指针。一个指针每次跑 1 条边，另一个指针一次跑 2 条边，然后相遇的点就在环上。

算法流程

先特判 n=4 和 n\in\mathbb{P}。

Pollard-Rho 需要一个伪随机函数 f(x,c,n)=x^2+c\bmod{n}。其中 x 表示上一个数，c 是我们生成的，用于保证随机性的数，n 是我们需要找因子的数。

可以发现这个函数最后会大概率生成一个混循环序列。如同希腊字母 \rho(\texttt{Rho}) 一般。

先选择一个随机数 c。两个指针从 0 出发，我们看成存在一个链表，其中存在边 (i,f(i,c,n))。然后在上面跑 Floyd 判圈算法，在 Floyd 中，如果一个指针在 t，一个指针在 k。若 \gcd(|t-k|,n)\neq1。则我们认为 \gcd(|t-k|,n) 是 n 的一个因数。如果找到了环，则重新选择一个 c，重复上述流程。

此时时间复杂度期望 O(n^{\frac{1}{4}}\log n)。

算法优化

上述算法在洛谷板题上只能获得 93 分（TLE 在了第 13 个点）。优化迫在眉睫。

我们发现求 \gcd 的 O(\log n) 需要被优化。我们可以固定一个 W，跳 W 次的时候统计 |t-k| 的乘积 p。最后和 n 取一次 \gcd。然后如果下一次 p 会到 0，那么也要跳出。因为后面都是 0。

这样子只要 W\gt \log n 就可以做到期望 O(n^{\frac{1}{4}})。我试了一下，貌似 W=256 表现不错。

另外还可以加一个记忆化，在后面的例题中有用。

参考实现

mt19937 engine(time(0));

inline int pr_rand(int x,int c,int n) { /*Pollard Rho 算法使用的伪随机数*/
    return M(M(((__int128)x)*x,n)-c,n);
}

unordered_map<int,int> prm;

int pollard_rho(int n) { /*Pollard Rho 算法求一个数的某一个质因子*/
    if (prm[n])return prm[n];
    if (n==4) return 2;
    if (MillerRabin(n)) return n;
    uniform_int_distribution<int> randint(3,n-1);
    while (1) {
        int c=randint(engine);
        int t=0,r=0,p=1,q=0;
        do{
            for(int i=1;i<=256;i++){
                t=pr_rand(t,c,n);
                r=pr_rand(pr_rand(r,c,n),c,n);
                int delta=(t-r)>0?(t-r):(r-t);
                if(t==r||(q=M(__int128(p)*delta,n))==0){
                    break;
                }
                p=q;
            }
            int d=__gcd(p,n);
            if(d>1) return prm[n]=d;
        }while(t!=r);
    }
}

P4718 【模板】Pollard's rho算法

简要题意

$1 \leq T \leq 350,1 \lt n \leq 10^{18}

思路

首先，我们先用一个 Miller-Rabin 算法来判断质数。我们可以用 Pollard-Rho 算法找出所有的因子（当然不用存起来，只需要用一个递归函数，最后 \max 统计答案）即可。由于唯一分解定理，这个算法是正确的。

注意我们需要加一个记忆化来保证复杂度，否则复杂度爆炸。

代码

关键代码如下（要完整代码的私信）：

unordered_map<int,int> mrm;

int max_factor(int n) { /*求一个数的最大质因子*/
    if (mrm[n]) return mrm[n];
    int factor=pollard_rho(n);
    if (factor==n) return mrm[n]=n;
    return mrm[n]=max(max_factor(factor),max_factor(n/factor));
}

AC Record

参考资料

Miller Rabin 算法详解 - 自为风月马前卒 - 博客园

算法学习笔记(55): Pollard-Rho 算法 - 知乎