bzoj 3667: Rabin-Miller算法

3667: Rabin-Miller算法

Time Limit: 60 Sec  Memory Limit: 512 MB

Description

Input

第一行：CAS,代表数据组数（不大于350），以下CAS行，每行一个数字，保证在64位长整形范围内，并且没有负数。你需要对于每个数字：第一，检验是否是质数，是质数就输出Prime
第二，如果不是质数，输出它最大的质因子是哪个。

Output

第一行CAS(CAS<=350，代表测试数据的组数)
以下CAS行：每行一个数字，保证是在64位长整形范围内的正数。
对于每组测试数据：输出Prime，代表它是质数，或者输出它最大的质因子，代表它是和数

Sample Input

6
2
13
134
8897
1234567654321
1000000000000

Sample Output

Prime
Prime
67
41
4649
5

HINT

数据范围：

保证cas<=350，保证所有数字均在64位长整形范围内。

Source

Rabin-Miller素数测试、Pollard-rho大因数分解 的模板题

附一个算法讲解，非常妙

http://www.cnblogs.com/Doggu/p/MillerRabin_PollardRho.html

orz 古大爷

#include<map>
#include<cmath>
#include<queue>
#include<cstdio>
#include<cstring>
#include<iostream>
#include<algorithm>
#include<cstdlib>
using namespace std;
#define inf 1000000007
#define ll long long
#define N 100010
inline ll rd()
{
	ll x=0,f=1;char ch=getchar();
	while(ch<'0'||ch>'9'){if(ch=='-')f=-1;ch=getchar();}
	while(ch>='0'&&ch<='9'){x=x*10+ch-'0';ch=getchar();}
	return x*f;
}
ll T,x,pri[10]={0,2,3,5,7,11,13};
ll mul(ll a,ll b,ll p)
{
	ll tmp=a*b-(ll)((long double)a/p*b+1e-8)*p;
	if(tmp<0) tmp+=p;return tmp;
}
ll ksm(ll a,ll b,ll p)
{
	ll sum=1;
	while(b)
	{
		if(b&1) sum=mul(sum,a,p);
		a=mul(a,a,p);b>>=1;
	}
	return sum;
}
bool MR(ll n)
{
	if(n==2) return 1;
	if(n<2||!(n&1)) return 0;
	for(int i=1;i<=6;i++)
	{
		if(n==pri[i]) return 1;
		if(ksm(pri[i],n-1,n)!=1) return 0;
	}
	ll t,a,tp;
	for(int i=1;i<=6;i++)
	{
		t=n-1;a=pri[i];
		while(!(t&1))
		{
			tp=ksm(a,t>>1,n);
			if(tp==n-1) break;;
			if(tp^1) return 0;
			t>>=1;
		}
	}
	return 1;
}
ll mx;
ll gcd(ll a,ll b){return b?gcd(b,a%b):a;}
ll f(ll x,ll n,ll c)
{
	ll y=mul(x,x,n)+c;
	if(y>=n) y-=n;
	return y;
}
ll rho(ll n)
{
	ll c=406,x,y,p;
	for(;;c++)
	{
		x=rand()%n;
		y=f(x,n,c);
		while(1)
		{
			p=gcd(abs(x-y),n);
			if(p^1) return p;
			x=f(x,n,c);y=f(f(y,n,c),n,c);
			if(x==y) break;
		}
	}
}
void sol(ll n)
{
	if(n==1||n<mx) return ;
	if(MR(n)){mx=max(mx,n);return;}
	ll tp=rho(n);
	sol(tp);sol(n/tp);
}
int main()
{
	srand(19740727);
	T=rd();
	while(T--)
	{
		x=rd();
		if(MR(x)){puts("Prime");continue;}
		mx=0;sol(x);
		printf("%lld\n",mx);
	}
	return 0;
}