二次筛法分解小记

command_block · 2020-02-22 09:34:23 · 个人记录

暂时搁置

我们现在有一个合数n,要想办法找到一个非平凡因子。

怎么构造呢?我们可以枚举b,然后查看n-b^2是否是完全平方数即可。

我们可以考虑枚举x,计算Q(x),然后查看是否是平方数,这做样并不会比费马分解聪明多少。

但是,我们如果能够找到若干个Q(x),使得其乘积是完全平方数,也能构造出一组解。

设我们有x_1x_2...x_m=u;\ Q(x_1)Q(x_2)...Q(x_m)=v^2

注意到在同余系下有x-n=x\pmod n

现在,u和v就是一组解了,我们最好祈祷它们是非平凡的……

现在问题来了,我们怎么找这些解?

对于若干个Q(x),我们并不能使用惊为天人的构造,只能够将其看做一般的数来处理。

选取子集乘积为平方数并不是一个简单的问题,方案数是指数级的。

注意到,我们如果将数质因数分解,每个素数的指数相加为偶数的话,就能得到平方数了。

问题在于,质因数分解并不是一件很轻松的事情。

我们观察得到,大多数数的素因子都是很小的,而分解的代价主要集中在大的素因子上。

对于那些最大素因子较大的Q(x),我们完全可以抛弃(反正也难以配对)。只分解那些小的。

这样会不会造成能够分解的Q(x)很多呢?事实上是不会的,原因见下文。

换句话说,我们钦定一个限制B,每次只用B以内的素数试除。

我们怎么分解呢?一个个试除也可以,但是太慢了。

需要判断 p|x^2-n 对哪些x生效。

等价于x^2=n\pmod p,我们需要求解n在模p下的二次剩余。

这里p都很小,可以直接暴力打表。

正常的求解方法请见二次剩余小记

我们求出x的通解kp+c之后,就类似埃氏筛把Q(kp+c)都提取一通。

最后查看每个Q是否被提取尽,如果没有被提取尽就丢掉。

我们的x应该如何取值呢?为了方便,我们要使得Q(x)的值为正,那么x>\sqrt{n}

为了方便改设Q(x)=(x+\lfloor\sqrt{n}\rfloor)^2-n,这样子x的取值就从1起步了。

能不能把x的值取很大呢?事实上,越大的数,含有大质因子的概率就越大(无法分解),我们要尽量让Q(x)的值小一些,来增加分解成功的概率。

现在问题就是 : 给出若干个01向量,要求选一些xor起来为0(线性相关).

我们知道,d维向量空间里面,最多能构造d个线性无关的向量。

我们任意找出d+k个向量,必然产生k组线性相关的关系。

大力高斯消元即可,采用类似线性基的方法来消元,记得使用bitset优化。

下面给出一些收集情况,选取一个素数集合,从N开始逐个分解,收集到|P|个能够分解的数。

事实上,观察(x+\lfloor\sqrt{n}\rfloor)^2-n=x^2+2x\lfloor\sqrt{n}\rfloor+\lfloor\sqrt{n}\rfloor^2-n,

它十分接近x(x+2\sqrt{n}),以至于误差不超过O(\sqrt{n})

可以想象,它不会出现质因子很大的情况,而是比较smooth

下面我们对Q(x)进行分解测试。注意要把注定无法分解的质数丢掉(二次剩余)

设分解次数为R,我们的复杂度是O(P^3/w+RloglogR)

一般来讲B越大,R越小,但是消元越慢,调参是一个玄学的事情。

这个算法在同样范围的数据中运行速度都差不多,所以可以对着随机数据调参。

为了节省内存,可以分段分解。

测试记录,可以看到喜提最劣解。

留一个测试代码,跑起来奇慢无比……

