关于 Miller-Rabin 与 Pollard Rho 算法

· · 个人记录

P4718

Miller_Rabin

用于检测大数素性（ \sqrt{n} \ge 1e8 ）.

对于素数 P ,有费马小定理：

枚举 \lbrack 1,P ) 中的任意 a ,若均满足其逆定理 a^{P-1} \equiv 1 \: ( mod \: p ) ,则 P 大概率 为素数。

上述即费马素性检验。

存在卡迈克尔数/费马伪素数也满足上式 (如 561=3 * 11 * 17 ) ,故需要优化费马素性检验。

二次探测定理：

设当前选取底数为 a , u \cdot 2^t=n-1 \: ( u 为奇数),对于 x=u \cdot 2^y \: (y \in \lbrack 0,t \rbrack) ，满足 a^{2x} \equiv 1 \: (mod \: n) 且 a^x \not\equiv \pm1 \: (mod \: n) ，则 n 为合数，否则大概率为质数。

bool check(ll a,ll n){
    ll u=n-1,t=0;
    while(!(u&1))u>>=1,t++;
    ll x=qpow(a,u,n);
    if(x==1)return false;
    for(ll i=1,tmp;i<=t;i++,x=tmp){
        tmp=qmul(x,x,n);//原题需要转__int128类型
        if(x!=1&&x!=n-1&&tmp==1)
            return true;
    }
    return x!=1;
}
bool Miller_Rabin(ll n){
    if(n==2)return true;
    if(n<2||!(n&1))return false;
    for(ll i=1;i<=Test;i++){
        if(check(rd()%(n-1)+1,n))//mt19937
            return false;
    }
    return true;
}

Miller-Rabin的错误率是 4^{-Test} ，为提高程序效率取了 10 。（参考 这里）

Pollard Rho

用 O(n^{\frac{1}{4}}) 的时间分解大数。

考虑随机函数 f(x)=(x^2+c) \: mod \: n ，初始值为 x_0 ，令 x_1= f(x_0) ,然后再将 x_1 代入函数，数列轨迹形似 \rho 状，故称Pollard Rho.

对于 {x} ，该数列基本随机，取相邻两项取差作gcd.

基于Floyd算法的优化

数列有环存在，可以用一种Floyd的判圈方法：

设 t=r=0 ,2v_t=v_r ，则在第 i 时刻，x_t 对应 x_i , x_r 对应 x_{2i} ，两者相遇时结束。

基于Floyd判圈方法的Pollard Rho算法：

ll f(ll x,ll c,ll n){
    return (qmul(x,x,n)+c)%n;
}
ll Pollard_Rho(ll n){
    ll c=rd(n-1)+1;
    ll t=f(0,c,n),r=f(t,c,n);
    while(t!=r){
        ll d=gcd(abs(t-r),n);
        if(d>1)return d;
        t=f(t,c,N),r=f(f(r,c,n),c,n);
    }
    return N;//多测保证正确性 
}

路径倍增优化与乘法累积

上述优化瓶颈为 gcd 的 log 复杂度，由于当 (a,b)>1 时，一定有 (ac \: mod \: b ,b)>1，可以通过乘法累积优化复杂度。

考虑倍增思想，在路径中取一段 \lbrack l=2^{k-1},r=2^k \rbrack ,以所有 |x_i-x_l| \: (l < i \le r) 为样本,以避免在环上过久或gcd的复杂度过大。

需要规定一个样本数限制相乘个数，这里取 127=2^7-1（据说经实验得出迭代7次作为上界为较优方案，具体不详。）

ll f(ll x,ll c,ll n){
    return (qmul(x,x,n)+c)%n;
}
ll Pollard_Rho(ll n){
    if(!(n&1))return 2;
    ll x=rd()%(n-1)+1,y=x,c=rd()%(n-1)+1,ret;
    for(ll i=1;;i<<=1,y=x,ret=1){
        for(ll j=1;j<=i;j++){
            x=f(x,c,n);
            ret=qmul(ret,abs(x-y),n);
            if(!(j&127)){
                ll d=__gcd(ret,n);
                if(d>1)return d;
            }
        }
        ll d=__gcd(ret,n);
        if(d>1)return d;
    }
}

粗略大概就是这么多，放个模板代码。

#include<bits/stdc++.h>
#define ll long long
#define Test 10
using namespace std;
mt19937 rd(233);
ll read(){
    ll x=0;
    char ch=getchar();
    while(!isdigit(ch))ch=getchar();
    while(isdigit(ch))x=x*10+(ch^48),ch=getchar();
    return x;
}
ll qmul(ll k,ll b,ll p){
    return (__int128)k*b%p;
}
ll qpow(ll k,ll b,ll p){
    ll val=1;k%=p;
    while(b){
        if(b&1)val=qmul(val,k,p);
        k=qmul(k,k,p),b>>=1;
    }
    return val;
}
ll T,n,ans;
bool check(ll a,ll n){
    ll u=n-1,t=0;
    while(!(u&1))u>>=1,t++;
    ll x=qpow(a,u,n);
    if(x==1)return false;
    for(ll i=1,tmp;i<=t;i++,x=tmp){
        tmp=qmul(x,x,n);
        if(x!=1&&x!=n-1&&tmp==1)
            return true;
    }
    return x!=1;
}
bool Miller_Rabin(ll n){
    if(n==2)return true;
    if(n<2||!(n&1))return false;
    for(ll i=1;i<=Test;i++){
        if(check(rd()%(n-1)+1,n))
            return false;
    }
    return true;
}
ll f(ll x,ll c,ll n){
    return (qmul(x,x,n)+c)%n;
}
ll Pollard_Rho(ll n){
    if(!(n&1))return 2;
    ll x=rd()%(n-1)+1,y=x,c=rd()%(n-1)+1,ret;
    for(ll i=1;;i<<=1,y=x,ret=1){
        for(ll j=1;j<=i;j++){
            x=f(x,c,n);
            ret=qmul(ret,abs(x-y),n);
            if(!(j&127)){
                ll d=__gcd(ret,n);
                if(d>1)return d;
            }
        }
        ll d=__gcd(ret,n);
        if(d>1)return d;
    }
}
void find(ll n){
   if(n<=ans)return;
    if(Miller_Rabin(n)){
        ans=max(ans,n);
        return;
    }
   ll p=n;
   while(p>=n)p=Pollard_Rho(n);
   while(n%p==0)n/=p;
   find(p),find(n);
}
int main(){
    T=read();
    while(T--){
        n=read();
        ans=1,find(n);
        if(ans==n)printf("Prime\n");
        else printf("%lld\n",ans);
    }

    return 0;
}