米勒来宾素数测试与Pollard rho算法

mod998244353 · 2021-03-11 17:59:40 · 个人记录

可能更好的阅读体验

米勒来宾素数测试

此处只考虑如何判断奇素数

根据费马小定理：

当p为（奇）素数，\gcd(a,p)=1时，a^{p-1}\equiv 1\pmod p

就可以写一个用快速幂判断一个数是不是奇素数

#include<bits/stdc++.h>
using namespace std;
const int plist[]= {2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71,73,79,83,89,97};
//素数表
const int L=25;
int fastpow(int a,int k,int mod) {
    int base=1;
    while(k) {
        if(k&1) base=base*1ll*a%mod;
        a=a*1ll*a%mod;
        k>>=1;
    }
    return base;
}
bool isprime(int p) {
    if(p<=100) {//在表内就直接判断
        for(int i=0; i<L; ++i) {
            if(plist[i]==p)return true;
        }
        return false;
    }
    for(int i=0; i<L; ++i) {//否则看看是否符合费马小定理
        if(fastpow(plist[i],p-1,p)!=1)return false;
    }
    return true;
}
int n;
int main() {
    scanf("%d",&n);
    printf(isprime(n)?"YES":"NO");
    return 0;
}

但这个效率低，而且错误的概率很大。

我们可以发现，当p为奇素数时，p-1为偶数，就可以写成2^r\times m的形式（m为奇数）。

当a^m\equiv 1\pmod p时，一定通过测试。

否则当p为奇素数时，必定有一个数k(0\leq k<r)使a^{2^k\times m}\equiv p-1\pmod p，这样再进行平方（指数乘二）后才能是1。

所以我们可以用k从0枚举到r-1，若有一个k使a^{2^k\times m}\equiv p-1\pmod p则通过测试，否则p不为素数。

同时计算a^{2^k\times m}不用快速幂，直接用(a^{2^{k-1}\times m})^2计算即可

const int plist[]= {2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71,73,79,83,89,97};
//素数表
const int L=25;
inline ll mul(const ll&a,const ll&b,const ll&mod) {  //可能爆ll的计算a*b%mod
    ll r=((long double)a*b/mod+0.5);  //注：用这个不要用太大的模数
    r=a*b-r*mod;  
    return r<0?r+mod:r;  
}//具体为什么可以这样，搜索“骆可强 快速乘”
inline ll fastpow(ll a,ll k,const ll&mod) {  //快速幂
    ll ans=1;  
    while(k) {  
        if(k&1)  
            ans=mul(ans,a,mod);  
        a=mul(a,a,mod);  
        k>>=1;  
    }  
    return ans;  
}
bool mr(const ll&n,const ll&a) {  //用a来hack掉n，是合数就返回true
    ll t=n-1;  
    const ll p=t;  
    while(!(t&1))t>>=1; //除到t为奇数
    ll buf=fastpow(a,t,n);  
    if(buf==1||buf==p)return 0; //通过测试
    while((t<<=1)<p) { //和上面1到r-1是一样的 
        buf=mul(buf,buf,n);  
        if (buf==p)return 0; //通过测试
    }  
    return 1;  
}  
bool ptest(ll n) {  
    if(n<2ll) 
        return 0;  
    if(!(n&1))return n==2; //偶数的情况
    for(int i=0; i<L; ++i) {//小范围的判断整除
        if(n==plist[i])return true;
        if(!(n%plist[i]))return false;
    }
    return n<=10000ll?1:!mr(n,2)&&!mr(n,61)&&!mr(n,233);//用这三个素数就够了  
}  

时间复杂度O(\log n)

Pollard rho质因数分解

P4718 【模板】Pollard-Rho算法

小结论1（生日悖论）

随便找23个人,他们之中有任意两人生日相同的概率已经大于50\%

证明：

不相同的概率：\dfrac{365}{365}\times\dfrac{364}{365}\times\cdots\times\dfrac{365-23+1}{365}

相同的概率即是1减去上式，手工算/程序算可以算出是大于50\%的。

扩展：

随机选出[1,n]中的整数，期望写到\sqrt{n}个整数时出现重复。

算法思路

先通过一个方法找到一个因子s，之后递归s和\dfrac{n}{s}，递归到质数时更新答案即可。

现在要考虑怎么找一个因子。

假设合数n有一个因子q(1<q<n)，则我们一定可以找到两个数x,y使x\equiv y\pmod q且x\not\equiv y\pmod n

也就是说，我们可以用\gcd(|x-y|,n)来求出一个对应的因子q。

问题便转换成求一个(x,y)。

考虑一个数列f_0=1,f_i=F(f_{i-1})\mod n

其中F某种单变量映射,结果和只和参数相关，例如F(x)=x^2+c（c可以随机取一个）

这个函数还算随机（即算个伪随机数），所以根据上面的小结论，大约到\sqrt{n}项就可能有环。

每个f_i的值即可表示为一个类似于\rho(rho)的东西：（假设q(1<q<n)为n的一个因数）

第一张是模q意义下的，其中编号为i个节点表示f_i\mod q的值

第二张是模n意义下的，其中编号为i个节点表示f_i\mod n的值

可以发现f_{11}\equiv f_{4}\pmod q,f_{11}\not\equiv f_4\pmod n，也就是出现环了，我们就可以通过\gcd(|f_{11}-f_{4}|,n)求出q。

所以当且仅当f_{11}\equiv f_4\pmod n时匹配失败。