#include<algorithm>
#include<cstring>
#include<cstdio>
#include<bitset>
#include<cmath>
#include<ctime>
#define i128 __int128
#define ll long long
using namespace std;
template<class T>
void read(T &x){
  char c;
  for (c=getchar();c<'0'||c>'9';c=getchar());
  for (x=0;c<='9'&&c>='0';c=getchar())x=x*10+(c&15);
}
char gu[100];
template<class T>
void write(T x){
  int tot=0;
  if (!x)putchar('0');
  while(x){
      gu[++tot]=x%10+'0';
      x/=10;
  }for (int i=tot;i;i--)
    putchar(gu[i]);
}
i128 n1,n2,N,sft;
inline i128 mul(i128 a,i128 b){
  i128 r=((__float128)a/N*b+0.4);
  r=a*b-r*N;
  return r<0?r+N:r;
}
#define MaxP 1250
struct BufData
{int p,c;}b[MaxP*2];
int p[MaxP],bcnt,pcnt;
void getPrime()
{
  int lim=5*pow(N,0.05)*pow(log2((double)N),0.333);printf("%d\n",lim);
  ll sav=sqrtl(N);
  for (int i=2;pcnt<=lim;i++){
    bool flag=1;
    for (int j=2;j*j<=i;j++)
      if (i%j==0)flag=0;
    if (flag){
      int t=N%i;flag=0;
      for (int j=0;j<i;j++)
        if (j*j%i==t){
          b[++bcnt]=(BufData){i,((j-sft)%i+i)%i};
          flag=1;
        }
      if (flag)p[++pcnt]=i;
    }
  }
}
i128 gcd(i128 a,i128 b){
  if(a==0)return b;
  if(b<0)b=-b;
  i128 t;
  while(b){t=a%b;a=b;b=t;}
  return a;
}
#define F(x) (((x)+sft)*((x)+sft)-N)
#define MaxT 270000
int Blo,scnt;
i128 Q[MaxT];
bitset<MaxP> s[MaxP],sa[MaxP],ts,tsa;
bool e[MaxP];
struct Data
{int c[MaxP];i128 x,p;}a[MaxP];
int tc[MaxP];
i128 powM(i128 a,int t)
{
  i128 ans=1;
  while(t){
    if (t&1)ans=mul(ans,a);
    a=mul(a,a);
    t>>=1;
  }return ans;
}
void tfac(bitset<MaxP> &s)
{
  memset(tc,0,sizeof(int)*(pcnt+2));
  i128 u=1,v=1;
  for (int i=1;i<=scnt;i++)
    if (s[i]){
      u=mul(u,a[i].x);
      v=mul(v,a[i].p);
      for (int j=1;j<=pcnt;j++)
        tc[j]+=a[i].c[j];
    }
  for (int i=1;i<=pcnt;i++)
    if (tc[i])v=mul(v,powM(p[i],tc[i]/2));
  i128 d=gcd(N,u-v);
  if (d!=1&&d!=N){
    u=d;v=N/d;if (u>v)swap(u,v);
    write(u);printf(" ");write(v);
    printf("\n%d\n",clock());
    exit(0);
  }
  d=gcd(N,u+v);
  if (d!=1&&d!=N){
    u=d;v=N/d;if (u>v)swap(u,v);
    write(u);printf(" ");write(v);
    printf("\n%d\n",clock());
    exit(0);
  }
}
void addset(i128 num,i128 x,i128 sp)
{
  ts.reset();tsa.reset();
  for (int i=1;i<=pcnt;i++){
    tc[i]=0;
    while(num%p[i]==0)
      {num/=p[i];tc[i]++;}
    if (tc[i]&1)ts[i]=1;
  }tsa[++scnt]=1;
  a[scnt].x=x;a[scnt].p=sp;
  memcpy(a[scnt].c,tc,sizeof(int)*(pcnt+2));
  bool flag=1;
  for (int i=pcnt;i;i--)
    if (ts[i])
        if (e[i]){ts^=s[i];tsa^=sa[i];}
      else {
        s[i]=ts;sa[i]=tsa;
        e[i]=1;flag=0;break;
      }
  if (flag){tfac(tsa);scnt--;}
}
void upd(int c,int p)
{
  for (int i=c;i<Blo;i+=p){
    if (Q[i]<2)continue;
    while(Q[i]%p==0)Q[i]/=p;
    if (Q[i]==1)addset(F(i),i+sft,1);
  }
}
void getPset()
{
  for (;;sft+=Blo){
    for (int i=0;i<Blo;i++)Q[i]=F(i);
    for (int i=1;i<=bcnt;i++){
      upd(b[i].c,b[i].p);
      b[i].c=((b[i].c-Blo)%b[i].p+b[i].p)%b[i].p;
    }printf("%d %d\n",scnt);
  }
}
int main()
{
  //n1=10007ll; n2=50021ll;
  //n1=1000000000039ll; n2=500000000023ll;
  n1=1000000000000037ll; n2=500000000000057ll;
  //n1=973806882626147ll;  n2=981039245468473ll;
  //n1=1000000000000000003ll; n2=4000000000000000631ll;
  //read(N);
  N=n1*n2;
  sft=sqrtl(N);
  Blo=max(1000,min((int)pow(N,0.2),262000));
  getPrime();
  printf("%d\n",clock());
  getPset();
  return 0;
}

