椭圆曲线质因数分解2

Posted 2020-12-26 tmzbot

#include <cstring>
#include <cstdio>
#include <cmath>
#include <utility>
#include <algorithm>
#include <chrono>
#include <random>
#include <vector>
#include <stdint.h>
typedef uint64_t u64;
typedef __uint128_t u128;
typedef __int128_t i128;
namespace Timer{
    template<class T,class...Args>
    std::pair<u64,T> Time(T(*func)(Args...ar),Args...ar){
        using hrc=std::chrono::high_resolution_clock;
        hrc::time_point start=hrc::now();
        T out=func(ar...);
        return std::pair<u64,T>(u64(std::chrono::duration_cast<std::chrono::nanoseconds>(hrc::now()-start)),out);
    }
}
inline u128 getint() {
    u128 ret=0;
    bool ok=0,neg=0;
    for(;;) {
        int c=getchar();
        if(c>=‘0‘&&c<=‘9‘) ret=(ret<<3)+ret+ret+c-‘0‘,ok=1;
        else if(ok)return neg?0-ret:ret;
        else if(c==‘-‘) neg=1;
    }
}
void printint(u128 n) {
    const u64 ten18=u64(1e18);
    if (n>=ten18) printf("%llu%018llu",u64(n/ten18),u64(n%ten18));
    else printf("%llu",u64(n));
}
#define rep(i,a,n) for (int i=a;i<n;i++)
struct u256 {
    u256() {}
    u256(u128 l,u128 h):lo(l),hi(h) {}
    static u256 mul128(u128 a,u128 b) {
        u64 a_hi=a>>64,a_lo=u64(a);
        u64 b_hi=b>>64,b_lo=u64(b);
        u128 p01=u128(a_lo)*b_lo;
        u128 p12=u128(a_hi)*b_lo+u64(p01>>64);
        u64 t_hi=p12>>64,t_lo=p12;
        p12=u128(a_lo)*b_hi+t_lo;
        u128 p23=u128(a_hi)*b_hi+u64(p12>>64)+t_hi;
        return u256(u64(p01)|(p12<<64),p23);
    }
    u128 lo,hi;
};
struct Mont{
    Mont(u128 n):mod(n) {
        inv=n;
        rep(i,0,6) inv*=2-n*inv;
        r2=-n%n;
        rep(i,0,4) if ((r2<<=1)>=mod) r2-=mod;
        rep(i,0,5) r2=mul(r2,r2);
    }
    u128 reduce(u256 x) const {
        u128 y=x.hi-u256::mul128(x.lo*inv,mod).hi;
        return i128(y)<0?y+mod:y;
    }
    u128 reduce(u128 x) const { return reduce(u256(x,0)); }
    u128 init(u128 n) const { return reduce(u256::mul128(n,r2)); }
    u128 mul(u128 a,u128 b) const { return reduce(u256::mul128(a,b)); }
    u128 mod,inv,r2;
};
// the Min-25 montgomery form manipulator
u128 ctz(u128 x){int a=__builtin_ctzll(u64(x>>64))+64,b=__builtin_ctzll(u64(x));return u64(x)?b:a;}
u128 gcd(u128 a,u128 b) {
    if (b==0) return a;
    int shift=ctz(a|b);
    b>>=ctz(b);
    while (a) {
        a>>=ctz(a);
        if (a<b) std::swap(a, b);
        a-=b;
    }
    return b<<shift;
}
i128 invert(i128 a,i128 b){
    if(!a||!b)return 0;
//  putchar(‘:‘),printint(u128(a)),putchar(‘ ‘),printint(u128(b)),putchar(‘
‘);
    bool truth=0;
    if(a<0)a=-a,truth=1;
    i128 b_or=b,alpha=1,beta=0;
    while(!(a&1)){
        if(alpha&1)alpha+=b_or;
        alpha>>=1,a>>=1;
    }if(b>a)std::swap(a,b),std::swap(alpha,beta);
    while(b&&(a^b)){
        a-=b;alpha-=beta;
        while(!(a&1)){
            if(alpha&1)alpha+=b_or;
            alpha>>=1,a>>=1;
        }if(b>a)std::swap(a,b),std::swap(alpha,beta);
    }
    if(a==b)b=0,alpha-=beta,std::swap(alpha,beta);
