Miller_Rabin & Pollard_Rho

hubingshan · 2023-09-18 15:57:02 · 个人记录

int a[]={0,2,3,5,7,11,13,17,19,23,29,31,37};
int mul(int x,int y,int mod){
    return (__int128)x*y%mod;
}
int ksm(int x,int y,int mod){
    int ans=1;
    while(y){
        if(y&1) ans=mul(ans,x,mod);
        x=mul(x,x,mod),y>>=1;
    }
    return ans;
}
int gcd(int x,int y){
    if(y==0) return x;
    return gcd(y,x%y);
}
int MR(int x){
    if(x<3||x%2==0) return x==2;
    int y=x-1,cnt=0,num;
    while(y%2==0) cnt++,y>>=1;
    for(int i=1;i<=12;i++){
        num=ksm(a[i],y,x);
//      printf("%lld %lld %lld\n",num,y,x);
        if(num==1||num==x-1||num==0) continue;
        for(int j=1;j<=cnt;j++){
            num=mul(num,num,x);
            if(num==x-1&&j!=cnt){
                num=1;break;
            }
            if(num==1) return 0;
        }
//      printf("!!%lld %lld\n",i,num);      
        if(num!=1) return 0;
    }
    return 1;
}
int f(int x,int c,int mod){
    return (mul(x,x,mod)+c)%mod;
}
int PR(int x){
    int s=0,t=0,c=rand()%(x-1)+1;
    for(int len=1,num=1;;len*=2,num=1,s=t){
        for(int i=1;i<=len;i++){
            t=f(t,c,x);
            num=mul(num,abs(t-s),x);
            if(i%127==0){
                int p=gcd(num,x);
                if(p>1) return p;
            }
        }
        int p=gcd(num,x);
        if(p>1) return p;
    }
}
int fenjie(int n){
    if(n==1) return 1;
    if(n<=ans) return ans;
    if(MR(n)) return ans=max(ans,n);
//  printf(":: %lld\n",n);
    int x=n;
    while(x==n) x=PR(n);
    while(n%x==0) n/=x;
    return max(fenjie(x),fenjie(n));
}