天天看點

POJ 1811 Prime Test(Pollard rho整數分解+miller_rabin素數測試)

不說什麼了,miller_rabin和pollard_rho都看的雲裡霧裡的,直接抄的模闆,無顔見江東父老 =_=

代碼抄自HIT《數論及應用》

#include <cstdio>
#include <cstdlib>
using namespace std;
typedef long long LL;
const int C =201;
#define MAX ((LL)1<<61)
#define Times 11
LL mini;
LL gcd(LL a,LL b)
{
    return b==0?a:gcd(b,a%b);
}
LL random(LL n)
{
    return (LL)((double)rand()/RAND_MAX*n+0.5);
}
LL Multi(LL a,LL b,LL m)
{
    LL ret=0;
    while(b)
    {
        if(b&1)
            ret=(ret+a)%m;
        b>>=1;
        a=(a<<1)%m;
    }
    return ret;
}
LL quick_mod(LL a,LL b,LL m)
{
    LL ret=1;
    while(b)
    {
        if(b&1)
            ret=Multi(ret,a,m);
        a=Multi(a,a,m);
        b>>=1;
    }
    return ret;
}
bool Witness(LL a,LL n)
{
    LL m=n-1;
    int j=0;
    while(!(m&1)) j++,m>>=1;
    LL x=quick_mod(a,m,n);
    if(x==1||x==n-1) return false;
    while(j--)
    {
        x=x*x%n;
        if(x==n-1) return false;
    }
    return true;
}
bool miller_rabin(LL n)
{
    if(n==2) return true;
    if(n<2||!(n&1)) return false;
    for(int i=1;i<=Times;i++)
    {
        LL a=random(n-2)+1;
        if(Witness(a,n)) return false;
    }
    return true;
}
LL pollard_rho(LL n,int c)
{
    LL x,y,d,i=1,k=2;
    x=random(n-1)+1;
    y=x;
    while(1)
    {
        i++;
        x=(Multi(x,x,n)+c)%n;
        d=gcd(y-x,n);
        if(1<d&&d<n) return d;
        if(y==x)
            return n;
        if(i==k)
        {
            y=x;
            k<<=1;
        }
    }
}
void find(LL n,int k)///遞歸分解因子
{
    if(n==1) return ;
    if(miller_rabin(n))
    {
        if(n<mini)
            mini=n;
        return ;
    }
    LL p=n;
    while(p>=n)
        p=pollard_rho(p,k--);
    find(p,k);
    find(n/p,k);
}
int main()
{
    int total;
    scanf("%d",&total);
    while(total--)
    {
        LL n;
        scanf("%lld",&n);
        mini=MAX;
        if(miller_rabin(n))
            printf("Prime\n");
        else {
           if(n&1){
                find(n,C);
                printf("%lld\n",mini);
            }
            else printf("2\n");
        }
    }
    return 0;
}