冒泡潜水
填一个坑。这个鸽了三年多了。
Preface
众所周知，Pollard-Rho 能以 \(O(n^{1/4}\log n)\) 的时间复杂度（通过一定优化能做到 \(O(n^{1/4})\)）找到合数 \(n\) 的一个非平凡因子。但事实上，我们有更快的方法做到这件事。
本篇文章将介绍一种亚指数级复杂度分解质因数的算法：二次筛法（Quadratic Sieve, QS）。
注：
- 本篇文章中的复杂度分析将不包含对大整数进行运算的开销。
- 本篇文章偏教程向，会混杂很多比较 trivial 的内容。
- 其实我自己写的代码基本是贺的 zzt 的提交并加以修改。
Smooth Numbers
数论中，一般认为素数的性质是比较好的。比较自然地，研究合数时，我们会希望研究那些“比较素”的合数。
对这样比较好的合数的刻画，容易会想到，使用素因子较少的合数。但这显然是荒谬的，因为质因数分解并不容易：对于大素数 \(p,q\)，合数 \(n=pq\) 很难被逆向表示为 \(p,q\)。
因而，我们转向另一个方向：研究素因子较小的合数，这便引出了 smooth numbers。
称一个数是 \(y\text{-smooth}\) 的，当且仅当其所有质因子不超过 \(n\)。
容易发现，smooth number 的一个很好的性质是，它们容易被分解。我们只需简单筛出不超过阈值 \(y\) 的素数 \(p_1,p_2,\dots,p_k\)，便能通过简单的试除法（因为 \(y\) 相对较小，这些开销是可以承受的）确定 \(n\) 的分解。
于是，我们会萌生一个想法：对于不 smooth 的 \(n\)，我们能否通过某些方法将其转化为比较 smooth 的数，再对其进行操作？想必，转化之后的操作应当是简单的。
事实上，数论中许多算法，如 Index Calculus 和本篇文章所要介绍的 Quadratic Sieve 都基于这样的转化。
Smooth Number 的密度
Smooth number 的密度由 Dickman 函数给出：
当 \(x\to\infty\) 时，有：
其中 \(\rho(u)\) 是 Dickman-de Brujin 函数。这表明，当 \(x\) 足够大时，smooth number 的密度不依赖于具体的 \(x,y\)，而只和 \(u=\frac{\ln x}{\ln y}\) 有关。
Dickman 证明了如下估计：
证明涉及微分学，这里不作详细展开。
Quadratic Sieve
Dixon's Random Square Method
我们的算法基于以下思想：若能找到 \(\alpha^2\equiv \beta^2\pmod N\) 且 \(\alpha\not\equiv\pm \beta\)，那么 \(\gcd(\alpha+\beta,N)\) 能够给出一个 \(N\) 的非平凡因子。问题在于，我们如何生成这样的 \(\alpha, \beta\)。
一个比较简单的想法是随机。我们随机这样的一个 \(\alpha_i\) 并计算 \(\gamma_i:=\alpha_i^2\bmod N\)，这是一个模 \(N\) 意义下的二次剩余。我们希望找到 \(\gamma_i\) 的一个非平凡平方根 \(\beta_i\)。假定 \(\gamma_i\perp N\)，那么 \(N\) 的每个质因子对应了 \(x^2\equiv \gamma_i\) 的两个根，因而找到的 \(\beta_i\equiv\pm \alpha_i\) 的概率只有 \(\frac{1}{2^{\omega(N)-1}}\le \frac 12\)。
理想情况下，若 \(\gamma_i\) 是比较 smooth 的，我们可以对其分解 \(\gamma_i=\prod p_j^{e_{i,j}}\)。若 \(\forall j\) 有 \(2\mid e_{i,j}\)，那么我们可以轻易地令 \(\beta_i=\prod p_j^{e_{i,j}/2}\)，这便得到一组 \(\alpha,\beta\)。
当然，我们运气应当不致于这么好，不仅 \(\gamma_i\) 是 smooth 的，而且 \(e_{i,j}\) 都是偶数。前者我们只能寄希望于概率，但是后者我们可以想想办法。考虑选取多个 \(\alpha_i\) 并计算 \(\gamma_i\)，再令向量 \(\boldsymbol v_i=(e_{i,1}\bmod 2,e_{i,2}\bmod 2,\dots,e_{i,k}\bmod 2)\in \mathbb Z_2^n\)。如果我们能找到一些线性相关的 \(\boldsymbol v_i\)，即存在指标集 \(I\) 使得 \(\sum_{i\in I}\boldsymbol v_i=\boldsymbol 0\)，那么可以取
一般的实现中，可以令 \(\beta_i:=\prod p_j^{\left\lfloor e_{i,j} / 2\right\rfloor}\)，那么就有：
假设选取的 smooth 阈值为 \(B\)，那么我们就是有若干个 \(\pi(B)\) 维的向量。只要我们找到 \(\pi(B)+k\) 个向量，我们就能通过高斯消元找出 \(k\) 组 \(\alpha,\beta\)，每组至少有 \(\frac 12\) 概率成功分解，那么总的分解成功概率就是 \(1-2^{-k}\)。
复杂度分析
我被一篇傻逼文献骗了好久
显然该方法的复杂度依赖于阈值 \(B\) 的选取。为方便计算，令 \(k = \pi(B),u=\frac{\ln N}{\ln B}\)。
\(\gamma_i=\alpha_i^2\bmod N\) 是模 \(N\) 的二次剩余，它大约在 \(1\sim N\) 中均匀分布，因此我们可以用前面的 smooth number 密度公式估算它 smooth 的概率，即 \(u^{-u}\)。我们需要找到 \(O(k)=O\left(\frac{B}{\ln B}\right)\) 个向量。每次尝试需要 \(k\) 次试除。最后，我们需要进行高斯消元。故可得到：
这种式子谁爱算谁算，我扔给 AI 了。算出来 \(T\) 最小时有：
The L-notation
细心的读者容易发现上面的 \(B\) 和 \(T\) 有相似的结构。事实上，这种结构的渐进表示在 smooth number 相关算法中相当常见。
一般地，我们定义 L-notation：
该函数的增长量级主要由 \(\alpha\) 控制：
- \(\alpha=1\) 时，\(L_n[1,c]=O\left(n^{c+o(1)}\right)\)（指数级复杂度）
- \(\alpha=0\) 时，\(L_n[0,c]=O\left((\ln n)^{c+o(1)}\right)=O(\text{polylog})\)
- \(\alpha\in(0,1)\) 时，\(L_n[\alpha,c]\) 表示亚指数级复杂度。
于是我们可以重新表示上面结果中的 \(B,T\)：
读者可以自行尝试证明该算法复杂度的另一种方法：设 \(B=L_N[1/2,c]\) 并带入以上复杂度式子，解出 \(c\)。
Optimized Dixon
朴素的 Dixon 方法使用了随机选取的 \(\alpha\)，因而 \(\gamma:=\alpha^2\bmod N\) 的 smoothness 似乎没有保证。我们会希望使用更好的 \(\alpha\) 来优化这一点。
一个比较显然的想法是选择 \(\alpha_i=\left\lceil \sqrt N\right\rceil+i\)。当 \(i\) 不太大时，这样得到的 \(\gamma_i\sim2i\sqrt N\) 比较小，应当会比较 smooth。
假设我们需要 \(M\) 个这样的 \(\alpha_i\)，那么应当有 \(\gamma_i \sim 2M\sqrt N\)，以这个数为上限套用之前的 smooth number 密度公式，计算得到，当 \(B=L_N[1/2,1/2]\) 时，有 \(T\approx L_N[1/2,3/2]\)。
QS
以上的算法在 \(N\approx 10^{30}\) 时就略显吃力。我们当然不满足于此。
注意到上面算法的复杂度瓶颈在于，对于每个 \(\alpha_i\)，我们都必须通过试除法确定 \(\gamma_i\) 是否是 smooth 的。由于选取合适的 \(\alpha_i\) 需要大量的尝试，这样反复试除的开销会变得不可接受。考虑如何减少试除次数。
我们都知道求质数可以用筛法，其实求 smooth 数也可以。更进一步，我们还可以筛多项式（因为 \(p\mid f(c)\Leftrightarrow p\mid f(kp+c)\)），比如一般的埃氏筛就是在筛 \(f(i)=i\)（如果跳过偶数可看作在筛 \(f(i)=2i+1\)）。
考察这样一个多项式：\(f(x)=x^2-N\)。当 \(i\) 较小时，有 \(f\left(\left\lceil\sqrt N\right\rceil+i\right)=\alpha_i^2\bmod N\)。我们希望筛一下它来优化上面的 Optimized Dixon 中的试除法。
类似埃氏筛地，我们枚举小质数 \(p\)。若 \(p^k\mid f(i)\) 就将 \(p^k\) 从 \(f(i)\) 的值中除去，筛完一遍变成 \(1\) 的就是 smooth 的 \(f(i)\)。若我们有 \(c\) 是 \(N\) 在模 \(p\) 意义下的二次剩余，那么当且仅当 \(x\equiv \pm c\pmod p\) 时会有 \(p\mid f(x)\)，枚举每个 \(x\) 并筛去 \(p\) 即可。
优化
对 QS 的改进方案中，提升较大的是选择更好的多项式去筛，这点留待后文详叙。这里讲几个小优化。
- \(f(i)\) 和 \(f(-i)\) 可以一起筛。注意消元的时候要把负号消掉。
- 注意到 \(p\mid f(i)\) 当且仅当 \(N\) 是模 \(p\) 意义下的二次剩余，因此在筛小素数时可以算一下 Legendre 符号 \(\left(\frac Np \right)\)，如果不为 \(1\) 这个素数就是没用的。
- 如果 \(f(i)\) 除掉小质数还剩一个不太大的数（一般的阈值是 \(B^2\)），可以把剩下这个值扔到
std::unordered_map里。如果后面找到的 \(f(j)\) 剩下同样的数，那么可以使用 \(f(i)f(j)\)，这个剩下的值就被平方了，乘到 \(\beta\) 里即可。（Lenstra 等提出，数据范围足够大时，可以尝试将除剩下的数分解为两个不大于 \(B^2\) 的数的乘积，并放到图上求解[1]） - 为了节省空间，我们在筛 \(f\) 时通常无法直接算出 \(\beta_i\) 和 \(\boldsymbol v_i\)，而需要对筛出的 smooth 的 \(f\) 再进行一次试除。这就使得筛 \(f\) 时的一堆大整数除法有点意义不明。可以考虑使筛 \(f\) 的条件变松：我们计算其对数，如果 $\log f(i) $ 减去 $ \sum_{p_j\mid f(i)}\log p_j$ 比较小，那么这个 \(f(i)\) 大概会比较 smooth，再对其进行试除即可。
- 朴素的高斯消元会成为算法的复杂度瓶颈。注意到 \(\boldsymbol v_i\) 中至多有 \(O(\log N)\) 个数非零，这说明我们最终要进行高斯消元的矩阵是一个稀疏矩阵。可以使用 Block Wiedemann 或 Block Lanczos 算法（Montgomery 给出过一种为素因数分解特化的在 \(GF(2)\) 上的 Block Lanczos[2]）将消元部分复杂度降至 \(O(k^2\log N)\)。实现中可以使用
std::bitset优化的在线高斯消元，性能不会有明显劣势。
复杂度分析
同 Dixon，我们需要筛出 \(O(k)\) 个数，每个数成功概率为 \(u^{-u}\)，筛每个数的成本为 \(\sum_{p\le B}p^{-1}=O(\log\log B)\)。高斯消元复杂度 \(O(k^2\log N)\)。故有：
当 \(T\) 最小时，有：
MPQS
Multiple Polynomial Quadratic Sieve，简称 MPQS （标题冒泡潜水的由来）。
朴素 QS 的问题是，当 \(N\) 较大时，我们需要筛的 \(f(i)\) 个数（设为 \(M\)）会较大，而 \(f(i)\sim 2M\sqrt N\) 也会较大，这使得 \(f\) 较 smooth 的猜想不再成立。为了解决这个问题，Montgomery 提出：我们可以筛多个多项式。
具体而言，我们希望用一个线性函数 \(ax+b\) 替代之前的 \(x\)，即考虑 \((ax+b)^2-N\)。如果 \(b^2\equiv N\pmod a\)，那么 \((ax+b)^2-n\) 可以写成 \(af(x)=a(ax^2+2bx+c)\) 的形式，其中 \(b^2-ac=N\)。为了方便处理，我们设 \(a\) 是一个 smooth 数的平方（比如 \(p^2\)，也有文献建议取一系列小素数乘积的平方），那么我们只需要去筛 smooth 的 \(f(x)\)。
假设我们只去筛 \(|x|\le M\) 的 \(x\)。与之前同样地，我们希望使 \(|f(x)|\) 较小以使得其较 smooth。根据中学数学知识，我们会希望抛物线对称轴 \(x=-\frac ba\) 靠近所筛区间 \([-M,M]\) 的中点 \(x=0\)，这要求 \(b\) 较小。另外，我们希望区间端点和中点处 \(|f(x)|\) 值相近，即 \(|f(M)|\approx |f(0)|\)，解出 \(a\approx\frac{\sqrt{2N}}{M}\)。
实现中，我们可以从 \(\frac{(2N)^{1/4}}{M^{1/2}}\) 开始枚举小素数 \(p\)，令 \(a=p^2\)，取 \(b\) 为 \(b^2\equiv N\pmod a\) 的较小解，再算出 \(c\)。然后按照 QS 的方法筛 \((ax+b)^2-N\) 即可。
MPQS 相比于 QS 复杂度不变，但有较大常数优势。zzq 说有分析认为提升倍数是 \(\frac 12\sqrt{\ln N\ln\ln N}\)，我不知道是怎么算出来的。
HMPQS
The Hypercube variation of Multiple Polynomial Quadratic Sieve，简称 HMPQS。
Peralta 提出，MPQS 中切换多项式还要重算一遍 \(b,c\) 之类的，比较麻烦。考虑生成一些有联系的多项式来消灭掉这个代价。其实和 Hypercube 没啥必然联系，Peralta 为了容易理解随便举的例子
我们发扬光大一下 MPQS 中的一个细节。令 \(t\) 为 \(n\) 个素数的乘积（\(t:=\prod_{i=1}^nq_i\)），且 \(t\sim N^{1/4}\)。假设 \(a=t^2\)，\(b\) 为方程 \(b^2\equiv N\pmod {t^2}\) 的一个解。显然这里 \(b\) 有 \(2^n\) 个解。我们希望这样所生成出的 \(2^n\) 个 \((ax+b)^2\) 有比较好的性质。
首先考虑如何表示 \(b\)。由 CRT，只需解出 \(\alpha_i^2\equiv N\pmod {q_i^2}\)，那么就可以令 \(b=\sum \delta_i\alpha_i\beta_i\)，其中 \(\delta_i=\pm1\)，\(\beta_i\) 是 CRT 算出的系数。\(\alpha_i\) 有两个解，此处我们选择使得 \(\gamma_i:=\alpha_i\beta_i\bmod t^2\le\frac{t^2}2\) 的解。
然后考虑如何遍历所有的 \(b\)。令数列 \(s_1,\dots,s_{2^n}\) 遍历所有的 \(b\) 解。我们希望遍历过程中每次只更改一个 \(\delta_i\)，那么可以有递推关系：
其中 \(\mu_i=\pm 1\)，代表 \(\delta_{k_i}\) 的改变。Peralta 的叙述中把这个递推关系理解为 \(n\) 维 Hypercube 上的一个 Hamilton 路径，\(\mu_i\) 代表第 \(k_i\) 维坐标的改变；当然其实还有其他理解方法，比如相信大家都知道格雷码，原理是类似的。\(\omega_i\) 可以理解为对 \(t^2\) 取模的一个修正系数，使得 \(s_{i+1}\in(0,t^2)\)。由之前对 \(\alpha_i\) 的选择，有 \(\omega_i\in\{\pm 1,0\}\)。一般的实现中可以直接取 \(\omega_i\equiv 0\)，此时仍有 \(s_i\in(-nt^2,nt^2)\)。
当然我们并不关心具体值。在筛法中，\(s_i\) 作为多项式的系数，我们关心其模 \(p\) 意义下的值。Peralta 提出，为了方便计算 \(s_i\bmod p\)，我们可以打表一下其差分。令：
注意到 \(\mu_i\) 有两种取值，\(\omega_i\) 有三种取值，\(\gamma_i\) 有 \(n\) 种取值，那么 \(\Delta_i\) 至多有 \(6n\) 种取值（如果固定 \(\omega_i=0\) 则只有 \(2n\) 种）。这个可以打表算出来。那么转移就只需极小的代价。
现在可以构造我们需要的多项式了。对 MPQS 中的多项式略作变换，除掉 \(a\) 这个平方因子，就得到我们构造的多项式：
我们仍须验证这样生成的 \(f_i(x)\) 不太大，以保证其 smooth 的概率。假设我们要筛的区间是 \(x\in[-M,M]\)，我们取 \(t\le\frac{N^{1/4}}{\sqrt M}\)。由于 \(s_i^2<t^4<N\)，通过简单放缩可以得到 \(\frac{s_i^2-N}{t^2}<0\) 且 \(\left|\frac{s_i^2-N}{t^2}\right|\le M\sqrt N\)。二次项 \(x^2t^2\le M^2\left(\frac{N^{1/4}}{\sqrt M}\right)^2=M\sqrt N\)。一次项 \(2s_ix\le 2M\sqrt N\)。因此我们大概有 \(f_i(x)=O(M\sqrt N)\)，这说明其值不太大，也就保证其大概比较 smooth。
在筛法中，令 \(p\mid f_i(x)\)，容易得到 \(x\equiv(-s_i\pm\sqrt N)t^{-2}\pmod p\)。这就得到筛 \(f_i(x)\) 的两个起点：
\(s_i\bmod p\) 容易由上述递推算出。如果预先算出 \(\sqrt N\bmod p, t^{-2}\bmod p,2\sqrt Nt^{-2}\bmod p\)，就可以很快算出 \(D_i,E_i\)。这就节省了切换多项式的时间。
还剩下最后一个细节：如果 \(p\mid t\)，那么以上同余式的推导不成立。事实上，当 \(p\mid t\) 时，当且仅当 \(x\equiv -\frac{(s_i^2-N)/t^2}{2s_i}\pmod p\) 时有 \(p\mid f_i(x)\)，即某个 \(x\) 使得 \(p\mid f_i(x)\) 的概率减半为 \(1/p\)。这提示我们，\(t\) 的因数尽量不要选取 smooth check 所使用的小素数，以增大 \(f_i(x)\) smooth 的概率。实现中，我们可以选择素数 \(q\in[B,2B]\) 来生成 \(t\)，这样大约有 \(n\approx\log_B N\approx2\sqrt\frac{\log N}{\log\log N}\)。
我并没有写这个的代码
References:
-
Prime Numbers: A Computational Perspective (zzq 那篇的 reference)
-
A Quadratic Sieve on the n-Dimensional Cube (Peralta)
浙公网安备 33010602011771号