bzoj 3667 Rabin-Miller算法

#include<iostream>
#include<cstdio>
#include<algorithm>
#include<cstring>
#include<cmath>
#define ll long long
using namespace std;
int n;
ll x,mx;
ll gcd(ll a,ll b)
{
    if(!b)return a;
    return gcd(b,a%b);
}
i64 mul(i64 a,i64 b,i64 c){
    if(c<=2000000000ll)return a*b%c;
    i64 r=a*b-i64(ld(a)/c*b)*c;
    if(r>=c||r<=-c)r%=c;
    return r>=0?r:r+c;
}
ll pow(ll a,ll b,ll p)
{
    ll ans=1;a%=p;
    while(b)
    {
        if(b&1)ans=mul(a,ans,p);
        b>>=1;
        a=mul(a,a,p);
    }
    return ans;
}
bool check(ll rd,ll n,ll r,ll s)
{
    ll ans=pow(rd,r,n),p=ans;
    for(int i=1;i<=s;i++)
    {
        ans=mul(ans,ans,n);
        if(ans==1&&p!=1&&p!=n-1)return 1;
        p=ans;
    }
    if(ans!=1)return 1;
    return 0;
}
bool MR(ll n)
{
    if(n<=1)return 0;
    if(n==2)return 1;
    if(n%2==0)return 0;
    ll r=n-1,s=0;
    while(r%2==0)r/=2,s++;
    for(int i=0;i<10;i++)
    {
        if(check(rand()%(n-1)+1,n,r,s))return 0;
    }
    return 1;
}
ll rho(ll n,ll c)
{
    ll k=2,x=rand()%n,y=x,p=1;
    for(ll i=1;p==1;i++)
    {
        x=(mul(x,x,n)+c)%n;
        if(y>x)p=y-x;
        else p=x-y;
        p=gcd(n,p);
        if(i==k)y=x,k+=k;
    }
    return p;
}
void solve(ll n)
{
    if(n==1)return ;
    if(MR(n)){mx=max(n,mx);return;}
    ll t=n;
    while(t==n)t=rho(n,rand()%(n-1)+1);
    solve(t);
    solve(n/t);
}
int main()
{
    scanf("%d",&n);
    while(n--)
    {
        scanf("%lld",&x);
        mx=0;
        solve(x);
        if(mx==x)puts("Prime");
        else printf("%lld
",mx);
    }
    return 0;
}