- 需要注意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;}
}
}