日渐成熟的Pollard_Rho板子

namespace PR{
	map<ll,int >book;
	ll mul(ll a,ll b,ll MOD)
	{return (a*b-(ll)((long double)a/MOD*b)*MOD+MOD)%MOD;}
	ll pp[]={2,3,5,7,11,13,17};
	ll ksm(ll a,ll b,ll MOD)
	{
		ll ans=1;
		for(;b;b>>=1,a=mul(a,a,MOD))
			if(b&1LL)ans=mul(ans,a,MOD);
		return ans;
	}
	bool ip(ll x)
	{
		if(x==1)return 0;
		for(int i=0;i<7;i++)
			if(x==pp[i])return 1;
			else if(x%pp[i]==0)return 0;
		ll tmp=x-1;
		int cnt=0;
		while((tmp&1)^1)++cnt,tmp>>=1;
		for(int i=0;i<7;i++){
			ll val=ksm(pp[i],tmp,x),nxt;
			for(int op=1;op<=cnt;op++){
				nxt=mul(val,val,x);
				if(nxt==1&&val!=1&&val!=x-1)return 0;
				val=nxt;
			}
			if(val!=1)return 0;
		}
		for(int i=0;i<3;i++){
			ll val=ksm(rand()%(x-1)+1,tmp,x),nxt;
			for(int op=1;op<=cnt;op++){
				nxt=mul(val,val,x);
				if(nxt==1&&val!=1&&val!=x-1)return 0;
				val=nxt;
			}
			
			if(val!=1)return 0;
		}
		return 1;
	}
	ll Abs(ll x)
	{return x<0?-x:x;}
	ll Found(ll x,ll c)
	{
		ll a=rand()%(x-1)+1,b=a,res;
		for(int i=2,j=2;;i++){
			a=(mul(a,a,x)+c)%x;
			if(a==b)return -1;
			if((res=__gcd(x,Abs(a-b)))>1&&res<x)return res;
			if(i==j){
				b=a;
				j<<=1;
			}
		}
	}
	void Div(ll x)
	{
		
		if(x==1||ip(x)){
			book[x]++;
			return;
		}
		ll t,c=107;
		while((t=Found(x,c--))<0);
		Div(t),Div(x/t);
	}
}
posted @ 2019-05-12 19:46  蒟蒻小果冻  阅读(211)  评论(0编辑  收藏  举报