51nod 1162 质因子分解

https://www.51nod.com/onlineJudge/questionCode.html#!problemId=1162

数据范围大约是2^97，需要高精度计算

可以使用pollard-rho算法，期望需要$O((logn)^{1/4})$次整数操作

需要采用 蒙哥马利约化(Montgomery reduction)避免大量取模，并每隔一段时间检查一次gcd。

#include<bits/stdc++.h>
typedef __int128 num;
typedef long long i64;
const int B=50000,K=20;
bool np[B+7];
int ps[B],pp=0,pps[B],fp=0,qp=0;
num fs[111],P,qs[111];
num read(){
num x=0;
int c=getchar();
while(!isdigit(c))c=getchar();
while(isdigit(c))x=x*10+c-48,c=getchar();
return x;
}
void pr(num x){
if(x<0)putchar('-'),x=-x;
int ss[55],sp=0;
do ss[++sp]=x%10;while(x/=10);
while(sp)putchar(48+ss[sp--]);
putchar(32);
}
num rnd(){
static std::mt19937 mt(time(0));
num z=mt();
z=z<<32^mt();
z=z<<32^mt();
return z%(P-1)+1;
}
const unsigned X=(1<<27)-1;
inline num fix1(num a){return a>=P?a-P:a;}
inline num fix2(num a){return a<0?a+P:a;}
inline num mm(num a,num b,num C=0){
unsigned b0=b&X;b>>=27;
unsigned b1=b&X;b>>=27;
unsigned b2=b&X;b>>=27;
unsigned b3=b&X;
num c=a*b3%P;
c=((c<<27)+a*b2)%P;
c=((c<<27)+a*b1)%P;
c=((c<<27)+a*b0+C)%P;
return c;
}
num iP,RR,iR;
#define HI(x) u64((x)>>52)
#define LO(x) u64(x&((1ll<<52)-1))
typedef unsigned __int128 u128;
typedef unsigned long long u64;
const num R=num(1)<<104,R1=R-1;
__attribute__((optimize("Ofast")))
inline u128 mul(u128 a,u128 b){
u128 m=a*b*iP&R1;
u128 t=(
(u128)HI(a)*HI(b)+
(u128)HI(m)*HI(P)
);
u128 t1=(
(u128)HI(a)*LO(b)+
(u128)LO(a)*HI(b)+
(u128)HI(m)*LO(P)+
(u128)LO(m)*HI(P)
);
u128 t0=(
(u128)LO(a)*LO(b)+
(u128)LO(m)*LO(P)
);
t+=HI(HI(t0)+LO(t1))+HI(t1);
return fix1(t);
}
inline u128 mulRR(u128 a){return mul(a,RR);}
num exgcd(num a,num b,num&x,num&y){
if(!b){x=1,y=0;return a;}
num g=exgcd(b,a%b,y,x);
y-=a/b*x;
return g;
}
void setP(num x){
P=x;
exgcd(P,R,iP,iR);
iP=R1&-iP;
iR=fix2(iR);
RR=mm(R%P,R%P);
}
num pw(num a,num n){
num v=1;
for(;n;n>>=1,a=mm(a,a))if(n&1)v=mm(v,a);
return v;
}
bool mr(){
num s=P-1;
int r=0;
for(;~s&1;s>>=1,++r);
num a=pw(rnd(),s);
for(;;){
if(a==1||a==P-1)return 1;
if(!--r)return 0;
a=mm(a,a);
}
}
bool isp(){
if(P<=B)return !np[P];
if(~P&1)return 0;
for(int i=0;i<15;++i)if(!mr())return 0;
return 1;
}
num gcd(num a,num b){return b?gcd(b,a%b):a;}
num abs(num x){return x>0?x:-x;}
void pollard_rho(num n){
setP(n);
if(isp()){
fs[fp++]=n;
return;
}
const int mod=4095;
for(;;){
num a,b,c,prod=mulRR(1);
a=b=mulRR(rnd());
c=mulRR(rnd());
for(i64 ii=1,k=2;;){
a=fix1(mul(a,a)+c);
prod=mul(prod,abs(a-b));
if(!prod){
prod=abs(a-b);
num g=gcd(P,mul(prod,1));
if(g==n)break;
if(g!=1){
qs[qp++]=g;
qs[qp++]=n/g;
return;
}
}
if(++ii==k)k<<=1,b=a;
if(!(ii&mod)){
num g=gcd(P,mul(prod,1));
if(g==n)break;
if(g!=1){
qs[qp++]=g;
qs[qp++]=n/g;
return;
}
}
}
}
}
void factor(num x){
for(int i=0;i<pp;++i){
for(int p=ps[i];x%p==0;fs[fp++]=p,x/=p);
}
if(x>1)qs[qp++]=x;
while(qp)pollard_rho(qs[--qp]);
std::sort(fs,fs+fp);
for(int i=0;i<fp;++i)pr(fs[i]);
}
void init(){
for(int i=2;i<=B;++i){
if(!np[i])ps[pp++]=i;
for(int j=0;j<pp&&i*ps[j]<=B;++j){
np[i*ps[j]]=1;
if(i%ps[j]==0)break;
}
}
for(int i=0;i<pp;++i){
int p=ps[i],s=1;
while(s<=B/p)s*=p;
pps[i]=s;
}
}
int main(){
init();
factor(read());
return 0;
}

View Code