专栏文章

笔记 · 素性判断等（Miller-Rabin，Pollard-Rho）

aaleph_
个人记录参与者 1已保存评论 0

文章操作

快速查看文章及其快照的属性，并进行相关操作。

当前评论
0 条
当前快照
1 份
快照标识符
@mio5ygdu
此快照首次捕获于
2025/12/02 13:54
3 个月前
此快照最后确认于
2025/12/02 13:54
3 个月前
查看原文

Miller-Rabin

两个引理

引理 11：费马小定理

设 p∈Primep\in\text{Prime}，a∈Za\in\mathbb{Z}，且 gcd⁡(a,p)=1\gcd(a,p)=1，那么有
ap−1≡1(modp)a^{p-1}\equiv 1\pmod p
证明：
对于两个整数 a≠b(modp)a\ne b\pmod p，和一个和 pp 互素的整数 cc，则有 ac≠bc(modp)ac\ne bc\pmod p.
考虑两组数：1,2,…,p−11,2,\dots,p-1；a,2a,…,(p−1)aa,2a,\dots,(p-1)a。根据上述结论，两者均为模 pp 的完全剩余系。
因此：
1⋅2⋅3⋯(p−1)≡a⋅2a⋅3a⋯(p−1)a(modp)1\cdot2\cdot3\cdots(p-1)\equiv a\cdot 2a\cdot 3a\cdots (p-1)a\pmod p
(p−1)!≡(p−1)!⋅ap−1mod  p(p-1)!\equiv (p-1)!\cdot a^{p-1}\mod p
∵gcd⁡((p−1)!,p)=1\because \gcd((p-1)!,p)=1
∴ap−1≡1(modp)\therefore a^{p-1}\equiv 1\pmod p

引理 22：二次探测定理

如果 p∈Primep\in\text{Prime}，且 0<x<p0<x<p，则方程 x2≡1(modp)x^2\equiv 1\pmod p 的解为 x=1,p−1x=1,p-1.
证明：
易知 x2−1≡0(modp)x^2-1\equiv 0\pmod p
∴(x+1)(x−1)≡0(modp)\therefore (x+1)(x-1)\equiv 0\pmod p
∴p∣(x+1)(x−1)\therefore p|(x+1)(x-1)
∵p∈Prime,0<x<p\because p\in\text{Prime},0<x<p
∴x=1,p−1\therefore x=1,p-1

算法流程

费马小定理的缺陷是，有一些特殊的整数，虽然 an−1≡1(modn)a^{n-1}\equiv 1\pmod n，但 nn 是合数。我们可以考虑通过随机化的方式，来规避这种数，这也是 Miller-Rabin 算法的核心思想。
  • 对于 0,1,20,1,2 和偶数，可以直接判断。
  • 设要测试的数是 NN。我们取一个较小的质数 aa。设 N−1=2s⋅tN-1=2^s\cdot t（其中 tt 是奇数）。这有什么用呢？后面会说到。
  • 对于 x=atx=a^t，我们不断去平方他，进行二次探测检验。根据二次探测定理，如果 x2≡1(modN)x^2\equiv 1\pmod N，但 x≠1x\ne 1 且 x≠N−1x\ne N-1，那么 NN 就一定不是质数，退出。
  • 我们最多可以平方 ss 次，显然这期间 xx 满足 0<x<N0<x<N。最后 x=a2s⋅t=aN−1x=a^{2^s\cdot t}=a^{N-1}，再使用费马小定理进行验证。如果 aN−1≢1(modN)a^{N-1}\not\equiv 1\pmod N，那么 NN 肯定不是质数，退出。
  • 注意，即使通过了上述检验，NN 也可能不是质数。那么我们就多取几个质数 aa，确保了检验的准确性。

注意事项

  • aa 取遍小于等于 3030 的所有质数，就可以保证 int 范围内的数不会出错。
  • 注意求平方的过程中，可能会爆 long long，需要用到龟速乘。怕超时的话，直接开 __int128
  • 对于单个数 nn，复杂度为 O(klog⁡3n)\mathcal{O}(k\log^3 n)，其中 kk 为检验次数，此处 k=10k=10。可以保证结果大概率正确。

代码

CPP
#include <bits/stdc++.h>
#define int long long
using namespace std;
int prime[10]={2,3,5,7,11,13,17,19,23,29};
int qmul(int a,int b,int mod){
    int ans=0;
    while(b){
        if(b&1) ans=(ans+a)%mod;
        a=(a+a)%mod;
        b>>=1;
    }
    return ans;
}
int qpow(int a,int b,int mod){
    int ans=1;
    while(b){
        if(b&1) ans=ans*a%mod;
        a=a*a%mod;
        b>>=1;
    }
    return ans;
}
bool Miller_Rabin(int n){
    if(n==2) return 1;
    if(n<2||!(n&1)) return 0;
    //特判0,1,2和偶数
    int s=0,t=n-1;
    while(!(t&1)){
        s++;
        t>>=1;
    }
    for(int i=0;i<10&&prime[i]<n;i++){
        int a=prime[i];
        int x=qpow(a,t,n);
        for(int j=1;j<=s;j++){
            int y=qmul(x,x,n);
            if(y==1&&x!=1&&x!=n-1) return 0;//二次探测检验
            x=y;
        }
        if(x!=1) return 0;//此时x=a^{n-1}，用费马小定理
    }
    return 1;
}
signed main(){
    int n;
    cin>>n;
    cout<<(Miller_Rabin(n)?"YES\n":"NO\n");
}