这个情况特别小，因为q的环一般远远比n的环要小，所以更换另一个c再匹配的复杂度可以忽略不计。

接下来就是要考虑判断环。

追及法求环

我们可以用两个变量x1,x2，初始时x1=1，x2=2，每次x1往后跳一个，x2往后跳两个（即v_{x2}=2v_{x1}）

也就是看f_{x1},f_{x2}就是看f_i和f_{2i}。

那么如果q有环,前者一定会被后者追上（即多跑整数个环）。

如果f_i=f_{2i}=0，则匹配失败。

否则如果\gcd(n,|f_{2i}-f_{i}|)>1则匹配成功。

根据小结论，这个是期望\sqrt{q}的。

若失败，则重置F(x)=x^2+c的c即可。

小优化

由于q\leq\sqrt{n}，得\sqrt{q}\leq n^{0.25}，所以我们要计算n^{0.25}次\gcd

也就是说，时间复杂度为O(n^{0.25}\log n)。

我们发现当\gcd(a,n)>1,\gcd(b,n)\geq1时，\gcd(ab\mod n,n)>1

我们可以一次存下z要测试的个值，同时求出它们的积\mod n，然后求与n的\gcd

大于1时，这z个值一定有一个使与n的\gcd大于1，找到它即可。

否则重新再找z个值，不断重复上面的操作。

一般z取128较合适，这样我们就优化掉了一个\log

时间复杂度O(n^{0.25})

#include <bits/stdc++.h>
using namespace std;//以下大段为米勒来宾素数测试
const int plist[]= {2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71,73,79,83,89,97};
const int L=25;
typedef long long ll;
inline ll read() {
    bool f=1;
    ll x=0;
    int c=getchar();
    while(c<=47||c>=58)f&=c!=45,c=getchar();
    while(c>=48&&c<=57) {
        x=(x<<3)+(x<<1)+(c&15);
        c=getchar();
    }
    return f?x:-x;
}
inline ll mul(const ll&a,const ll&b,const ll&mod) {
    ll r=((long double)a*b/mod+0.5);
    r=a*b-r*mod;
    return r<0?r+mod:r;
}
inline ll fastpow(ll a,ll k,const ll&mod) {
    ll ans=1;
    while(k) {
        if(k&1)
            ans=mul(ans,a,mod);
        a=mul(a,a,mod);
        k>>=1;
    }
    return ans;
}
bool mr(const ll&n,const ll&a) {
    ll t=n-1;
    const ll p=t;
    while(!(t&1))t>>=1;
    ll buf=fastpow(a,t,n);
    if(buf==1||buf==p)return 0;
    while((t<<=1)<p) {
        buf=mul(buf,buf,n);
        if (buf==p)return 0;
    }
    return 1;
}
bool ptest(ll n) {
    if(n<2ll)
        return 0;
    if(!(n&1))return n==2;
    for(int i=0; i<L; ++i) {
        if(n==plist[i])return true;
        if(!(n%plist[i]))return false;
    }
    return n<=10000ll?1:!mr(n,2)&&!mr(n,61)&&!mr(n,233);
}//上面都是米勒来宾素数测试
inline ll gcd(ll a,ll b) {//gcd
    a=abs(a),b=abs(b);
    if(!a)return b;
    for(; b; swap(a,b))
        a%=b;
    return a;
}
const int lim=128;//z
ll sav[150];//存z个值的数组
int tot;
ll prho(const ll&n,const ll&c) {//追及法求环
    ll x1=c+1,x2=mul(x1,x1,n)+c,buf=1;//先取第一个fi
    tot=0;
    if(x1>=n)x1-=n;
    if(x2>=n)x2-=n;
    while(1) {
        buf=mul(x1-x2,buf,n);//乘上要测试的数的值
        ++tot;
        sav[tot]=x1-x2;//存下要测试的数
        if(tot==lim) {//够z个数了
            if(gcd(buf,n)>1)break;//满足要求
            tot=0;//重新开始存
        }
        x1=mul(x1,x1,n)+c;
        x2=mul(x2,x2,n)+c;
        x2=mul(x2,x2,n)+c;//往后跳
    }
    for(int i=1; i<=tot; ++i) {
        buf=gcd(sav[i],n);
        if(buf>1)return buf;//返回值
    }
    return n;//不成功
}
void solve(const ll&n,ll&ans) {
    if(n<=ans)return ;//剪枝
    if(ptest(n)) {//是素数
        ans=max(ans,n);
        return ;
    }
    const ll p=n-1;
    ll sav=prho(n,rand()%p+1);
    while(sav==n)sav=prho(n,rand()%p+1);//没成功就换一个c值
    solve(sav,ans);//递归
    solve(n/sav,ans);
}
ll pollard_rho(ll n) {
    ll ans=0;
    if(!(n&1)) {
        ans=2;
        while(!(n&1))n>>=1;
    }
    const int p=sqrt(min(10000ll,n));
    for(int i=3; i<=p; i+=2)
        if(!(n%i)) {
            ans=i;
            n/=i;
            while(!(n%i))n/=i;
        }//上面是先把小范围的素数判了
    solve(n,ans);
    return ans;
}
int main() {
    srand(233);
    for(int T=read(); T; --T) {
        ll n=read();
        ll ans=pollard_rho(n);
        if(ans==n)printf("Prime\n");
        else printf("%lld\n",ans);
    }
    return 0;
}