//  putchar(‘:‘),printint(u128(alpha)),putchar(‘ ‘),printint(u128(beta)),putchar(‘
‘);
//  putchar(‘:‘),printint(u128(-alpha)),putchar(‘ ‘),printint(u128(-beta)),putchar(‘
‘);
    if(truth)alpha=b_or-alpha;
    alpha=alpha%b_or;
    if(alpha<0)alpha+=b_or;
    if(a!=1)return 0;
    return alpha;
}
//invert and gcd
u64 sqrt_approx(u64 x){
    u64 approx=sqrt(double(x));
    return (approx+x/approx)>>1;
}
u64 sqrt(u64 x){
    u64 approx=sqrt(double(x));
    u64 apt=(approx+x/approx)>>1;
    approx=apt*apt;
    if(approx>x)return apt-1;
    if(x-approx>=2*apt-1)return apt+1;
    return apt;
}
u128 sqrt(u128 r){
    if(!(r>>64))return sqrt(u64(r));
    int cnt=(((64-__builtin_clzll(u64(r>>64)))+1)|1)^1;
    u128 approx=u128(sqrt_approx(u64(r>>cnt)))<<(cnt/2);
    approx=(approx+r/approx)>>1;
    u128 apt=u128(u64(approx))*u128(u64(approx));
    return approx-((r-apt)>>127);
}
// fast int128 square root
#define ModularManipulate     u128 n=Modu->mod;     const auto add=[&] (u128 x,u128 y) { return (x+=y)>=n?x-n:x; };     const auto sub=[&] (u128 x,u128 y) { return i128(x-=y)<0?x+n:x; };     const auto mul=[&] (u128 x,u128 y) { return Modu->mul(x,y); };     const auto get=[&] (u128 x)        { return Modu->reduce(x); };     const auto set=[&] (u128 x)        { return Modu->init(x); };     const auto dbl=[&] (u128 x)        { return (x<<=1)>=n?x-n:x; };

u128 invert(u128*inv,u128*lis,int len,Mont*Modu){
    ModularManipulate
    for(int i=1;i<len;++i)
        inv[i-1]=lis[i],
        lis[i]=mul(lis[i],lis[i-1]);
    u128 invt=u128(invert(get(lis[len-1]),n));
//  printint(get(lis[len-1])),putchar(‘ ‘),printint(invt),putchar(‘
‘);
    if(!invt){
        while(~--len){
            u128 factor=gcd(lis[len],n);
            if(factor==1)break;
            if(factor<n)return factor;
        }return 1;
    }invt=set(invt);
    for(int i=len-1;i;--i)
        inv[i]=mul(invt,lis[i-1]),
        invt=mul(invt,inv[i-1]);
    inv[0]=invt;
    return 0;
}
const int maxn=10010;
// invert a list of u128 in parallel, while returning 1 indicates failure, returning 0 indicates inverted,
// returning other indicates successful factorization
struct affine{u128 x,y,c;};

affine tempaff[10][maxn];
u128 tempui[10][maxn];//

