BZOJ3667: Rabin-Miller算法

Time Limit: 60 Sec  Memory Limit: 512 MB
Submit: 1650  Solved: 570
[Submit][Status][Discuss]

Description

 

Input

第一行：CAS,代表数据组数（不大于350），以下CAS行，每行一个数字，保证在64位长整形范围内，并且没有负数。你需要对于每个数字：第一，检验是否是质数，是质数就输出Prime 
第二，如果不是质数，输出它最大的质因子是哪个。 

Output

第一行CAS(CAS<=350，代表测试数据的组数) 
以下CAS行：每行一个数字，保证是在64位长整形范围内的正数。 
对于每组测试数据：输出Prime，代表它是质数，或者输出它最大的质因子，代表它是和数 

Sample Input

6
2
13
134
8897
1234567654321
1000000000000

Sample Output

Prime
Prime
67
41
4649
5

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
");
        else       printf("%lld
",mx);
    }
    return 0;
}