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;
}