Pollard-Rho

算法概述

该算法能快速找到大整数的一个非1、非自身的因子。
原理这篇文章都讲透了
首先，可以考虑生成 [1,n][1,n] 的随机数 xx，计算 d=gcd⁡(x,n)d=\gcd(x,n)，若 d>1d>1，那么就找到了，dd 就是一个 nn 的一个因子。
最坏情况下，n=p2n=p^2，pp 是一个素数。那么当 x=p,2p,3p,⋯ ,(p−1)px=p,2p,3p,\cdots,(p-1)p 时有 d>1d>1，一共有 p−1≈np-1\approx \sqrt{n} 个这样的 xx，所以最差期望时间复杂度是 O(nlog⁡n)\mathcal{O}(\sqrt{n}\log n) 的。（log是算gcd带来的）

生日悖论

对于一个 [1,N][1,N] 内整数的理想随机数生成器，生成序列中第一个重复的数字前期望有 πN2\sqrt{\frac{\pi N}{2}} 个数。(生日悖论)
对于 n=p2n=p^2，在 [1,n−1][1,n-1] 间生成的随机数模 pp，大致是 [0,p−1][0,p-1] 间的随机数。因此期望在生成大约 p=n14\sqrt{p}=n^{\frac{1}{4}} 个数后，可以出现两个在模 pp 意义下相同的数。那么这两个数的差 dd，就满足 d≡0(modp)d\equiv 0\pmod p，也就满足 gcd⁡(d,n)>1\gcd(d,n)>1，这时就找到了。但问题在于，我们没法两两进行比较验证，否则时间复杂度就会退回，还是没有进步。

伪随机数生成器

接下来是 Pollard-Rho 算法的精髓——使用一种特别的伪随机数生成器来生成 [0,n−1][0,n-1] 内的伪随机数序列。
设序列第一个数是 xx，f(x)=(x2+c)mod  nf(x)=(x^2+c)\mod n
则 x,f(x),f(f(x)),…x,f(x),f(f(x)),\dots 是生成的伪随机数序列。
观察这个伪随机数生成器，刚开始的时候是随机的，到后面就会进入一个环。因为下一个数是依赖于上一个数的，而可生成的数又是有限的，因此会有环。所以这个伪随机数生成器生成的数，可以看作是一条链，最终连向一个环。看起来就像字母 ρ\rho。
我们可以使用龟兔赛跑算法（Floyd判环法）。
  • 设置两个变量 t,rt,r。每次判断是否有 ∣t−r∣≡0(modn)|t-r|\equiv 0\pmod n，即判断 gcd⁡(∣t−r∣,n)>1\gcd(|t-r|,n)>1。有的话，我们就找到了。这个期望显然是 n14n^{\frac{1}{4}}。
  • 如果没有，那么令 t=f(t),r=f(f(r))t=f(t),r=f(f(r))，rr 跑的快，所以如果最后没有找到答案，最终会在环上追上 tt。这时退出，换一个 cc 重新生成伪随机数。
注意到这个伪随机数生成器有一个性质。
对于 ∣i−j∣≡0(modp)|i-j|\equiv 0\pmod p，有
∣f(i)−f(j)∣=∣i2−j2∣=∣(i−j)(i+j)∣≡0(modp)|f(i)-f(j)|=|i^2-j^2|=|(i-j)(i+j)|\equiv 0\pmod p。
如果 i,ji,j 是环上满足条件的两个数，它们之间的距离为 dd（两者之间间隔的边数），那么它们的后继 f(i),f(j)f(i),f(j) 的距离也为 dd，也满足条件。同理，f(f(i)),f(f(j))f(f(i)),f(f(j)) 也满足条件……
也就是说，如果环上两个相距为 dd 的数满足条件，那么环上所有相距为 dd 的数都满足条件。而在Floyd判环的过程中，由于两个指针的移动速度不同，每次移动都是在检查一个新的距离 dd，而这样就无需两两比较了。
这个算法的复杂度依赖于该伪随机数生成器的随机程度。如果它是足够随机的，那么期望时间复杂度就是 O(n14log⁡n)\mathcal{O}(n^{\frac{1}{4}}\log n)。

分块优化

