POJ 1811 Prime Test (miller_rabin + pollard_rho)

题目：http://poj.org/problem?id=1811

题意：判断一个数N（2 <= N < 2⁵⁴）是否为素数，若为素数输出“Prime”,否则输出最小的素因子

^{分析：先用miller_rabin随机算法判断N是否为素数，若是合数，用pollard_rho算法进行整数分解。}

^代码：

#include <iostream>
#include <cstdio>
#include <vector>
#include <algorithm>
using namespace std;

#define TIMES 5
typedef long long LL;

LL minf;

LL random(LL n)  
{  
    return (LL)((double) rand()/RAND_MAX*n+0.5);  
}

LL gcd(LL a,LL b)
{
	return b==0?a:gcd(b,a%b);
}

LL multi(LL a,LL b,LL mod) //a*b%mod
{
	LL ans=0;
	while(b)
	{
		if(b&1)
			ans=(ans+a)%mod;
		b>>=1;
		a=(a<<1)%mod; 
	}
	return ans;
}

LL kpow(LL a,LL n,LL mod) //a^n%mod  
{  
    LL ans=1;  
    while(n)  
    {  
        if(n&1)  
            ans=multi(ans,a,mod);  
        n>>=1;  
        a=multi(a,a,mod);
    }  
    return ans;  
}  
  
bool witness(LL a,LL n)  //二次探测 
{  
    LL m=n-1;  
    LL j=0;  
    while(! (m&1))  
    {  
        j++;  
        m>>=1;  
    }  
    LL x=kpow(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 false;  
    if(n==2) return true;  
    if(!(n&1)) return false;  
    for(LL 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<minf)
			minf=n;
		return ; 
	}
	LL p=n;
	while(p>=n)
		p=pollard_rho(p,k--);
	find(p,k);
	find(n/p,k);
}

int main()
{
	LL n;
	int ncase;
	scanf("%d",&ncase);
	while(ncase--)
	{
		scanf("%lld",&n);
		if(miller_rabin(n))
		{
			printf("Prime\n");
			continue;
		}
		minf=1ll<<60;
		find(n,201);
		cout<<minf<<"\n";
	}
	return 0;
}

附上kuangbin 效率更高的代码链接： http://www.cnblogs.com/kuangbin/archive/2012/08/19/2646404.html