二次筛法分解小记
command_block
2020-02-22 09:34:23
**暂时搁置**
我们现在有一个合数$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$都很小,可以直接暴力打表。
正常的求解方法请见[二次剩余小记](https://www.luogu.com.cn/blog/command-block/er-ci-sheng-yu-xiao-ji)
我们求出$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$ |
| $168$ | $10^{12}$ | $39419$ |
| $168$ | $10^{15}$ | $548451$ |
| $489$ | $10^{12}$ | $27137$ |
| $489$ | $10^{15}$ | $205057$ |
| $489$ | $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$越小,但是消元越慢,调参是一个玄学的事情。
这个算法在同样范围的数据中运行速度都差不多,所以可以对着随机数据调参。
为了节省内存,可以分段分解。
[测试记录](https://loj.ac/submission/750169),可以看到喜提最劣解。
留一个测试代码,跑起来奇慢无比……
```cpp
#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`,距离我们的最终目标更是遥远。
[测试记录](https://loj.ac/submission/750582)
```cpp
#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$的时效,然后决定使用哪个。
但是这样根本做不到比原来快,因为带的常数太大了。
考虑其他的方法,比如