pku1811 Prime Test（Miller_Rabin + Pollard_Rho）

2 <= N < 2^54,巨恶心...
一些数据，只有前4个是Prime
18
2
3
5
18014398509481931
4
6
8
10
1000009
561
1105
1729
294409
56052361
118901521
172947529
216821881
228842209

 

//390-650MS:

//rand() % X + Y = [Y,X+Y)

#include <stdio.h>
#include <time.h>
#include <stdlib.h>

typedef long long LL;

#define MAXFACN 100
#define TIMES 10

LL fac[MAXFACN];
int ftop;

LL gcd(LL a,LL b)
{
    LL t;
    if(a<b)
    {
        t=a;
        a=b;
        b=t;
    }
    while(b)
    {
        t=a%b;
        a=b;
        b=t;
    }
    return a;
}

/*LL muti_mod( LL a, LL b, LL m) //recursion is too slow,but can also solve it
{
    if( !b ) return 0;
    return ( muti_mod( a, b >> 1 , m ) * 2 + a * ( b & 1 ) ) % m;
}*/
LL muti_mod(LL a,LL b,LL m)//s=(a*b)(mod m) (n<2^62)
{
  //  return ( (a%m) * (b%m) ) % m; // this will overflow!
    a %= m;
    b %= m;
    LL s = 0;
    while(b)
    {
        if(b & 1)
        {
            s += a;
            if(s >= m) s -= m;
        }
        a <<= 1;
        if(a >= m) a -= m;

        b >>= 1;
    }
    return s;
}

inline LL mod_exp(LL a,LL b,LL m)// s=a^b(mod m)
{
    bool bit[64];
    int i=0;
    while(b)
    {
        bit[i++] = b&1;
        b>>=1;
    }
    LL s = 1;
    for(i--; i>=0; i--)
    {
        s = (muti_mod(s,s,m));
        if(bit[i]) s = muti_mod(s,a,m);
    }
    return s;
    /*   a=a%m;
       LL s=1;
       while(b>0)
       {
           if(b&1) s=muti_mod(s,a,m);
           a = muti_mod(a,a,m);
           b>>=1;
       }
       return s;*/
}

bool Witness(LL a,LL n)
{
    LL u=n-1;
    int t=0;
    while(!(u&1))
    {
        t++;
        u>>=1;
    }
    LL y,x=mod_exp(a,u,n);
    while(t--)
    {
        y = muti_mod(x,x,n);
        if(y==1 && x!=1 && x!=n-1) return true;
        x = y;
    }
    if(y!=1) return true; // composite
    return false; //may be prime
}

bool Miller_Rabin(LL n,int T=TIMES)
{
    if(n<2) return false; //composite
    if(n==2) return true; // prime
    if( !(n&1) ) return false;
    while(T--)
    {
        LL a = rand() % (n-1) +1;
        if(Witness(a,n)) return false;
    }
    return true; //may be prime
}

LL Pollard_Rho(LL n,int c)
{
    LL x,y,d,i=1,k=2;
    y = x = rand() % (n-1) + 1;
    while(1)
    {
        i++;
        x = (muti_mod(x,x,n) + c) % n;
        if(y>x) d = gcd(y-x,n);
        else d = gcd(x-y,n);
        if(1<d && d<n)
        {
     //       printf("%I64d = %I64d * %I64d\n",n,d,n/d);
            return d;
        }
        if(y==x) return n;
        if(i==k)
        {
            y=x;
            k<<=1;
        }
    }
}

void findfac(LL n)
{
//   if(n==1) return;
    if(Miller_Rabin(n))
    {
        fac[++ftop]=n;
        return;
    }
    LL p=n;
    while(p>=n) p=Pollard_Rho(p, rand() % (n-1) + 1 );
    findfac(p);
    findfac(n/p);
}

int main()
{
#ifndef ONLINE_JUDGE
    freopen("tdata.txt","r",stdin);
#endif
    srand(unsigned(time(0)));
    int n;
    LL p;
    while(scanf("%d",&n)!=EOF)
    {
        while(n--)
        {
            scanf("%I64d",&p);
            if(Miller_Rabin(p)) printf("Prime\n");
            else
            {
                ftop=-1;
                findfac(p);
                LL ans = p;
                for(int i=0; i<=ftop; i++)
                    if(ans > fac[i]) ans=fac[i];
                printf("%I64d\n",ans);
            }
        }
    }
    return 0;
}