关于Miller-Rabin与Pollard-Rho算法的理解(素性测试与质因数分解)
阅读原文时间:2023年07月10日阅读:1

前置

费马小定理(即若P为质数,则\(A^P\equiv A \pmod{P}\))。

欧几里得算法(GCD)。

快速幂,龟速乘。

素性测试

素性测试是OI中一个十分重要的事,在数学毒瘤题中有着举足轻重的地位。

常见的素性测试如下:

int check(int N){
    for(int i=2;i*i<=N;i++)
        if(N%i==0)return 0;
    return 1;
}

以上是一个\(O(\sqrt{N})\)的算法,虽然不优,但在绝大多数情况下是可以的。

但是,假如\(N\)的范围达到了\(1e18\),以上算法很明显就不行了,我们得考虑更优的算法。

引入Miller-Rabin算法。

定理1:如果\(P\)为一个大于2的素数,那么方程\(X^2\equiv1 \pmod{P}\)的解只有1或者-1。

(是真得水的) 证明如下:

由\(X^2\equiv1 \pmod{P}\),得\((X^2-1)\equiv0 \pmod{P}\),

即\({(X+1)(X-1)}\equiv0 \pmod{P}\)。

则\(P \mid (X+1)\)或\(P\mid(X-1)\),

即\((X+1)\equiv0 \pmod{P}\)或\((X-1)\equiv0 \pmod{P}\)。

则\(X\equiv-1 \pmod{P}\)或\(X\equiv1 \pmod{P}\),即证。

(注:由于P>2,不能同时满足\(X\equiv-1 \pmod{P}\)与\(X\equiv1 \pmod{P}\)。)

定理1的逆否定理即为:

如果有\(P>2\)且\(X^2\equiv1 \pmod{P}\)且\(X\not\in\{1,-1\}\),则\(P\)不为质数。

证明如下:

假设\(P\)为质数,则\(P\)满足定理1,即方程\(X^2\equiv1 \pmod{P}\)的解只有1或者-1。

与条件矛盾,故\(P\)不为质数。

基于以上定理,我们可以开始了。

假如我们要验证\(P(P>2)\)是否是一个质数。

先设\(P\)是一个素数。

那么对于所有满足\(A^2\equiv1 \pmod{P}\)的整数\(A\),都有\(A\in\{1,-1\}\)。

但如果我们枚举\(A\)来暴力判断,则复杂度又上升到了\(O(N)\)。

(不妨令\(A<P\),以下\(A\)均小于\(P\))

考虑做判断时,如何快速的找到满足\(A^2\equiv1 \pmod{P}\)的整数\(A\)。

我们考虑将\((P-1)\)分解成\(2^S·D\)的形式,其中\(D\)为奇数。

(由于\(P>2\)且\(P\)为素数,则\((P-1)\)为偶数,即\(S\ge1\))

则根据费马小定理,有

\[A^{P-1}\equiv 1\pmod{P}
\]

\[A^{2^S·D}\equiv 1\pmod{P}
\]

\[({A^{2^{S-1}·D}})^2\equiv 1\pmod{P}
\]

则此时\(A^{2^{S-1}·D}\)又变成了在\(A^2\equiv1 \pmod{P}\)的新的\(A\)值。

不妨设\(A\)值数组为\(\{A_0,A_2,A_3,A_4….A_{S-1}\}\),其中\(A_i=A^{2^i·D}\)。

(注:不难发现对于任意满足\(0<i<S\)的\(i\),都有\(A_i=2·A_{i-1}\))

这样的话对于每个初始值\(A\),我们可以知道所构成的\(A\)数组值。

此时我们可以考虑从\(A_0\)走到\(A_{S-1}\)(即从\(A^D\)到\(A^{2^{S-1}·D}\))。

①考虑初始情况的\(A_0\)(即为\(A^D\)),

它如果满足\(A_0\equiv1 \pmod{P}\)或\(A_0\equiv-1 \pmod{P}\)。

则对于任意\(A_i\)(\(0<i< S\)),都有\(A_i\equiv1 \pmod{P}\),即通过本次测试。

(注:通过测试只能说明它可能为素数,而没通过则说明它一定不是素数。)

②考虑过程中的情况:设有一满足\(0<t<S\)的\(t\)。

  1. 若有\(A_t\equiv-1\pmod{P}\),则说明

    对于任意满足\(t<i<S\)的\(i\),都有\(A_i\equiv1\pmod{P}\),即通过此次测试。

  2. 若有\(A_t\equiv1\pmod{P}\),则说明\({A_{t-1}}^2\equiv1\pmod{P}\),

    而\(A_{t-1}\)又不满足\(A_{t-1}\equiv-1\pmod{P}\)或\(A_{t-1}\equiv-1\pmod{P}\),

    因为如果\(A_{t-1}\)满足以上条件,则在上一次或\(A_0\)时就已经判了。

    故存在一个\(A\),使得\(A^2\equiv1\pmod{P}\)且\(A\not\in\{1,-1\}\)(定理1的延伸)。

    故未通过此次测试。

  • 注:若\(P(P>2)\)为素数: 则\(A^2\equiv1\pmod{P}\)只有\(A\)满足

    \(A\equiv1\pmod{P}\)或\(A\equiv-1\pmod{P}\)时满足(由定理1得)。

    但\(A^2\equiv-1\pmod{P}\)对于\(A\)的取值却没有特殊性

