二次筛选法

二次筛选（Quadratic Sieve）算法是一个整数分解算法，在实际用途中为已知第二快的方法（目前第一快为普通数域筛选法）。但对于大约 100 位数以内的整数，它仍然是最快的算法，而且比起普通数域筛选法来说简洁得多。

这是一个通用的整数分解算法，意即其运算时间完全取决于欲分解的整数本身位数的大小，而不是在于特殊结构或特性。

二次筛选法是由卡尔·帕梅朗斯在1981年所发明，并作为理查德·施罗佩尔的线性筛法之改良版。

## 基础目标

此算法试图去建立一个模 $n$ （$n$为欲分解的数）下的平方同余，这往往即是$n$的因数分解。算法有两个阶段：“数据收集”，在此阶段收集可能可以找到一个平方同余的资料；以及“数据处理”，它把所有收集的数据放进一个矩阵里，并解决、获得一个平方同余数。

数据收集的阶段可以很轻易地使用多个处理器去平行化。但数据处理阶段需要大量的内存，并且在多个运算节点之间有效地平行化相当困难，也可能每个运算节点的内存不够足以储存整个阵列。而Block Wiedemann算法可以使用在一些可以保存阵列的系统。

要找到一个平方同余，一个较天然的方法便是随机挑选数字，将其平方，并希望模$n$之后的非负余数是一个完全平方数。例如：$80^2$ 模 $5959 = 441$，同时也是$21^2$。

这种方法对很大的$n$值而言，可以找到一个同余的平方数的情况很罕见。但是当真的找到了一个时，在大多数情况下，同余数为非平凡解而整数分解便完成了。这大致上即是费马因式分解法（Fermat's factorization method）的核心。

而二次筛选法改良自狄克森因式分解法。

一般来说，二次筛选法的执行时间（去质数分解一个整数$n$时）为

$$e^{(1 + o(1))\sqrt{\ln n \ln\ln n}} =L_n\left[1/2,1\right]$$

参见 L-符号。

上式常数 e 为自然对数之底数。

## 解决方法

令 $x$ 模 $y$ 表示为 $x$ 除以 $y$ 之后所剩的余数。为了分解整数，费马因式分解法牵涉到需要寻找一个数字 $a$（$n^{1/2}<a<n-1$），使得$a^2\pmod{n}$是一个完全平方数。但这些$a$值相当难找到。二次筛选法包括对于好几个$a$值去计算了$a^2\pmod{n}$，然后在$a$值与$a^2\pmod{n}$的集合中找到一个子集，当中的元素之乘积为完全平方数。而产生出一个平方同余。

例如：$41^2$ 模 $1649 = 32$、$42^2$ 模 $1649=115$ 以及 $43^2$ 模 $1649$ 为 $200$。在这些数字（32、115、200）当中皆无完全平方数，但存在一乘积$32\times 200=6400=80^2$是一个平方数。模1649 之后，这个乘积$32\times 200=(41^2)\times(43^2)=(41\times 43)^2=114^2$（因为 $41\times 43$ 模 $1649 = 114$）。$32\times 200 = 80 ^2$ 的观察因而给出了一个平方同余：$114^2\equiv 80^2$ (模 $1649$)。

但是，如何将以下问题解决呢？“给予一组数字，找到一个子集使其乘积是平方数。”该解决方案使用了指数向量的概念。而指数向量，例如根据算术基本定理，504 可分解为 $2^3 3^2 5^0 7^1$。

这表示可以借由指数向量 $(3,2,0,1)$，代表 $2,3,5,7$ 在因式分解的指数值。490 会同样可分解为指数向量 $(1,0,1,2)$。将这些数字相乘相当余把其指数向量的对应值一一相加：$504\times 490$ 得一向量 $4,2,1,3$。

有一些数字为平方数，其满足每个在其指数向量的各个数字为偶数。例如，向量 $(3,0,0,1)$、$(1,2,0,1)$ 之和为 $(4,2,0,2)$，因此 56 乘以 126 是一个平方数。找寻一个平方数只需对于向量里数字之的奇偶性之知识，所以有可能将整个向量简化为模 2 的形式并作模 2 下的加法: $(1,0,0,1)+(1,0,0,1)=(0,0,0,0)$。在实作中，这相当地有效率，因其可以表示为一位元集（bitsets）且模 2 之加法将变为位元运算互斥或（XOR）。

