模板:pollard_rho

  • 需要注意mr筛的次数稍微有点大,尽量保证正确性,sd要赋初值。
ll mul(ll x, ll y, ll p){
	ll t=(x*y-(ll)((long double)x/p*y)*p)%p;
	if (t<0) return t+p; return t;
}
ull sd;ll rd(){sd^=sd>>13,sd^=sd<<7,sd^=sd>>23;return sd>>1;}
ll gcd(ll x,ll y){return (x%y==0)?y:gcd(y,x%y);}
namespace miller_rabin{
	const int a[6]={2,3,5,7,11,61};
	ll ksm(ll x,ll y,ll p){ll s=1;for(;y;y/=2,x=mul(x,x,p)) if (y&1) s=mul(s,x,p); return s;}
	ll p,q,c;
	int chk(ll x){
		if (ksm(x,p-1,p)!=1) return 1;
		x=ksm(x,q,p);
		if (x==p-1||x==1) return 0;
		int k=0;
		while (k<=c&&x!=p-1) k++,x=mul(x,x,p);
		if (k>c) return 1;
	}
	int mr(ll p0){
		p=p0; if (p==1||!(p&1)) return p==2;
		q=p-1,c=0;
		while (q&1) q>>=1,c++;
		for(int t=0;t<6;t++) if (chk((a[t]-1)%(p-1)+1)) return 0;
		for(int t=0;t<8;t++) if (chk(rd()%(p-1)+1)) return 0;
		return 1;
	}
}
using miller_rabin::mr;

void getd(ll P){
	if (P==1) return;
	if (mr(P)) {
		for(int i=1;i<=tot;i++) if (p[i]==P) {c[i]++;return;}
		tot++,p[tot]=P,c[tot]=1;
		return;
	}
	while (1){
		ll c=rd()%(P-1)+1,t1=rd()%(P-1)+1,t2=(mul(t1,t1,P-1)+c)%(P-1)+1,g=1;
		int cnt=0;
		while (t1!=t2){
			ll t=abs(t2-t1);
			g=mul(g,t,P);
			if (!g){
				ll d=gcd(t,P);
				getd(d),getd(P/d);
				return;
			}
			if (++cnt==127){	
				cnt=0;
				ll d=gcd(g,P);
				if (d>1) {getd(d),getd(P/d);return;}
			}
			t1=(mul(t1,t1,P-1)+c)%(P-1)+1;
			t2=(mul(t2,t2,P-1)+c)%(P-1)+1;
			t2=(mul(t2,t2,P-1)+c)%(P-1)+1;
		}
		ll d=gcd(g,P);
		if (d>1) {getd(d),getd(P/d);return;}
	}
}
上一篇:matplotlib 画极坐标下圆孔应力集中现象的应力云图


下一篇:miller_rabin与pollard_rho