Dixon's & Quadratic Sieve

Sya_Resory · 2022-08-31 22:13:56 · 个人记录

~~开个坑。~~ 现在只差诈骗一个人过来帮我写代码了！

Prologue

两种基于 Smooth Number 理论的质因数分解算法。其思想都是寻找一个模 N 意义下二次剩余的非平凡平方根。

在需分解质因数的数 N 较大时优于 Pollard Rho，但是其实在 OI 界不怎么有用（？）

Smooth Numbers

称一个数是 y\text{-smooth} 的，当且仅当这个数的所有质因子小于 y。

Smooth Numbers 主要应用于一些质因数分解算法和求解离散对数算法中。

Lemma（Smooth Number 的密度）

令 \psi(x,y) 为 1\sim x 中所有 y\text{-smooth} 的数的个数，u\approx\frac{\log x}{\log\log x}，则有 \frac{\psi(x,x^{1/u})}x\approx u^{-u}。

证明比较麻烦。给出一个比较宽松的界的证明：

令 y^u=x，p_1,\dots,p_k 为小于 y 的素数，那么有 k\approx \frac y{\ln y}。那么对于所有 e_1,\dots,e_k 满足 \sum e_i\le u，都有 \prod p_i^{e_i} 为小于 x 的 y\text{-smooth} 的数。插板法可得 \psi(x,y)\ge \binom{k+u}u\ge(\frac ku)^u\approx \frac x{(u\ln y)^u}。于是有 \frac{\psi(x,x^{1/u})}x\ge(u\ln y)^{-u}。

Dixon's random square method

Dixon's 和 Quadratic Sieve 都基于以下思想：若 \alpha^2\equiv \beta^2\pmod N，且 \alpha\not\equiv\pm\beta，则 \gcd(\alpha+\beta,N) 为 N 的一个非平凡因子。考虑如何构造出符合条件的 \alpha 和 \beta。

考虑随机一个 \alpha 并令 \gamma=\alpha^2\bmod N，那么这个 \gamma 大概率是比较 \text{smooth} 的（证明在复杂度分析里）。考虑定一个阈值 y，令 p_1,\dots,p_k 为小于 y 的质数。如果 \gamma 是 y\text{-smooth} 的，并且质因数分解 \gamma=\prod p_i^{e_i} 后所有 e_i 都是偶数，那么可以取 \beta=\prod p_i^{e_i/2}。当 \gamma\perp N 时，由于 N 质因数分解后的每个因子 p^k 都对应 x^2\equiv\gamma 的两个根，因此 \gamma 在模 N 意义下有 2^{\omega(N)} 个平方根。于是，\alpha\equiv\pm\beta 的概率只有 \frac 1{2^{\omega(N)-1}}，这是可以接受的概率。

可惜并不是所有 \gamma 都满足 e_i 为偶数。考虑随机选取多个 \alpha_i，令 \gamma_i=\alpha_i^2\bmod N=\prod p_j^{e_{i,j}}。将指数看作向量，即令 v_i=(e_{i,1}\bmod 2,\dots,e_{i,k}\bmod 2)。要找到和为 0 的若干个向量，就相当于在 k 维空间中找出若干个线性相关的向量。因此只需要找到 k+1 个 \alpha 并高斯消元即可。

假设找到了 I\subset\{1,\dots,k+1\} 使得 \sum_{i\in I}v_i=0，那么我们只需要令 \alpha=\prod_{i\in I}\alpha_i，\gamma=\prod_{i\in I}\gamma_i=\prod p_i^{e_i}，\beta=\prod p_i^{e_i/2} 即有较大概率找出 N 的一个非平凡因子。

代码不知道为什么随机生成数的过程寄了（/oh），所以就不放完整的了。放个主要部分代码：

inline void Solve(u128 N) {
    for(;;) {
        int cnt = 0;
        for(;cnt <= pcnt;) {
            u128 x = randint(N),y = Mul(x,x); vec z;
            // check smoothness
            for(int i = 1;i <= pcnt;i ++) for(;!(y % pri[i]);y /= pri[i]) z[i] = !z[i];
            if(y == 1) a[++ cnt] = x,v[cnt] = z;
        }
        if(gauss()) return ; // gaussian elimination
    }
}

Lemma（\gamma=\alpha^2 \bmod N 大概率比较 smooth）：

令 k 为比阈值 y 小的质数个数，p 为第 k 大的质数，r=\lfloor\log_pN\rfloor。那么：
S=\sum_{\alpha<N,\gamma=\alpha^2\bmod N}[\gamma\text{ is }y\text{-smooth}]\ge\frac{k^{2r}}{(2r)!}

Proof：

考虑这样一个集合：T=\{\alpha\mid\alpha=\prod_{i=1}^kp_i^{e_i},\sum e_i=r\}，即所有质因数次数之和为 r 的 y\text{-smooth} 数。~~（我也不知道为什么文献上只数次数之和为 r 的）~~