于是此问题变化为：“给予一个 0, 1向量构成的集合，找到一个子集，其中所有向量之和为模 2 的零向量。”而这是一个线性代数的问题；且该解答为线性相依的。

线性相依是线性代数中的一个定理：当有比每个向量中含有的元素还要多的向量时，这种相依关系必然存在。而它可以被高效率地找到，例如：把所有向量一列一列地排在一个矩阵里，然后使用高斯消去法。比起实数来说，此方法尤其容易套用到模 2 后的整数上。而此算法所需的平方数即是那些向量所对应的数字之积。

然而，纯粹地只去将一堆随机数字平方并模 $n$ 会产生很大量的、不同的质因数，也因此会产生出很长的向量以及一个非常大的矩阵。解决方法是去找到一些特别的数字 $a$，使得 $a^2$ 模 $n$ 之值只由很小的质因数组成（它们都是光滑数）。此种数字很难被找到，但是若仅使用光滑数将可以保持向量和矩阵之尺寸更小、更容易处理。

而二次筛选法使用一种之后会提及名为“筛法”（sieving）的技巧去找寻光滑数，也就是此算法的命名由来。

## 算法

总地来说，二次筛算法基本有以下主要的步骤：

- 选择一个光滑数之上界 $B$。以质数计数函数 $\pi(B)$ 表示小于 $B$ 的质数之数量，其将控制之后向量的长度以及需要的向量之数量。
- 使用筛法找到 $\pi(B)+1$ 个数字 $a_i$ 使得 这样 $b_i=(a_i^2$ 模 $n)$ 为一个 $B$-光滑数。
- 将 $b_i$ 作质因数分解生成一个指数向量（每个数字都要模 $2$）。
- 套用线性代数的概念找到一个子集，其中的每个向量之和为一零向量。把这些向量所对应的 $a_i$ 相乘、$b_i$ 相乘并模 $n$：便得到一个 $B$-光滑数 $b^2$.
- 现在，我们得到方程式 $a^2=b^2 ($ 模 $n)$因为从步骤 4 得到 $(a^2$ 模 $n)$ 的两个平方根，一个即是对整数 $b^2$ 取平方根，也就是 $b$；另一个则是步骤 4 得到的 $a$ 本身。
- 因此现在，我们掌握了所需的恒等式：$(a+b)(a-b)\equiv0 \pmod n$ $(a+b)(a-b)\equiv0 \pmod n$ $(a+b)(a-b)\equiv0\pmod{n}$。计算 $n$ 与 $a-b$ （或是 $a+b$）的 GCD （最大公因数，Greatest Common Divisor）。此将产生 $n$ 的其中一个因数，尽管有可能是一个平凡（trivial）因数 (即为 $n$ 或是 1)。如果此因数是平凡的，使用不同的线性相依或是不同的 $a$ 值去再次尝试。

本文的剩余部分将解释这个基本算法的细节和延伸。

## 二次筛法（QS）如何最佳化找寻同余

二次筛选法试图找到一整数对 $x$ 和 $y(x)$ （其中 $y(x)$ 为 $x$ 的函数）其满足比 $x^2\equiv y ($模 $n)$ 还要弱得多的条件。它选择一些质数作为一集合作为“因数基底”，并试图找到 $x$，使得 $y(x)=x^2$ 模 $n$ 之值的质因数只会在此因数基底。此时可称 $y$ 值：对于此因数基底是光滑的。

$y(x)$ 的其中一个值之因式分解（为因数基底的一部分），跟 $x$ 一起，被称为“关系”（relation）。二次筛选法借由采取接近 $n$ 的平方根之 $x$ 值，以加速寻找这类“关系”的过程。这将确保 $y(x)$ 会较小，因而具有更大的可能性是光滑的。

$$y(x)=\left(\left\lceil\sqrt{n}\right\rceil+x\right)^2-n\hbox{ (where }x\hbox{ is a small integer)}$$

$$y(x)\approx 2x\left\lceil\sqrt{n}\right\rceil$$