当然还有几个优化。

对于Q(x)值为负的,也不是不可利用,只是要把-1当做素数提取出来。

具体的写法可以查看代码。这样能够提速40%左右。

然后,分解完毕后某些1<Q(x)<p_{max}^2,这肯定剩下了一个大质数。

对于单独的这样一个Q(x),我们无计可施,毕竟向矩阵里面加入一行成本太高。

但是,我们如果收集到了多个Q(x)分解成p,我们就可以把它们两两乘起来,这样p就被平方了。

这样又能够提速40%左右。而且因为分解概率增大所以矩阵大小不需要开那么大,节省了空间。

最终得到的代码,虽然快了许多,单组最差500ms,但还是打不过优秀的Prho,距离我们的最终目标更是遥远。

测试记录

#include<algorithm>
#include<cstring>
#include<cstdio>
#include<bitset>
#include<cmath>
#include<ctime>
#include<map>
#define i128 __int128
#define ll long long
using namespace std;
template<class T>
void read(T &x){
  char c;
  for (c=getchar();c<'0'||c>'9';c=getchar());
  for (x=0;c<='9'&&c>='0';c=getchar())x=x*10+(c&15);
}
char gu[100];
template<class T>
void write(T x){
  int tot=0;
  if (!x)putchar('0');
  while(x){
      gu[++tot]=x%10+'0';
      x/=10;
  }for (int i=tot;i;i--)
    putchar(gu[i]);
}
i128 n1,n2,N,sft;
inline i128 mul(i128 a,i128 b){
  i128 r=((__float128)a/N*b+0.4);
  r=a*b-r*N;
  return r<0?r+N:r;
}
#define MaxP 1250
struct BufData
{int p,c;}b[MaxP*2];
int p[MaxP],bcnt,pcnt;
void getPrime()
{
  int lim=8*pow(N,0.05)*pow(log2((double)N),0.333);
  //printf("%d\n",lim);
  ll sav=sqrtl(N);
  for (int i=2;pcnt<=lim;i++){
    bool flag=1;
    for (int j=2;j*j<=i;j++)
      if (i%j==0)flag=0;
    if (flag){
      int t=N%i;flag=0;
      for (int j=0;j<i;j++)
        if (j*j%i==t){
          b[++bcnt]=(BufData){i,j};
          if (j)b[++bcnt]=(BufData){i,i-j};
          flag=1;break;
        } 
      if (flag)p[++pcnt]=i;
    }
  }
}
i128 gcd(i128 a,i128 b){
  if(a==0)return b;
  if(b<0)b=-b;
  i128 t;
  while(b){t=a%b;a=b;b=t;}
  return a;
}
#define F(x) (((x)+sft)*((x)+sft)-N)
#define MaxT 270000
int Blo,scnt;
i128 Q[MaxT];
bitset<MaxP> s[MaxP],sa[MaxP],ts,tsa;
bool e[MaxP];
struct Data
{int c[MaxP];i128 x,p;}a[MaxP];
int tc[MaxP];
i128 powM(i128 a,int t)
{
  i128 ans=1;
  while(t){
    if (t&1)ans=mul(ans,a);
    a=mul(a,a);
    t>>=1;
  }return ans;
}
void check(i128 x)
{
  x=gcd(N,x);
  if (x!=1&&x!=N){
    i128 u=x,v=N/x;
    if (u>v)swap(u,v);
    write(u);printf(" ");write(v);
    //printf("\n%d\n",clock());
    exit(0);
  }//else puts("!!!");
}
void tfac(bitset<MaxP> &s)
{
  memset(tc,0,sizeof(int)*(pcnt+2));
  i128 u=1,v=1;
  for (int i=1;i<=scnt;i++)
    if (s[i]){
      u=mul(u,a[i].x);
      v=mul(v,a[i].p);
      for (int j=1;j<=pcnt;j++)
        tc[j]+=a[i].c[j];
    }
  for (int i=1;i<=pcnt;i++)
    if (tc[i])v=mul(v,powM(p[i],tc[i]/2));
  check(u-v);check(u+v);
}
void addset(i128 num,i128 x,i128 sp)
{
  ts.reset();tsa.reset();
  if (num<0){num=-num;ts[0]=1;}
  for (int i=1;i<=pcnt;i++){
    tc[i]=0;
    while(num%p[i]==0)
      {num/=p[i];tc[i]++;}
    if (tc[i]&1)ts[i]=1;
  }tsa[++scnt]=1;
  a[scnt].x=x;a[scnt].p=sp;
  memcpy(a[scnt].c,tc,sizeof(int)*(pcnt+2));
  bool flag=1;
  for (int i=pcnt;i>=0;i--)
    if (ts[i]){
      if (e[i]){ts^=s[i];tsa^=sa[i];}
      else {
        s[i]=ts;sa[i]=tsa;
        e[i]=1;flag=0;break;
      }
    }
  if (flag){tfac(tsa);scnt--;}
}
void upd(int c,int p)
{
  for (int i=c;i<Blo;i+=p){
    if (Q[i]<2)continue;
    while(Q[i]%p==0)Q[i]/=p;
  }
}
int aaa;
map<int,i128> si;
void secwork()
{
  for (int i=0;i<Blo;i++){
    Q[i]=F(i);
    if (Q[i]<0)Q[i]=-Q[i];
  }
  for (int i=1;i<=bcnt;i++)
    upd((b[i].c-sft-b[i].p)%b[i].p,b[i].p);
  for (int i=0;i<Blo;i++){
    if (Q[i]==1)
      addset(F(i),i+sft,1);
    else if (Q[i]<p[pcnt]*p[pcnt]){
      if (si.find(Q[i])!=si.end()){
        i128 sav=si[Q[i]];
        addset((F(i)/Q[i])*((sav*sav-N)/Q[i]),(i+sft)*(sav),Q[i]);
        aaa++;
      }else si[Q[i]]=i+sft;
    }
  }
}
void getPset()
{
  sft=sqrtl(N);
  for (i128 det=0;;det+=Blo){
    sft+=det;secwork();sft-=det;
    sft-=det+Blo;secwork();sft+=det+Blo;
    //printf("%d %d %lld\n",scnt,aaa,(long long)det);
  }
}
int main()
{
  //n1=10007ll; n2=50021ll;
  //n1=1000000000039ll; n2=500000000023ll;
  //n1=1000000000000037ll; n2=500000000000057ll;
  //n1=973806882626147ll;  n2=981039245468473ll;
  //n1=100000000000000003ll; n2=400000000000000013ll;
  //N=n1*n2;
  read(N);
  Blo=max(1000,min((int)pow(N,0.2),260000));
  getPrime();//printf("%d\n",clock());
  getPset();
  return 0;
}

