BZOJ3667: Rabin-Miller演算法

Description

Input

第一行：CAS,代表數據組數（不大於350），以下CAS行，每行一個數字，保證在64位長整形範圍內，並且沒有負數。你需要對於每個數字：第一，檢驗是否是質數，是質數就輸出Prime 
第二，如果不是質數，輸出它最大的質因數是哪個。 

Output

第一行CAS(CAS<=350，代表測試數據的組數) 
以下CAS行：每行一個數字，保證是在64位長整形範圍內的正數。 
對於每組測試數據：輸出Prime，代表它是質數，或者輸出它最大的質因數，代表它是和數 

Sample Input

Sample Output

HINT

數據範圍： 

保證cas<=350，保證所有數字均在64位長整形範圍內。 

Source

rho的裸題

註意這題卡$\log{n}$的快速乘

參考了一下黃學長的，Get到了$O(1)$快速乘的騷操作:grin:

詳細見代碼吧

#include<iostream>
#include<cstdio>
#include<cstring>
#include<cmath>
#include<cstdlib>
#define LL long long 
using namespace std;
const LL MAXN=2*1e7+10;
const LL INF=1e7+10;
inline char nc()
{
    static char buf[MAXN],*p1=buf,*p2=buf;
    return p1==p2&&(p2=(p1=buf)+fread(buf,1,MAXN,stdin),p1==p2)?EOF:*p1++;
}
inline LL read()
{
    char c=nc();LL x=0,f=1;
    while(c<'0'||c>'9'){if(c=='-')f=-1;c=nc();}
    while(c>='0'&&c<='9'){x=x*10+c-'0';c=nc();}
    return x*f;
} 
LL mx=0;
LL num[15]={2,3,5,7,11,13,17,19};
LL fastmul(LL a,LL b,LL p)
{
    LL tmp=(a*b-(LL)((long double)a/p*b+1e-8)*p);
    return tmp<0?tmp+p:tmp;
}
LL fastpow(LL a,LL p,LL mod)
{
    LL base=1;
    while(p)
    {
        if(p&1) base=(base*a)%mod;
        a=(a*a)%mod;
        p>>=1;
    }
    return base;
}
bool MR(LL n)
{
    if(n==2) return 1;
    for(LL i=0;i<8;i++) if(n==num[i]) return 1;
    if(n==1||( (n&1)==0)) return 0;
     
    LL temp=n-1,t=0;
    while( (temp&1)==0) temp>>=1,t++;
    
    for(LL o=0;o<8;o++)
    {
        LL a=num[o];
        LL now=fastpow(a,temp,n),nxt=now;
        for(LL i=1;i<=t;i++)
        {
            nxt=fastmul(now,now,n);
            if(nxt==1&&now!=1&&now!=n-1) return 0;
            now=nxt;
        }
        if(now!=1) return false;
    }
    return 1;
}
LL gcd(LL a,LL b)
{
    return b==0?a:gcd(b,a%b);
}
LL rho(LL n,LL c)
{
    LL x=rand()%n,y=x,k=2,p=1;
    for(LL i=1;p==1;i++)
    {
        x=( fastmul(x,x,n)+c )%n;
        p=gcd(abs(y-x),n);
        if(i==k) y=x,k+=k;
    }
    return p;
}
void find(LL now)
{
    if(now==1) return ;
    if(MR(now)) {mx=max(mx,now);return ;}
    LL t=now;
    while(t==now) 
        t=rho(now,rand()%(now-1)+1);//未找到因數之前無限遞歸 
    find(t);
    find(now/t);
}
int main()
{
    #ifdef WIN32
    freopen("a.in","r",stdin);
    #else
    #endif
    LL N=read();
    while(N--)
    {
        mx=0;
        LL p=read();
        find(p);
        if(mx==p) printf("Prime\n");
        else       printf("%lld\n",mx);
    }
    return 0;
}