BZOJ 3667: Rabin-Miller算法 (Pollard-Rho 模板)

说实话，我知道每一步都干啥，但我完全不知道为啥这么做，也不知道为什么是正确的，反正会用就行了~

#include <cmath>   
#include <cstdio>  
#include <algorithm>      
#define ll long long 
#define ull unsigned long long 
#define setIO(s) freopen(s".in","r",stdin)  
using namespace std;    
int array[20]={2,3,5,7,11,13,17,19}; 
ll Max;   
ll mult(ll x,ll y,ll mod) 
{
	ll tmp=(long double)x/mod*y; 
	return ((ull)x*y-(ull)tmp*mod+mod)%mod;        
}   
ll qpow(ll base,ll k,ll mod) 
{
	ll tmp=1; 
	for(;k;k>>=1,base=mult(base,base,mod)) if(k&1) tmp=mult(tmp,base,mod);   
	return tmp;   
}
ll F(ll x,ll c,ll mod) 
{
	return (mult(x,x,mod)+c)%mod;   
}
int isprime(ll x) 
{
	if(x<=1) return 1;  
	int k,i,j; 
	ll cur,a,pre; 
	for(k=0,cur=x-1;cur%2==0;cur>>=1) ++k;   
	for(i=0;i<8;++i) 
	{
		if(x==array[i]) return 1;     
		pre=a=qpow(array[i],cur,x);      
		for(j=1;j<=k;++j) 
		{ 
			a=mult(a,a,x);
			if(a==1&&pre!=1&&pre!=x-1) return 0;   
			pre=a;    
		}
		if(a!=1) return 0; 
	} 
	return 1;    
}
ll pollard_rho(ll x) 
{  
	ll s=0,t=0,c=rand()%(x-1)+1,val=1,d;    
	int k,step;  
	for(k=1;;k<<=1,s=t,val=1) 
	{ 
		for(step=1;step<=k;++step) 
		{
			t=F(t,c,x);                                
			val=mult(val,abs(s-t),x);             
			if(step%127==0) 
			{
				d=__gcd(val,x); 
				if(d>1) return d;   
			}
		}
		d=__gcd(val,x); 
		if(d>1) return d;    
	}
}
void solve(ll x) 
{ 
	if(x<=Max||x<2) return;      
	if(isprime(x)) 
	{ 
		Max=max(Max,x); 
		return;  
	} 
	ll p=x;
	for(;p>=x;) p=pollard_rho(x); 
	for(;x%p==0;) x/=p; 
	solve(x),solve(p);   
}
int main() 
{
	int T,i,j; 
	// setIO("input");  
	scanf("%d",&T); 
	for(i=1;i<=T;++i) 
	{
		ll a; 
		scanf("%lld",&a),Max=0,solve(a);   
		if(Max==a) printf("Prime\n");   
		else printf("%lld\n",Max);  
	}
	return 0; 
}