③考虑结尾情况:

若对于这一组\(A\)都没有一个\(A\)使得\(A^2\equiv1\pmod{P}\),

则在判\(A_{S-1}\)时,实际上判的是\({({A^{2^{S-1}·D}})}^2\not\equiv1\pmod{P}\)

即\(A^{P-1}\not\equiv1\pmod{P}\),故\(P\)不满足费马小定理,即\(P\)未通过测试。

为保证算法的正确性与高效性,

我们考虑随机化许多个整数\(A\)来多次判断使得\(P\)为素数的概率尽量大。

那么我们要随机多少次才能使得\(P\)素数的概率能让我们接受呢?

根据大佬们的研究表明,Prime是一个RP问题(通过大量数据来算概率),

而对于每次测试的期望失败概率是比\(\frac{1}{4}\)要小的。(通过试验得到的结论)

即我们做K次成功的测试然后判错\(P\)的概率即为\(4^{-K}\)。

这个失败率在代码中也是可以接受的。

(注:经过试验,当A值取\(\{2,3,5,7,11,13,17,19,23,29\}\)时,可以过long long级别的询问)

int MillerTest(LL N){//素数测试
    if(N<=1||N==4)return 0;//特判N比较小情况
    if(N<=3)return 1;
    LL D=N-1;
    while(D%2==0)D/=2;
    for(int i=1;i<=PK;i++){
        if(Prime[i]>=N)return 1;
        int Ok=0;
        LL A=Prime[i];
        LL x=quick_Pow(A,D,N);
        if(x==1||x==N-1)
            continue;
        while(D!=N-1){
            x=quick_Mul(x,x,N);//龟速乘防爆long long
            D*=2;
            if(x==1)return 0;
            if(x==N-1){Ok=1;break;}
        }
        if(!Ok)return 0;
    }
    return 1;
}

质因数分解

质因数分解是OI中一个十分重要的事,在数学毒瘤题中有着举足轻重的地位。

常见的质因数分解如下:

void divide(int N){
    for(int i=2;i*i<=N;i++){
        int cnt=0;
        while(N%i==0)cnt++,N/=i;
        if(cnt)printf("%d %d\n",i,cnt);
    }
    if(N!=1)printf("%d 1\n",N);
}

以上是一个\(O(\sqrt{N})\)的算法,虽然不优,但在绝大多数情况下是可以的。

但是,假如\(N\)的范围达到了\(1e18\),以上算法很明显就不行了,我们得考虑更优的算法。

引入Pollard-Rho算法。

说起Pollard-Rho算法,我们可能不太熟悉。

但提到随机化,大家应该多少都有些了解。

//随机化找N的一个因数
int i=0;do{
    i=rand()%(N-2)+2;
}while(N%i);
printf("I found %d \n",i);

我们不难发现:随机化一次找到某个因数的概率十分小。

当\(N=P·Q\)时(\(P,Q\)均为素数),成功的概率为\(\frac{2}{N-1}\)。

期望次数为\(O(N)\)级别,(甚至不如暴力)。

但是,如果我们加上一些优化呢?

考虑每次随机化。

其实对于随机的某个数\(i\),只要\(gcd(i,N)\not=1\)就可以找到一个因数。

(时间复杂度仿佛还是\(O(N)\)….)

我们考虑用生日悖论来优化随机化。

尽管听起来特别玄学,但是此做法还是有一定科学依据的。

生日悖论:

生日悖论,指如果一个房间里有23个或23个以上的人,那么至少有两个人的生日相同的概率要大于50%。这就意味着在一个典型的标准小学班级(30人)中,存在两人生日相同的可能性更高。对于60或者更多的人,这种概率要大于99%。

以上是生日悖论的一个例子,其函数图像大致如下。

考虑换一种思考方式,我们随机化\(N\)个数出来,得到\(N^2\)个差值。

然后我们用这些差值来执行优化1,时间复杂度将明显降低。

考虑生日相同的情况,其可以被描述为某两个人的生日值差值为0。

从以上描述差值为0的图像中,可以明显发现选取的人数\(N\)越大,则概率越大。

则可以发现\(N\)的升大对于所有\(K\)的概率均有提升。

(易发现\(N^2\)对于每个\(K\)的期望值为总值域\(W\)的根号\(O(\sqrt{W})\))

然后,我们使\(N\)值固定。

我们设\(P(K)\)表示以\(K\)为差值的概率。

则感性发现:\(P(1)>P(0),P(W-1)>P(W)\)

感性推出:实际的\(P\)图像其实是在前方激增,中间接近1,后方再次猛降。

故我们在判断时,得特判差值出现在前面的情况。

(注:\(N\)越大时,\(P\)图像的前方激增斜率便越大。)

