椭圆曲线质因数分解2

Upd: 写了一点资料. https://pan.baidu.com/s/1GyzysqQUVuUyI8Bsqs9O-g

#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('\n');
	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('\n');
//	putchar(':'),printint(u128(-alpha)),putchar(' '),printint(u128(-beta)),putchar('\n');
	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('\n');
	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:\n%d\n[",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\n",len*60);
		for(u32 i=0;i<len;++i){
			for(u32 j=0;j<30;++j)
				if((a[i]&(1<<j)))
					printf("%llu\n",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\n",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;
}