萌新求助ECM
  • 板块学术版
  • 楼主wkywkywky
  • 当前回复1
  • 已保存回复1
  • 发布时间2025/1/30 21:21
  • 上次更新2025/1/31 15:19:36
查看原帖
萌新求助ECM
133954
wkywkywky楼主2025/1/30 21:21

求大神教我椭圆曲线分解,我自己写的代码总是跑不出来了/kel

#include<bits/stdc++.h>
#define ll __int128
#define ld long double
using namespace std;
const int B=7000;
struct node{ll v;int k;}p[110];
struct point{ll x,y;};
int cnt,tot,prime[B],test[12]={2,3,5,7,11,13,17,19,23,29,31,37};
ll N;
bitset<B+10>mark;
void init(){
  int i,j,v;
  mark.set();
  for(i=2;i<=B;i++){
    if(mark[i])prime[++tot]=i;
    for(j=1;j<=tot&&prime[j]*i<=B;j++){
      v=i*prime[j];mark[v]=0;
	  if(i%prime[j]==0)break;
	}
  }
  return;
}
ll qmul(ll a,ll b,ll c){
  /*b%=c;if(b<0)b+=c;
  ll res=0;
  while(b){
    if(b&(ll)1)res=(res+a)%c;
    a=(a+a)%c;
    b>>=(ll)1;
  }*/
  return a*b%c;
}
ll ksm(ll a,ll b,ll c){
  ll res=(ll)1;
  while(b){
    if(b&(ll)1)res=qmul(res,a,c);
    a=qmul(a,a,c);
    b>>=(ll)1;
  }
  return res;
}
bool MR(ll x,ll n){
  ll b=n-(ll)1,res=ksm(x,n-(ll)1,n);
  if(res!=(ll)1)return false;
  while(b){
    if(b&(ll)1)return true;
    b>>=(ll)1;
    res=ksm(x,b,n);
    if(res!=(ll)1&&res!=n-(ll)1)return false;
    if(res==n-(ll)1)return true;
  }
  return true;
}
bool is_prime(ll n){
  int i;
  if(n==1ll)return 0;
  for(i=0;i<12;i++)if(n==(ll)test[i])return true;
  for(i=0;i<12;i++)if(MR(test[i],n)==false)return false;
  return true;
}
ll gcd(ll a,ll b){
  if(!b)return a;
  return gcd(b,a%b);
}
ll exgcd(ll a,ll b,ll &x,ll &y){
  if(b==0){x=1;y=0;return a;}
  ll d=exgcd(b,a%b,x,y);
  ll t=x;
  x=y;y=t-a/b*x;
  return d;
}
ll inv(ll a,ll b){
  ll x,y;
  exgcd(a,b,x,y);
  return x;
}
ll ABS(ll x){if(x>=0)return x;else return -x;}
inline ll Rand(){return (ll)rand()*rand()*rand()*rand();}
inline point point_add(point P,point Q,ll a,ll mod){
  point R;ll s;
  if(P.y==-1)return Q;
  if(Q.y==-1)return P;
  ll xp=P.x,yp=P.y,xq=Q.x,yq=Q.y;
  if(xp==xq&&yp!=yq){R.x=R.y=-1;return R;}
  if(xp==xq){
    ll d=gcd(2*yp,mod);
    if(d>1){R.x=d;R.y=-1;return R;}
	s=(ll)3*qmul(xp,xp,mod) %mod;
	s=(s+a)%mod;
    s=qmul(s,inv(2*yp,mod),mod);
  }
  else{
    ll d=gcd(ABS(xq-xp),mod);
    if(d>1){R.x=d;R.y=-1;return R;}
    s=qmul(yq-yp,inv(xq-xp,mod),mod);
  }
  R.x=(qmul(s,s,mod)-xp-xq)%mod;
  if(R.x<0)R.x+=mod;
  R.y=(qmul(s,xp-R.x,mod)-yp)%mod;
  if(R.y<0)R.y+=mod;
  return R;
}
inline point point_mul(point P,ll k,ll a,ll mod){
  point R;
  R.x=-1;R.y=-1;
  while(k){
    if(k&(ll)1){
	  R=point_add(R,P,a,mod);
	  if(R.x!=-1&&R.y==-1)return R;
	}
    P=point_add(P,P,a,mod);
    if(k!=1&&P.x!=-1&&P.y==-1)return P;
    k>>=(ll)1;
  }
  return R;
}
inline ll read(){
  ll sum=0;char ch=getchar();while(ch<'0'||ch>'9')ch=getchar();
  while(ch>='0'&&ch<='9')sum=sum*(ll)10+(ll)ch-(ll)48,ch=getchar();
  return sum;
}
void print(ll x){
  if(x>=(ll)10)print(x/(ll)10);
  printf("%d",(int)(x%(ll)10));
  return;
}
inline void ECM(ll n){
  //print(n);printf("\n");
  if(is_prime(n)){++cnt;p[cnt].v=n;p[cnt].k=1;return;}
  int i,j;ll a,x,y;
  for(i=1;i<=200;i++){
	a=Rand()%n;if(a==0)a=1;
    //x=Rand()%n;if(x==0)x=1;
    //y=Rand()%n;if(y==0)y=1;
	x=0;y=1;
	/*ll val=qmul(a,a,n);
	val=qmul(val,a,n);
	val=qmul(4,val,n);
	if((val+27)%n==0)continue;*/
	point P,Q,R;
    P.x=x;P.y=y;
    for(j=1;j<=tot;j++){
	  int w=(ld)log(B)/log(prime[j]);w=w+1;
	  ll v=ksm(prime[j],w,n);
	  R=point_mul(P,v,a,n);
	  if(R.x!=-1&&R.y==-1)break;
	  P=R;
	}
	if(R.x==1||R.x==n)continue;
	if(R.x!=-1&&R.y==-1){
	  ECM(R.x);ECM(n/R.x);
	  return;
	}
  }
  return;
}
bool cmp(node x,node y){return x.v<y.v;}
int main(){
	srand(time(0));
	init();
	int i,T;
	T=read();
	while(T--){
	  N=read();cnt=0;
	  ECM(N);
	  sort(p+1,p+cnt+1,cmp);
	  if(cnt==1)printf("Prime");
	  else print(p[cnt].v);
	  printf("\n");
	}
	return 0;
}
2025/1/30 21:21
加载中...