这意味着， $y$ 在 $2x\sqrt{n}$ 的数量级上.。然而，这也意味着 $y$ 的增长幅度与 $x$ 乘以 ($n$的平方根) 成正比。

另一个可以增加光滑的可能性是，即是单纯地增大因数基底的大小。然而，比起因数基底的质数数量，至少找到一个光滑的“关系”还是必要许多，其确保存在一个线性相依。

### “部分关系”以及循环

即使对于某些“关系”来说， $y(x)$ 并非光滑的。但如果两个 $y$ 刚好是由因数基底以外的相同质数之乘积，也可能可以合并这两个部分“关系”，以形成一个完整的“关系”。注：此形同于因数基底的扩展。

例如：如果因数基底为$\{2,3,5,7\}$和 $n=01$，存在“部分关系”（partial relations）：

$${21^2\equiv 7^1\cdot 11\pmod{91}}$$

$${29^2\equiv 2^1\cdot 11\pmod{91}}$$

将上面两式乘在一起：

$${(21\cdot 29)^2\equiv2^1\cdot7^1\cdot11^2\pmod{91}}$$

并将等号两边皆乘上 $\left(11^{-1}\right)^2$ 模 $91$。而 $11^{-1}$ 对 $91$ 取模为 $58$，所以：

$$(58\cdot 21\cdot 29)^2\equiv 2^1\cdot7^1\pmod{91}$$

$$14^2\equiv 2^1\cdot7^1\pmod{91}$$

即产生了一个完整的“关系”。这样的一个完整的“关系”（借由结合“部分关系”所获得的）称为循环。有时候，从两个“部分关系”形成的循环，可以直接导向一个平方同余，但是此情况非常罕见。

### 借由筛选来检查光滑度

有好几种方法可以 $y$ 值们的光滑度。最直觉的是借由试除法，尽管这样会增加数据收集阶段的运行时间。

另一个方法较能被接受的方法是椭圆曲线因式分解（ECM）。

而在实作中，称为筛选的方法比较会被经常使用。设 $f(x)$ 为多项式 $f(x)=x^2-n$ 我们得：

$$f(x)=x^2-n$$

$$f(x+kp)=(x+kp)^2-n$$

$$f(x+kp)=x^2+2xkp+(kp)^2-n$$

$$f(x+kp)=f(x)+2xkp+(kp)^2\equiv f(x)\pmod{p}$$