令 l=\omega(N)，那么每个模 N 意义下的二次剩余有 2^l 个根。考虑质因数分解 N=\prod_{i=1}^lp_i^{f_i}，并将每个 \alpha 根据是否是模 p_i^{f_i} 意义下的二次剩余分为 2^l 类。显然同一类的数乘起来一定是模 N 意义下的二次剩余。将这些类中的数与 T 的交记作 T_i(1\le i\le2^l)。

每两个同类中的数相乘可以生成一个平方。由于每个平方根有 2^l 个根，且每个根至多被计算 \binom{2r}r 次，因此有：

S\ge\frac{2^l}{\binom{2r}r}\sum|T_i|^2\ge\frac1{\binom{2r}r}\left(\sum|T_i|\right)^2\ge\frac{\binom{k+r-1}{r}^2}{\binom{2r}{r}}\ge\frac{k^{2r}}{(2r)!}

得证。大概是挺松的一个下界。

算法的时间复杂度就是高斯消元的复杂度和 check smoothness 的复杂度。高斯消元显然 O(k^3)。随机数的复杂度根据上面的 Lemma 瞎折腾一下得到是 O(k^2\ln^{2r+2}N)。当 r=\sqrt{\frac{2\ln n}{\ln\ln n}}，y\approx \exp(2^{-1/2}\sqrt{\ln N\ln\ln N}) 时，复杂度为 O(\exp(\sqrt{\frac 92\ln N\ln\ln N})\text{polylog})。

The Quadratic Sieve

Dixon's 的复杂度瓶颈（之一）是判断 \gamma_i 是否为 y\text{-smooth}，因为我们使用的方法是朴素的试除。考虑用非随机化的方法生成 \alpha_i 并加速 check smoothness 的过程。

考虑这样一个数列：\{x^2-N\}，x=\lfloor\sqrt N\rfloor+i，1\le i\le B，其中 B 是某个设定的块长。我们要从这个数列中找出 k+1 个 y\text{-smooth} 数。这个序列的 smoothness 比较玄学，但我们可以假设其中一个数 smooth 的概率与随机选取的数一样。

类似 Eratosthenes 筛的过程，考虑把这个数列中所有小于 y 的质因子筛掉。假设我们当前筛到 p，p\nmid N。令 c 为模 p 意义下 N 的一个二次剩余，则 p\mid (x^2-N) 当且仅当 x\equiv \pm c\pmod p。于是 i\bmod p 仅有两种取值，枚举 p 的倍数从而找到所有 i，并将 x^2-N 除掉 p 的最高次项即可。这样筛掉 p_1,\dots,p_k 之后，变成 1 的就是 y\text{-smooth} 的数。

代码鸽了。

时间复杂度分析：

由定义，x 大致的取值范围是 0\sim X，其中 X\approx 2N^{1/2}。因为我们假设 Smooth Number 在该序列中的出现概率与随机类似，因此要生成 k+1 个 y\text{-smooth} 数就约需要 B\approx(k+1)\frac{X}{\psi(X,y)} 个数。由 Smooth Number 的密度可知当 y=X^{1/u} 时，B=O(X^{1/u}u^u)。

类似 Eratosthenes 筛的复杂度（O(n\log\log n)=O(n\text{polylog})），这里的筛法的复杂度约为 O(B\log\log y)=O(X^{1/u}u^u\text{polylog})。当 u\approx\sqrt{\frac{2\ln X}{\ln\ln X}} 时该式最小，此时 y\approx \exp(2^{-1/2}\sqrt{\ln X\ln\ln X})，复杂度为 O(\exp(\sqrt{\ln N\ln\ln N})\text{polylog})。

在优化 1 中，高斯消元的复杂度可以降至 O(k^2\log N)=O(\exp(\sqrt{\ln N\ln\ln N}))，于是这个就不是复杂度瓶颈了。

一些优化

可以参考 zzq 的文章，但是 zzq 把 Quadratic Sieve 的筛法一笔带过了（/fn/fn/oh），zzq 还是太强了。

注意到高斯消元的矩阵中，每行最多只有 O(\log N) 个非零元素，即这是个稀疏矩阵。用 Berlekamp-Massey 科技可以优化到 O(k^2\log N)。实战中用 bitset 优化效果（大概？）也差不多。
我们筛的数列是 \{x^2-N\}，x=\lfloor\sqrt N\rfloor+i，但其实可以同时筛 x=\lfloor\sqrt N\rfloor-i，常数大概能减半。
计算二次剩余要用 Cipolla's Algorithm，常数比较大。考虑定义一个数的 smoothness 的时候只用 p=4k+3 的素数（即该数的质因子均小于 y 且形如 p=4k+3），这样计算二次剩余的时候直接 c=N^{(p+1)/4} 即可。常数大概变为 \frac 12\sim\frac 14。
可以不只筛一个多项式。这个还是看 zzq 的文章吧，比较玄学，而且比较难写。