BZOJ 3667 Rabin-Miller算法

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

~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

Miller-Rabin算法~

这个算法就是把一堆近似算法放在一起，各判断一次，把判错的可能性降到几乎为0。

具体做法详见http://blog.csdn.net/thy_asdf/article/details/51347390，我觉得就是代替高精度的那个取模乘法比较有用……抄起来真是累啊~

#include<cstdio>
#include<cstdlib>
#include<iostream>
#include<algorithm>
using namespace std;
#define ll long long
const ll a[]={2,3,5,7,11,13,17,19,23,29};

ll t,n,maxx;

ll read()
{
	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<<1)+(x<<3)+ch-'0';ch=getchar();}
	return x*f;
}

ll cheng(ll u,ll v,ll mod)
{
	ll res=(long double)u/mod*v+1e-8;
	res=u*v-mod*res;
	return res<0 ? res+mod:res; 
}

ll mi(ll u,ll v,ll mod)
{
	ll res=1;u%=mod;
	for(;v;v>>=1,u=cheng(u,u,mod)) if(v&1) res=cheng(res,u,mod);
	return res;
}

ll gcd(ll u,ll v)
{
	return v ? gcd(v,u%v):u;
}

bool chec(ll u,ll mod,ll k,ll v)
{
	ll now=mi(u,k,mod),las=now;
	for(int i=1;i<=v;i++)
	{
		now=cheng(now,now,mod);
		if(now==1 && las!=1 && las!=mod-1) return 0;
		las=now;
	}
	return now==1 ? 1:0;
}

bool mr(ll u)
{
	if(u<=1) return 0;
	ll res=u-1,now=0;
	while(!(res&1)) res>>=1,now++;
	for(int i=0;i<9;i++)
	{
		if(u==a[i]) return 1;
		if(!chec(a[i],u,res,now)) return 0;
	}
	return 1;
}

ll chan(ll u,ll v)
{
	ll k=2,x=rand()%u,y=x,p=1;
	for(ll i=1;p==1;i++)
	{
		x=(cheng(x,x,u)+v)%u;
		p=y>x ? y-x:x-y;p=gcd(p,u);
		if(i==k) y=x,k+=k;
	}
	return p;
}

void solve(ll u)
{
	if(u==1) return;
	if(mr(u))
	{
		maxx=max(maxx,u);return;
	}
	ll now=u;
	while(now==u) now=chan(u,rand()%(u-1));
	solve(now);solve(u/now);
}

int main()
{
	srand(20000619);t=read();
	while(t--)
	{
		n=read();maxx=0;solve(n);
		if(n==maxx) puts("Prime");
		else printf("%lld\n",maxx);
	}
	return 0;
}