应用优化1与优化2,考虑如何找到\(N^2\)个差值,

(可以认为生日悖论中的某个差值\(K\)是我们所需要的答案)

我们考虑随机出两个数\(X\),\(Y\),然后用\(abs(X-Y)\)作为差值。

然后令\(f(X)=(X^2+P) \mod W,X=f(X),Y=f(Y)\)。

这样的话,我们只用随机化两个初值,一个常数\(P\),就可以了。

但是易发现这样可能重复,即走进一个无解的循环里。

所以我们让\(X,Y\)的初值相等,每次随机的时候令\(Y=f(f(Y))\),

即出发点一样,\(Y\)的速度为\(X\)的两倍,若有无解循环,

则总有一次\(X,Y\)会相遇(复杂度至多为经过的点数量),此时我们就可以直接退出了。

(注:\(f(X)=(X^2+P) \mod W\)是一个先人传下的神奇函数)

(其神奇之处在于玄学般地大幅提升了代码效率)

其实随机这些差值的本质还是随机化,依旧是\(N^2\)个差值,但复杂度却玄学降低了。

???

故这个算法的复杂度十分看RP,

在一些数据的复杂度甚至超过了\(O(\sqrt{W}\)),

不过能解决一些long long的数据还是很强大了。

结合Miller-Rabin与Pollard-Rho算法。

我们对于一个整数\(N\),先用Miller-Rabin判断其是否为质数。

若是,则上传N。

若不是,则用Pollard-Rho找一个因数出来,

将其分解成更小的两个数继续递推。

//该代码只是找到了N的一个因数,可能并不是素数
LL PollardRho(LL N){//找因数
    if(N==1)return -1;
    for(int i=1;i<=PK;i++)
        if(N%Prime[i]==0)
            return Prime[i];//特判出奇迹
    LL x=(rand()%(N-2))+2,y=x;
    LL P=(rand()%(N-1))+1,D=1;
    while(D==1){
        x=(quick_Mul(x,x,N)+P)%N;
        y=(quick_Mul(y,y,N)+P)%N,y=(quick_Mul(y,y,N)+P)%N;
        D=gcd(abs(x-y),N);//x=y时,D值为N
    }
    return D;
}

例题及代码

板题:Prime Test

大意:找到\(N\)的最小的质因数。

#include<ctime>
#include<cstdio>
#include<algorithm>
using namespace std;
#define LL long long
const int MAXK=15;
const long long INF=1e18;
int PK=10,Prime[MAXK]={0,2,3,5,7,11,13,17,19,23,29};
//此行数对于long long范围内的素性测试已够用了
long long N;
LL Add(LL x,LL y,LL p){
    if(x+y<p)return x+y;
    return x+y-p;
}
LL gcd(LL x,LL y){//最大公因数
    if(x==0)return y;
    return gcd(y%x,x);
}
LL quick_Mul(LL x,LL y,LL p){//龟速乘
    x%=p;y%=p;
    LL ret=0;
    while(y){
        if(y&1)ret=Add(ret,x,p);
        x=Add(x,x,p);
        y>>=1;
    }
    return ret;
}
LL quick_Pow(LL x,LL y,LL p){//快速幂
    LL ret=1;
    while(y){
        if(y&1)ret=quick_Mul(ret,x,p);
        x=quick_Mul(x,x,p);
        y>>=1;
    }
    return ret;
}
int MillerTest(LL N){//素数测试
    if(N<=1||N==4)return 0;
    if(N<=3)return 1;
    LL D=N-1;
    while(D%2==0)D/=2;
    for(int i=1;i<=PK;i++){
        if(Prime[i]>=N)return 1;
        int Ok=0;
        LL A=Prime[i];
        LL x=quick_Pow(A,D,N);
        if(x==1||x==N-1)
            continue;
        while(D!=N-1){
            x=quick_Mul(x,x,N);
            D*=2;
            if(x==1)return 0;
            if(x==N-1){Ok=1;break;}
        }
        if(!Ok)return 0;
    }
    return 1;
}
LL PollardRho(LL N){//找因数
    if(N==1)return -1;
    for(int i=1;i<=PK;i++)
        if(N%Prime[i]==0)
            return Prime[i];//特判出奇迹
    LL x=(rand()%(N-2))+2,y=x;
    LL P=(rand()%(N-1))+1,D=1;
    while(D==1){
        x=(quick_Mul(x,x,N)+P)%N;
        y=(quick_Mul(y,y,N)+P)%N,y=(quick_Mul(y,y,N)+P)%N;
        D=gcd(abs(x-y),N);
    }
    return D;
}
LL Min_Div(LL N){//最小质因数
    if(N==1)return INF;
    if(MillerTest(N))return N;
    LL X=PollardRho(N);
    return min(Min_Div(X),Min_Div(N/X));
}
int main(){
    //srand(time(NULL));
    while(~scanf("%lld",&N)){
        printf("%lld\n",Min_Div(N));
    }
}

手机扫一扫

移动阅读更方便

阿里云服务器
腾讯云服务器
七牛云服务器

你可能感兴趣的文章