博客详情 - oirush

admin 的博客
二次筛法分解小记
admin LV 10 SU @ 2026-7-10 9:03:42

更新于 2025/12/12 00:04:05 作者

command_block

暂时搁置

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

$n\leq 10^{30}$ ,时间限制 $\texttt{100ms}$ ,允许__int128.

费马分解法

这是一个naive的分解方法。

如果能构造 $n=a^2-b^2$ ,就能由平方差公式得到 $n=(a+b)(a-b)$

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

$a$ 在最坏情况下需要从 $1$ 枚举到 $\sqrt{n}$ ,并没有什么复杂度上的提升。

甚至对于随机情况,由于最小因子通常较小,试除法往往能够更快地得到结果。

即使是两种方法并用,也无法将复杂度变得低于 $O(\sqrt{n})$ 。

起点式子

费马分解法虽然没有成功,但也给我们一些启示 : 合数分解可以被转化成简洁的构造。

考虑这样的一个同余式 : $x^2=y^2\pmod n$ (即 $x^2-y^2=kn$ )

也就相当于, $y^2-x^2=(x+y)(x-y)$ 是 $n$ 的倍数。

那么 $x+y$ 和 $x-y$ 中必定有 $n$ 的一部分因子,我们可以使用 $gcd$ 来提取。

问题在于,可能恰好提取到 $n$ 本身,这就是无用功。

事实上,如果 $n$ 是奇合数,至少有一半的解是非平凡的。证明需要群论,超出了本文的范围,不再赘述。

也就是说,我们如果用随机的手法构造出一组解,就有 $\frac{1}{2}$ 的概率得到分解,那么我们期望只需要尝试 $O(1)$ 次。

怎么构造呢?

设 $Q(x)=x^2-n$

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

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

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

$v^2=Q(x_1)Q(x_2)...Q(x_m)=(x_1^2-n)(x_2^2-n)...(x_m^2-n)$

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

$=x_1^2x_2^2...x_m^2=u^2$

现在, $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|$ 个能够分解的数。

$P$	起始数	分解次数
$62$	$10^{9}$	$5876$
$62$	$10^{12}$	$110593$
$168$	$10^{9}$	$4145$
	$10^{12}$	$39419$
	$10^{15}$	$548451$
$489$	$10^{12}$	$27137$
	$10^{15}$	$205057$
	$10^{18}$	$1990822$

事实上,观察$(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)$ 进行分解测试。注意要把注定无法分解的质数丢掉(二次剩余)

$P$	$\sqrt{N}$ 的数量级	分解次数
$177$	$10^{11}$	$40000$
$283$	$10^{12}$	$1400000$
$605$	$10^{15}$	$6000000$
$1171$	$10^{18}$	$22000000$

设分解次数为 $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^2-y^2=(x+y)(x-y)$ 能够整理出 $x^2=y^2\pmod n$ ,所以我们就是要构造两个 $\!\!\pmod n$ 相等的完全平方数。

而其他带有一次项的分解,如 $(x+y+...)(x-y+...)=x^2+cx-y^2-dy$

只能整理出 $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$ ,我们不能再舍去,常数大了一倍。

$Q(x)=(x+\lfloor\sqrt{cn}\rfloor)^2-cn=x^2+2\lfloor\sqrt{cn}\rfloor+O(\sqrt{cn})=O(x\sqrt{cn})$,

可以看到,其绝对值的增加速度取决于 $\sqrt{c}$ ,实际上,我们可以记录每一个 $Q$ 的时效,然后决定使用哪个。

但是这样根本做不到比原来快,因为带的常数太大了。

考虑其他的方法,比如

二次筛法分解小记

状态

开发

支持

关于

二次筛法分解小记

状态

开发

支持

关于

还没有账户？

登录