POJ 1811 Prime Test(Pollard rho整数分解+miller_rabin素数試験)

2264 ワード

何も言わないでmillerrabinと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
"); else { if(n&1){ find(n,C); printf("%lld
",mini); } else printf("2
"); } } return 0; }