BZOJ.3667.Rabin-Miller算法(MillerRabin PollardRho)

题目链接
Pollard_Rho:http://blog.csdn.net/thy_asdf/article/details/51347390

#include<cstdio>
#include<cctype>
#include<algorithm>
#define gc() getchar()
const int p[]={2,3,5,7,11,13,17,19};
typedef long long LL;
LL Ans;

inline LL read()
{
    LL now=0,f=1;register char c=gc();
    for(;!isdigit(c);c=gc()) if(c=='-') f=-1;
    for(;isdigit(c);now=now*10+c-'0',c=gc());
    return now*f;
}
inline LL Mult(LL a,LL b,LL p)//O(1)快速乘
{
    LL tmp=a*b-(LL)((long double)a/p*b+1e-8)*p;
    return tmp<0?tmp+p:tmp;
}
LL Fast_Pow(LL n,LL k,LL p)
{
    LL t=1;
    for(;k;k>>=1,n=n*n%p)
        if(k&1) t=t*n%p;
    return t;
}
bool Miller_Rabin(LL n)
{
    if(n==2) return 1;
    if(!(n&1)||n==1) return 0;
    for(int i=0;i<8;++i)
        if(n==p[i]) return 1;
        else if(!(n%p[i])) return 0;
    LL u=n-1,now,las; int t=0;
    while(!(u&1)) u>>=1,++t;
    for(int i=0;i<8;++i)
    {
        now=Fast_Pow(p[i],u,n);
        for(int j=1;j<=t;++j)
        {
            las=now, now=Mult(now,now,n);
            if(now==1&&las!=1&&las!=n-1) return 0;
        }
        if(now!=1) return 0;
    }
    return 1;
}
LL gcd(LL x,LL y)
{
    return y?gcd(y,x%y):x;
}
LL Rho(LL n,LL delta)
{//现要分解n，有两个随机数x,y，若p=gcd(x-y,n)!=1&&p!=n，那么p为n的一个约数...省略
    LL x=rand()%n,y=x,p=1; int k=2;//设定k为此次路径长
    for(int i=1;p==1;++i)
    {
        x=(Mult(x,x,n)+delta)%n;//随机函数f(x)=x*x+d
        p=gcd(std::abs(x-y),n);//多次生成随机数，直至找到p是n的一个因子
        if(i==k) y=x,k<<=1;//达到k次后把y赋值为x。路径每次倍长
    }
    return p;
}
void Find(LL n)
{
    if(n==1) return;
    if(Miller_Rabin(n)) {Ans=std::max(Ans,n);/*fac[++cnt]=n;*/ return;}
    LL t=n;
    while(t==n) t=Rho(n,rand()%(n-1)+1);
    //t=n说明这个随机函数会导致走到n的环上，再换一个重试即可
    Find(t), Find(n/t);
}

int main()
{
#ifndef ONLINE_JUDGE
    freopen("3667.in","r",stdin);
#endif

    int t=read();LL n;
    while(t--)
        n=read(),Ans=0,Find(n),Ans==n?puts("Prime"):printf("%lld\n",Ans);
    return 0;
}

原文地址：https://www.cnblogs.com/SovietPower/p/8448321.html