前面我们都是在做一些常数上边边角角的优化,现在我们来一些有理有据的东西。

我们所利用的式子就是x^2-y^2=(x+y)(x-y),有没有什么其他的分解方式呢?

当然有,比如a^3-b^3=(a-b)(a^2+ab+b^2),但是这要求我们构造两个完全立方数,难度更大了。

看来,我们要利用二次式来产生分解。

只能整理出x(x+c)=y(y+d)\pmod n,这要求我们在x(x+c)这种东西上面构造。这比纯粹的平方数还要困难得多……

这也就是说,我们还是要坚持x^2-y^2=(x+y)(x-y)不放,相信费马nb!

此外,就是我们的函数Q(x)=(x+\lfloor\sqrt{n}\rfloor)^2-n了,当x取到较大值的时候,Q(x)能够成功分解的概率会明显降低,大概和1/x有关。

我们如果有多个函数,能在绝对值小的范围里产生多一些数,能够分解的概率就大大提升了。

一个容易想到的方法是Q(x)=(x+\lfloor\sqrt{cn}\rfloor)^2-cn,但是尴尬的是,c为合数时会跟原先的Q(x)产生大量重复值,所以我们只能让c为质数。

在分解的时候,原来我们可以针对单一的平方式,我们可以舍去无二次剩余的|P|,现在由于多个Q,我们不能再舍去,常数大了一倍。