这个 log⁡\log 还是有点讨厌。怎么去除呢？我们发现瓶颈在于计算 gcd⁡\gcd。
注意到若 gcd⁡(d,n)>1\gcd(d,n)>1，那么一定有 gcd⁡(kdmod  n,n)>1\gcd(kd\mod n,n)>1，所以我们可以先把一些待选数乘起来，再统一和 nn 求公因数，这样就减少了求 gcd⁡\gcd 的次数。
所以我们可以每隔 1,2,4,8,…1,2,4,8,\dots 步求一次 gcd⁡\gcd。但这样到了后面就会走越来越多步才能退出，如果早就找到答案却迟迟不退出，显然不优，因此可以在倍增取距离的基础上，设定一个固定步数上限 CC。当然也可以简单地每隔 CC 步就计算一下。这样时间复杂度至此优化至 O(n14+n14log⁡nC)\mathcal{O}(n^{\frac{1}{4}}+\frac{n^\frac{1}{4}\log n}{C})。当 C>log⁡nC>\log n 时，时间复杂度为 O(n14)\mathcal{O}(n^{\frac{1}{4}})。

注意事项：

  • 为了方便，通常令 x,cx,c 任取。
  • 当 n=4n=4 时，无论 cc 取多少，都得不到解，需要特判。质数也直接特判返回。
  • 一些细节和优化见代码。

例题

模版：P4718
题意：q(1≤q≤350)q(1\le q\le 350) 次询问，每次给出一个正整数 n(2≤n≤1018)n(2\le n\le10^{18})，如果它是质数输出 Prime，否则输出它最大的质因数。

Solution：

基本上就是 Pollar-Rho 和 Miller-Rabin 板子。注意，本题中用龟速乘会超时，只能强转 __int128。
还需注意本题中要求找最大质因数。我们可以递归地解决该问题。对于当前数 nn，Pollard-Rho 算法每次会找到一个 nn 的一个非1、非自身因数 dd。我们可以递归地去解决数为 dd 和 nd\frac{n}{d} 的子问题。注意往下递归之前判一下 dd 或者 nd\frac{n}{d} 是否已经是质数了，因为Pollard-Rho算法规避找到自己的方法是“多试几次”，如果本来就是质数的话，就会一直循环循环……
容易证明，递归的底层一定会得到原数的质因数，因为如果当前数是合数就能继续向下递归。我们取最大值即可。

Code:

CPP
#include <bits/stdc++.h>
#define int long long
using namespace std;
const __int128 one=1;
int prime[12]={2,3,5,7,11,13,17,19,23,29,31,37};//n拓展到10^18之后，k=10不太够了，这里取k=12
int qpow(int a,int b,int mod){
    int ans=1ll;
    while(b){
        if(b&1) ans=one*ans*a%mod;
        a=one*a*a%mod;
        b>>=1;
    }
    return ans;
}
bool Miller_Rabin(int n){
    if(n==2) return 1;
    if(n<2||!(n&1)) return 0;
    int s=0,t=n-1;
    while(!(t&1)){
        s++;
        t>>=1;
    }
    for(int i=0;i<12&&prime[i]<n;i++){
        int a=prime[i];
        int x=qpow(a,t,n);
        for(int j=1;j<=s;j++){
            int y=one*x*x%n;
            if(y==1&&x!=1&&x!=n-1) return 0;
            x=y;
        }
        if(x!=1) return 0;
    }
    return 1;
}
const int C=128;
static std::mt19937 aleph;
int Pollard_Rho(int n){
    if(n==4) return 2;
    if(Miller_Rabin(n)) return n;
    std::uniform_int_distribution<int> rnd(3,n-1);
    int t=rnd(aleph),r=t,c=rnd(aleph);
    auto f=[=](int x){return (one*x*x%n+c)%n;};
    t=f(t),r=f(f(r));
    for(int lim=1;t!=r;lim=min(lim<<1,C)){
        int fac=1ll;
        for(int i=0;i<lim;i++){
            fac=one*fac*abs(t-r)%n;
            if(!fac) break;//小优化：如果当前乘积是0，后面也没必要乘了，怎么乘还是0，直接退出。
            t=f(t),r=f(f(r));
        }
        int d=gcd(fac,n);
        if(d>1) return d;
    }
    return n;
}
int ans=0;
void mx(int z){
    if(z>ans) ans=z;//优化，比max快
}
void dfs(int n){
    int d=Pollard_Rho(n);
    while(d==n) d=Pollard_Rho(n);//Pollard-Rho会随机找到一个因数，如果刚好找到了n，那就再做几次
    int d2=n/d;
    if(Miller_Rabin(d)) mx(d);//递归下去前判一下质数，不然死循环
    else dfs(d);
    if(Miller_Rabin(d2)) mx(d2);
    else dfs(d2);
}
int calc(int x){
    if(Miller_Rabin(x)) return x;
    ans=0;
    dfs(x);
    return ans;
}
void solve(){
    int n;
    scanf("%lld",&n);
    int x=calc(n);
    if(x==n) printf("Prime\n");
    else printf("%lld\n",x);
}
signed main(){
    srand(time(0));
    int T;
    cin>>T;
    while(T--){
        solve();
    }
}

评论

共 0 条评论，欢迎与作者交流。

正在加载评论...