u128 Add(affine*p1,affine*p2,int len,Mont*Modu){
    ModularManipulate
    u128*inv=tempui[0],*invr=tempui[1];
    for(int i=0;i<len;++i)
        inv[i]=sub(p1[i].x,p2[i].x);
    u128 k=invert(invr,inv,len,Modu);
    if(k)return k;
    for(int i=0;i<len;++i){
        k=mul(sub(p1[i].y,p2[i].y),invr[i]);
        p2[i].x=sub(sub(mul(k,k),p1[i].x),p2[i].x);
        p2[i].y=sub(mul(k,sub(p1[i].x,p2[i].x)),p1[i].y);
    }return 0;
}
u128 Addsub_x(affine*p1,affine*p2,u128*sum,u128*dif,int len,Mont*Modu){
    ModularManipulate
    u128*inv=tempui[0];
    for(int i=0;i<len;++i)
        sum[i]=sub(p2[i].x,p1[i].x);
    u128 k=invert(inv,sum,len,Modu),r;
    if(k)return k;
    for(int i=0;i<len;++i){
        k=mul(sub(p2[i].y,p1[i].y),inv[i]);
        r=mul(add(p2[i].y,p1[i].y),inv[i]);
        sum[i]=sub(sub(mul(k,k),p1[i].x),p2[i].x);
        dif[i]=sub(sub(mul(r,r),p1[i].x),p2[i].x);
    }return 0;
}
u128 pow(u128 base,u128 exp,Mont*Modu){
    ModularManipulate
    u128 ca[4];
    ca[0]=1;ca[1]=base;
    ca[2]=mul(base,base),ca[3]=mul(ca[2],base);
    bool f=0;
    for(int i=126;i>=0;i-=2){
        if(f)ca[0]=mul(ca[0],ca[0]),ca[0]=mul(ca[0],ca[0]);
        int q=(exp>>i)&3;
        if(q)f=1,ca[0]=mul(ca[0],ca[q]);
    }return ca[0];
}
u128 Double(affine*p1,int len,Mont*Modu){
    ModularManipulate
    u128*inv=tempui[0],*invr=tempui[1];
    for(int i=0;i<len;++i)
        inv[i]=dbl(p1[i].y);
    u128 k=invert(invr,inv,len,Modu);
    if(k)return k;
    for(int i=0;i<len;++i){
        u128 r=p1[i].x;
        k=mul(r,r);
        k=mul(add(dbl(k),add(k,p1[i].c)),invr[i]);
        p1[i].x=sub(mul(k,k),dbl(r));
        p1[i].y=sub(mul(k,sub(r,p1[i].x)),p1[i].y);
    }return 0;
}
u128 Sub(affine*p1,affine*p2,int len,Mont*Modu){
    ModularManipulate
    u128*inv=tempui[0],*invr=tempui[1];
    for(int i=0;i<len;++i)
        inv[i]=sub(p2[i].x,p1[i].x);
    u128 k=invert(invr,inv,len,Modu);
    if(k)return k;
    for(int i=0;i<len;++i){
        k=mul(add(p1[i].y,p2[i].y),invr[i]);
        p2[i].x=sub(sub(mul(k,k),p1[i].x),p2[i].x);
        p2[i].y=add(mul(k,sub(p1[i].x,p2[i].x)),p1[i].y);
    }return 0;
}
u128 NAFConv(u64 E){//NAF with a leading 01. So use with E<2^63
    u128 res=1;
    while(E){
        if(E&1)res=(res<<2)|(E&3),E-=2-(E&3);
        else res<<=2;
        E>>=1;
    }return res;
}
#define prr(x) printpoints(x,len,Modu)
void printpoints(affine*af,int len,Mont*Modu){
    ModularManipulate
    printf("Count:
%d
[",len);
    for(int i=0;i<len;++i){
        putchar(‘[‘);
        printint(get(af[i].x)),putchar(‘,‘),
        printint(get(af[i].y)),putchar(‘,‘),
        printint(get(af[i].c)),puts("],");
    }
    puts("]");
}
u128 FastMultiply(affine*p1,u64 d,int len,Mont*Modu){
    if(d==1)return 0;
    u128 Na=NAFConv(d);
    affine*tem=tempaff[0];
    std::copy(p1,p1+len,tem);
    Na>>=2;
//  prr(p1);
    while(Na!=1){
        int op=Na&3;
        u128 k=Double(p1,len,Modu);
//      puts("*2");
//      prr(p1);
        if(k)return k;
        if(op==1)k=Add(tem,p1,len,Modu);//,puts("+1"),prr(p1),prr(tem);
        else if(op==3)k=Sub(tem,p1,len,Modu);//,puts("-1"),prr(p1),prr(tem);
        if(k)return k;
        Na>>=2;
    }return 0;
}
u128 InitPoints(u128*param,affine*points,int len,Mont*Modu){
    ModularManipulate
    u128 five=set(5),two=set(2),one=set(1);
    for(int cn=0;cn<len;++cn){
        u128 sigma=param[cn];
        u128 u=sub(mul(sigma,sigma),five);
        u128 v=dbl(dbl(sigma));
        u128 i=mul(mul(u,u),dbl(dbl(mul(u,v))));
        points[cn].x=u;
        points[cn].y=v;
        points[cn].c=i;
        param[cn]=mul(i,v);
    }
    u128*inv=tempui[0];
    u128 ret=invert(inv,param,len,Modu);
    if(ret)return ret;
    for(int j=0;j<len;++j){
        u128 u=points[j].x,v=points[j].y,i=points[j].c;
        u128 in=inv[j];
        u128 t1=sub(v,u),t2=add(dbl(u),add(u,v));
        t1=mul(mul(t1,t1),t1);
        u128 a=sub(mul(mul(t1,t2),in),two);
        t1=mul(u,mul(i,in));
        u128 x0=mul(t1,mul(t1,t1));
        u128 b=mul(add(x0,a),x0);
        b=mul(add(b,one),x0);
        x0=mul(b,x0);
        u128 y0=mul(b,b);
        t1=get(a);
        while(t1%3)t1+=n;
        t1=set(t1/3);
        x0=add(x0,mul(t1,b));
        t2=mul(y0,sub(one,mul(t1,a)));
        points[j].x=x0;
        points[j].y=y0;
        points[j].c=t2;
    }return 0;
}
namespace Sieve{
    typedef unsigned int u32;
    typedef unsigned long long ull;
    const char pr60[]={2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59};
    const char masks[][4]={
        {3,7,11,13},
        {3,17,19,23},
        {2,29,31},
        {2,37,41},
        {2,43,47},
        {2,53,59}
    };
    const u32 segsize=65536;
    void Apply_mask(u32*a,u32*b,u32 l1,u32 l2){
        u32 t=0;
        for(u32 q=0,r=l1/l2;q<r;++q)
            for(u32 i=0;i<l2;++i)
                a[t++]|=b[i];
        for(u32 i=0;t<l1;++i)
            a[t++]|=b[i];
    }
    void Gen_mask_sub(u32*a,u32 l1,u32 b){
        u32 st=b>>1,rt=0;
        while(rt<l1){
            a[rt]|=1<<st;
            st+=b;
            if(st>=30)st-=30,++rt;
            if(st>=30)st-=30,++rt;
        }
    }
    void PrintMask(u32*a,u32 len){
        printf("Mask of len %u
",len*60);
        for(u32 i=0;i<len;++i){
            for(u32 j=0;j<30;++j)
                if((a[i]&(1<<j)))
                    printf("%llu
",i*60ull+j*2ull+1ull);
        }
    }
    u32 Gen_mask(u32*a,int id){
        int len=masks[id][0];
        u32 ll=1;
        for(int i=1;i<=len;++i)
            ll*=masks[id][i];
        memset(a,0,4*ll);
        for(int i=1;i<=len;++i)
            Gen_mask_sub(a,ll,masks[id][i]);
    //  PrintMask(a,ll);
        return ll;
    }
    const u32 mask=0x1a4b3496;
    const u32 pr60_m=0xdb4b3491;
    u32 pr[10000][4],prl;
    std::vector<u32> main(ull ma){
        ull tma,tmx;
        tma=(ma-1)/60+1;
        tmx=tma*60;//upper limit
        u32*sieve=new u32[tma];// getting a sieve ready
        u32*maske=new u32[7429];
        std::fill(sieve,sieve+tma,mask);
        for(int i=0;i<6;++i)
            Apply_mask(sieve,maske,tma,Gen_mask(maske,i));

        ull preseg=std::min(tmx,ull(sqrt(double(ma))/60)+1);
        u32 j=61;
        for(;ull(j)*j<=preseg*60;j+=2){
            u32 v=j/60,u=(j%60)>>1;
            if(!(sieve[v]&(1<<u))){
                v=j/30,u=j%30;
                u32 rt=j*3/60,st=(j*3%60)>>1;
                while(rt<preseg){
                    sieve[rt]|=1<<st;
                    rt+=v;
                    st+=u;
                    if(st>=30)st-=30,++rt;
                }
                pr[prl][0]=v;
                pr[prl][1]=u;
                pr[prl][2]=rt;
                pr[prl][3]=st;
                prl++;
            }
        } // Non-segmented sieve core
        if(preseg==tmx)goto end;
        for(u32 segl=preseg;segl<tma;segl+=segsize){
            u32 segr=std::min(segl+segsize,u32(tma));
            for(;ull(j)*j<=segr*60;j+=2){
                u32 v=j/60,u=(j%60)>>1;
                if(!(sieve[v]&(1<<u))){
                    v=j/30,u=j%30;
                    ull t=j*ull(j);
                    u32 rt=t/60,st=t%60>>1;
                    pr[prl][0]=v;
                    pr[prl][1]=u;
                    pr[prl][2]=rt;
                    pr[prl][3]=st;
                    prl++;
                }
            }
            for(int i=0;i<prl;++i){
                u32 v=pr[i][0],u=pr[i][1],rt=pr[i][2],st=pr[i][3];
                while(rt<segr){
                    sieve[rt]|=1<<st;
                    rt+=v;
                    st+=u;
                    if(st>=30)st-=30,++rt;
                }
                pr[i][0]=v;
                pr[i][1]=u;
                pr[i][2]=rt;
                pr[i][3]=st;
            }
        }
        end:
        sieve[0]=pr60_m;
        std::vector<u32> ret;
        ret.push_back(2);
        for(u32 i=0;i<tma;++i){
            for(u32 j=0;j<30;++j)
                if(!(sieve[i]&(1<<j)))ret.push_back(i*60+j*2+1);
        }
        return ret;
    }
}
u128 factor(u128 N,int SmoothBound,int curve_count,std::vector<unsigned int> primes){
    u128*param=tempui[2];
    Mont Mod(N);
    Mont*Modu=&Mod;
    ModularManipulate
    for(int i=0;i<curve_count;++i)
        param[i]=set(693595790+i);
    affine*plist=tempaff[2];
    u128 ret=InitPoints(param,plist,curve_count,&Mod);
//  printpoints(plist,curve_count,&Mod);
    if(ret)return ret;
    u128 std=u128(1)<<63;
    u64 qbound=u64(sqrt(N));
    u64 rbound=qbound+sqrt(qbound)+1;
    for(int i:primes){
        if(i>SmoothBound)break;
        u64 t=i,g=SmoothBound/i;
        while(t<=g)t*=i;
        ret=FastMultiply(plist,t,curve_count,Modu);
//      printf("Mult %llu
",t);
//      printpoints(plist,curve_count,&Mod);
        if(ret)return ret;
    }
    return 1;
}
void Test(u128 N,int SmoothBound,int curve_count,std::vector<unsigned int> primes){
    u128*param=tempui[2];
    affine*plist=tempaff[2];
    affine*plist2=tempaff[3];
    Mont Mod(N);
    Mont*Modu=&Mod;
    ModularManipulate
    for(int i=0;i<curve_count;++i)
        param[i]=set(6+i);
    u128 ret=InitPoints(param,plist,curve_count,&Mod);
    printpoints(plist,curve_count,Modu);
    std::copy(plist,plist+curve_count,plist2);
    Double(plist,curve_count,Modu);
    printpoints(plist,curve_count,Modu);
}
int main(){
    u128 inp=getint();
    std::vector<unsigned int> pr=Sieve::main(1000000);
    u128 p=factor(inp,50000,160,pr);
    u128 q=inp/p;
    if(q<p)std::swap(p,q);
    printint(p),putchar(‘ ‘),printint(q);
//  Test(inp,200,1,pr);
    return 0;
}