因此解决出 $f(x)\equiv 0 ($模 $p)$对于某个 $x$ 值，将产生出一整个序列，当中的每个数值 $y~(y=f(x)$ 皆可被 $p$ 整除。此问题便是对某个质数取模下找到一个平方根，对其存在着高效率的算法，例如谢克斯–托内里算法的。（这便是二次筛法的名称来由： $y$ 是一个 $x$ 的二次多项式且筛选过程中的运算类似埃拉托斯特尼筛法。）

筛选一开始将一个大阵列 $A[]$ 每个“元”（entry）的每个字节设为零。对于每一个 $p$，去解决模$p$ 下的二次方程式并得到两个根 $\alpha$ 和 $\beta$，然后在每个 $y(x)=0$ 模 $p$“元”之中加入一个近似于 $\log(p)$ 之值……也就是 $A[kp+\alpha]$ 和 $A[kp+\beta]$。为了辨识数字是否可被因数基底中的质数之平方所整除，解决几个模 ( $p$ 的小次方) 下的二次方程式也是必要的。

在因数基底的尾端，任何 $A[]$ 有包含一个值超过大约为 $\log(x^2-n)$ 的临界值，将会对应到一个 $y(x)$ 值，其由因数基底的部分组成。那些包含了确定 $y(x)$ 可以被哪些质数整除的资讯已经遗失掉了，但是因为其只包含一些小的因数，而且已知有很多优良的算法可以去分解那些已知只有小因数的数字。例如小质数的试除法、SQUFOF、波拉德 ρ，以及ECM，以上是经常作为一起使用的方法。

基本上很多 $y(x)$ 值都会是可行的，因此因式分解过程的尾声不需要是完全可信的；通常此过程大约有 5% 的输入会出现异常，此时需要做少量的额外筛选。

## 基本筛选的例子

以下例子将演示没有作对数优化或是质数次方的标准的二次筛法。令要分解的数为 $N=15347$，因此平方根 $N$ 无条件进位为124。由于 $N$ 很小，因此基本的多项式即足够了： $y(x)=(x+124)^2-15347$。

### 数据收集

因为 $N$ 为小数字，所以只需 4 个质数。满足在模 $p$ 下 15347 有一平方根的前 4 个质数 $p$ 为 2、17、23 以及 29(换句话说，对这些质数来说，15347是一个模这些数字的二次剩余)。这些质数将是筛选的基础。

现在我们要建造出我们的筛选 $V_X$从 $Y(X) = (X + \lceil\sqrt{N}\rceil)^2 - N = (X+124)^2-15347$并开始对基底里每个质数进行筛选，以下选择筛出 $0\leq X<100$ 的那些 $Y(X)$：

$$\begin{align}V &= \begin{bmatrix} Y(0) & Y(1) & Y(2) & Y(3) & Y(4) & Y(5) & \cdots & Y(99) \end{bmatrix} \\ <pre> & =\begin{bmatrix} 29 & 278 & 529 & 782 & 1037 & 1294 & \cdots & 34382 \end{bmatrix}\end{align}$$

下一步即是去作筛选的动作。对于我们的质数基底 $\lbrace 2, 17, 23, 29\rbrace$中的每一个质数 $p$ 值去解决以下方程式：

$$Y(X) \equiv (X + \left\lceil\sqrt{N}\right\rceil)^2 - N \equiv 0 \pmod{p}$$

找到阵列 $V$ 之中可被 $p$ 所整除的那些“元”。

对于 $p=2$ 解出 $(X + 124)^2 - 15347 \equiv 0 \pmod{2}$ 得到了 $X \equiv \sqrt{15347}-124 \equiv 1 \pmod{2}$。

所以，从 $X=1$ 开始每次 $+2$，每个“元”可被 2 整除。把那些元除以 2 之后得到：

$$V = \begin{bmatrix} 29 & 139 & 529 & 391 & 1037 & 647 & \cdots & 17191 \end{bmatrix}$$

同理，对于剩下的质数 $p$ $\lbrace 17, 23, 29\rbrace$ 方程式$X \equiv \sqrt{15347} - 124 \pmod{p}$ 也解决了。

值得注意的是，对于每一个 $p>2$，因为有两个模平方根，因此得到 2 个线性方程式。

$$\begin{align} <pre> X & \equiv \sqrt{15347} - 124 & \equiv 8 - 124 & \equiv 3\pmod{17} \\ & & \equiv 9 - 124 & \equiv 4\pmod{17} \\ X & \equiv \sqrt{15347} - 124 & \equiv 11 - 124 & \equiv 2\pmod{23} \\ & & \equiv 12 - 124 & \equiv 3\pmod{23} \\ X & \equiv \sqrt{15347} - 124 & \equiv 8 - 124 & \equiv 0\pmod{29} \\ & & \equiv 21 - 124 & \equiv 13\pmod{29} \\ </pre> \end{align}$$

每个方程式 $X \equiv a \pmod{p}$ 导致 $V_x$ 从 $x=a$ 和之后每一次递增一个 $p$ 值的那些项次皆可被 $p$ 整除。把 $V$ 中的 $a$、$a+p$、$a+2p$、$a+3p$等等的位置除以 $p$ , 如此对于每个在基底中的质数可以找到为相异质数的乘积（一次方）之光滑数。

$$V = \begin{bmatrix} 1 & 139 & 23 & 1 & 61 & 647 & \cdots & 17191 \end{bmatrix}$$

在 $V$ 之中的值等于一的那些“元”皆对应到一个光滑数。因为 $V_0$, $V_3$， $V_{71}$ 等于一，因此对应到：

$X+124$｜$Y$｜因数

$124$｜$29$｜$2^0\cdot17^0\cdot23^0\cdot29^1$

$127$｜$782$｜$2^1\cdot17^1\cdot23^1\cdot29^0$

$195$｜$22678$｜$2^1\cdot17^1\cdot23^1\cdot29^1$

### 矩阵处理

由于根据 $Y \equiv Z^2 \pmod{N}$的性质我们已经找到平滑数 $Y$，而算法接着的剩余部分等同于狄克森因式分解法中的任何变体。

将方程式中的一个子集里的指数乘积

$$\begin{align} 29 &= 2^0 \cdot 17^0 \cdot 23^0 \cdot 29^1 \\ 782 &= 2^1 \cdot 17^1 \cdot 23^1 \cdot 29^0 \\ 22678 &= 2^1 \cdot 17^1 \cdot 23^1 \cdot 29^1 \\ \end{align}$$

转为一个矩阵形式 （在模 2 下）得到以下方程式：

$$S \cdot \begin{bmatrix} 0 & 0 & 0 & 1 \\ 1 & 1 & 1 & 0 \\ 1 & 1 & 1 & 1 \end{bmatrix} \equiv \begin{bmatrix} 0 & 0 & 0 & 0 \end{bmatrix} \pmod{2}$$

此方程式可由零空间（null space）的概念所给出一个解，为：

$$S = \begin{bmatrix}1 & 1 & 1 \end{bmatrix}$$

因此三个方程式的乘积产生了一个平方数（模 $N$ 之下）：

$$29 \cdot 782 \cdot 22678 = 22678^2$$

以及

$$124^2 \cdot 127^2 \cdot 195^2 = 3070860^2$$

所以此算法找到了

$$22678^2 \equiv 3070860^2 \pmod{15347}$$

测试其结果得到 $\gcd(3070860 - 22678, 15347)=103$，为 15347 的一个非平凡因数，而另一个为149。

而以上恰好显示出，二次筛法只适用于 $n$ 值较大时。对于例如像 15347 这类的小数字，此算法显得过犹不及。试除法或是波拉德 ρ都可以在少量许多的计算之下找到一个因数。

## 倍数多项式

在实际用途上，有许多相异的多项式用在 $y$ 上，因为仅仅一个多项式通常不足以产生出对于因数基底的光滑数对 $(x,y)$。使用的多项式使用必须要有一个特别形式，因为它们需要为模 $n$. 下的平方数。多项式必定会与原始的 $y(x)=x^2-n$ 有类似的形式：

$$y(x)=(Ax+B)^2-n \qquad A,B\in\mathbb{Z}$$

假设 $B^2-n$ 是 $A$ 的一个倍数，则 $B^2-n = AC$ 且多项式 $y(x)$ 可以被写作 $y(x) = A\cdot(Ax^2+2Bx+C)$。而如果 $A$ 为一个完全平方数，则只需考虑 $(Ax^2+2Bx+C)$的部分。

此方式（称为 MPQS，倍数多项式二次筛选法（Multiple Polynomial Quadratic Sieve））非常适合平行运算，因为每一个处理因式分解的处理器可以单纯的给入 $n$、因数基底以及多项式的集合，且直到运算完多项式之前都不须跟中央处理器作任何传输。

## 大质数

### 单一的大质数

如果在除以所有小于 $A$ 的因数之后，剩余的数字（余因子）小于 $A^2$，那么这个余因子必为质数。实际上，借由对于余因子去排序“关系”表，则它可以添加进因数基底里。如果 $y(a)=7\times 11\times 23\times 137$ 且 $y(b)=3\times 5\times 7\times 137$, 则 $y(a)\times y(b)==3\times 5\times 11\times 23\times 7^2\times 137^2$。此可以降低以上完整执行因式分解的筛选阵列之“元”的临界值。

### 更多的大质数

甚至更进一步去降低临界值，并且使用一个高效处理将 $y(x)$ 之值分解为一些更大的质数之积（ECM 适合处理这样子的东西）可以找到因数大多在因数基底，但有两个甚至三个大质数的“关系”。循环的寻找过程因此允许一个共享好几个质数的“关系”集合，合并成为单一的“关系”。

## 实际例子下的参数

为了展示在一个有包含多个多项式以及大质数优化下的实作方式去跑实际例子，会有的典型参数选取，将一个 267 位元的半质数输入进 msieve 中，产出了以下的参数：

- 试除因数分解截止于：27 位元
- 筛选区间（对于每个多项式）：393216（12 个大小为 32768 的区块）
- 光滑数之上界：1300967 (共 50294 个质数)
- 对于多项式 A 的系数之因数数量：10 (见上面倍数多项式条目)
- 大质数之上界：128795733 (26 位元) (见上面大质数条目)
- 光滑数的发现数：有 25952 为直接筛出，另外的 24462 为借由合并那些有大质数的数字所得出
- 最终矩阵的大小：50294 × 50414，借由过滤法减少到 35750 × 35862
- 非平凡的线性相依之发现数：15
- 总执行时间 (在 1.6 GHz UltraSparc III 上)：35 分 39 秒
- 最大内存使用量：8 MB

## 整数分解的纪录

直到发现普通数体筛选法（number field sieve, 简称 NFS）之前，二次筛法（QS）曾是已知渐近最快的通用整数分解算法。现在，伦斯特拉椭圆曲线因式分解具有跟 QS 有相同的渐近运行时间（在 $n$ 由两个相同大小级别的质数相乘所得的情况下），但在实际情况中，QS 速度更快，因为它采用的是单精度浮点数操作而不是椭圆曲线所使用的高精度计算。

在 1994 年的 4 月，RSA-129 的因数分解借由 QS 完成了。其为一个由两个大质数相乘的129 位数数字一个因数为 64 位长而另一个为 65 位。此因数分解的因数基底包含了 524339 个质数。数据收集阶段花了 5000 个 MIPS 年，其完成于互联网上的分散式计算。数据收集总量为 2 GB。数据处理花了45个小时在 Bellcore ( 现为 Telcordia 科技公司) 的 MasPar （大规模的平行化）超级电脑。这曾是最大的、借由通用算法的公开分解，直至 NFS 被用于分解 RSA-130，于 1996 年 4 月 10 日 完成。所有自此以后分解的 RSA号码 皆使用 NFS。

目前 QS 的纪录是 $2^{803}-2^{402}+1$ 的一个 135 位数长之余因子，其为 $2^{1606}+1$的一个 Aurifeuillian因数，在 2001 年分解为 66 位以及 69 位数长的质因数。

## 实作

- PPMPQS and PPSIQS
- mpqs
- SIMPQS 是由 William Hart 编写的自我初始化（self-initializing）的倍数多项式二次筛选法的快速实作。其提供大质数的优化变体，并使用 Jason Papadopoulos' block Lanczos 程式码于线性代数阶段.。SIMPQS 可以使用 qsieve 指令在 SageMath 电脑代数套件上存取，或是原始来源里下载。SIMPQS 被优化用于 Athlon 和 Opteron 机器上，但仍可在最常见的 32、 64 位元的结构上运行。而其完全是由 C 语言编写而成的。
- 由 Dario Alpern 所提供的 factoring applet，其在特定状况之下会使用二次筛选法。
- PARI/GP 电脑代数套件包含着自我初始化（self-initializing）的倍数多项式二次筛选法的一个实作并有着大质数的优化变体。其源自于 Thomas Papanikolaou 以及 Xavier Roblot 的一个编写给 LiDIA 计划的筛选法。自我初始化的方法是基于一个 Thomas Sosnowski 的一篇论文上的一个点子。
- 一个二次筛选法的变体开放于 MAGMA 电脑代数套件。其基于 1995 年 Arjen Lenstra 一个使用于他自己的“透过电子邮件分解整数”计划的一次实作。
- msieve，一个支援单个或双个大质数的倍数多项式二次筛选法之实作，由 Jason Papadopoulos 所编写。源代码以及 Windows 的二进制档案皆是公开的。
- YAFU，由 Ben Buhrow 所编写，与 msieve 相似但是对于现今大多的处理器来说更快。其使用 Jason Papadopoulos' block Lanczos 程式码。源代码以及 Windows、Linus 的二进制档案皆是公开的。
- Ariel，一个用于教学用途的二次筛选法 Java 简易实作。
- java-math-library 包含着也许是编写于 Java 最快的二次筛选法（PSIQS 4.0 的后继者）。
- Java QS，一个开源的 Java 计划包含着 QS 的基本实作。由 Ilya Gazman 于 2016 年 2 月 4 日所释出。
- C Quadratic Sieve – 使用C语言编写的二次筛选法算法的公共领域实现。于2022年发布，支持对最多330位的数进行批量分解，并以JSON或CSV格式输出结果。

## 外部链接

- Reference paper "The Quadratic Sieve Factoring Algorithm" by Eric Landquist