19、大规模并行椭圆曲线分解算法的实现与优化

大规模并行椭圆曲线分解算法的实现与优化

1. 引言

分解算法的理论和实践研究对分析众多知名的公钥密码系统至关重要。分解算法通常分为通用算法和特殊用途算法。通用算法的预期运行时间仅取决于被分解数 (n) 的大小，而特殊用途算法的预期运行时间还与 (n) 的（未知）因子的性质有关。特殊用途算法包括试除法、Pollard 的 rho 方法、Pollard 的 (p - 1) 方法和椭圆曲线方法。

算法名称	运行时间特点
试除法	找到因子 (p) 的时间近似与 (p) 成线性关系
Pollard 的 rho 方法	近似与 (\sqrt{p}) 成线性关系
Pollard 的 (p - 1) 方法	近似与 (p - 1) 的最大质因子成线性关系
椭圆曲线方法	找到 (p) 的预期时间为 (O((\log n)^2L_n[\sqrt{2}]))，在最坏情况下 (p \approx \sqrt{n}) 时，方法的预期时间为 (L_n[1])，是关于 (n) 的次指数时间

从安全角度看，上述指数时间的方法通常无需担忧，但有时会采取预防措施，避免 Pollard 的 (p - 1) 方法意外成功。然而，椭圆曲线方法无法通过此类预防措施避免，因为只要该方法恰好找到一个接近 (p) 且仅含小因子的数，就有可能找到 (p)。由于该方法包含多个独立试验，所以总有发现 (p) 的可能性。

本文探讨了在单指令多数据（SIMD）机器上实现椭圆曲线方法的两种方案。尽管大规模并行在该领域并不新鲜，但据我们所知，此前尚未有人考虑过在 SIMD 机器上实现椭圆曲线方法。SIMD 机器相对多指令多数据（MIMD）机器成本较低，却能实现出色的椭圆曲线分解性能。

2. 硬件介绍

我们使用的 16K MasPar 是一款 SIMD 机器，大致由前端、阵列控制单元（ACU）和一个 (128 \times 128) 的处理单元（PE）阵列组成。通过掩码或条件语句可选择和更改 PE 阵列中的活动处理器子集，即活动集。

由于它是 SIMD 机器，指令按顺序执行，涉及并行数据的指令由活动集中的所有处理器同时执行，而 PE 阵列中的其他处理器则处于空闲状态。涉及单数据（即非并行数据）的指令在前端或 ACU 上执行，在本文描述中，前端和 ACU 作用相同。

每个 PE 每秒大约可执行 (2 \times 10^5) 次 32 位整数加法，可视为 0.2 MIPS 的处理器。每个 PE 拥有 64KB 内存，整个 PE 阵列则有 1GB 内存。每个处理器可与北、东北、东、东南、南、西南、西和西北方向的邻居进行通信，采用环形环绕方式。实际上，处理器可向这八个方向上任意距离的处理器发送数据，中间的处理器也可能获取传输数据的副本。此外，还有一个全局路由器允许任意处理器与其他处理器通信，但我们未使用该功能。

每个作业大小在 8K 到 1M 之间，反映了其使用的 PE 内存量。只有占用 PE 内存总和不超过 64K 的作业才会以轮询方式调度，每个作业执行 10 秒后被抢占，其他作业则需等待。这意味着作业不会从 PE 内存中交换出去。

我们的实现使用了 MasPar 并行应用语言 MPL，从我们的角度来看，它是 C 语言的简单扩展。

3. 椭圆曲线方法

椭圆曲线方法由多个独立试验组成。每次试验随机选择一个模 (n) 的椭圆曲线 (E) 以及与 (E) 相关的群 (G) 中的一个点 (x)。群 (G) 中的运算（以乘法表示）包含多个模 (n) 的整数加法、减法、乘法和求逆运算，每次运算时间复杂度为 (O((\log n)^2))。若在群运算过程中需要对某个满足 (\gcd(n, y) \neq 1) 的整数 (y) 进行模 (n) 求逆，群运算将失败，此时就找到了 (n) 的一个因子。

利用群运算将点 (x) 提升到一个巨大的幂 (k)，(k) 是小于某个界限 (B_1) 的所有素数幂的乘积。若该计算因群运算失败而无法完成，则此次试验成功，找到了 (n) 的一个因子；若成功计算出 (x^k \in G)，则试验失败。若 (p) 是 (n) 的最小素因子，且 (\psi_1 = L_p[\sqrt{m}])，则每次试验分解 (n) 的概率为 (L_p[-\sqrt{m}])，这意味着分解 (n) 所需的独立试验次数预计为 (L_p[\sqrt{m}])。一次试验的时间复杂度为 (O((\log n)^2B_1))，因此总预期运行时间为 (O((\log n)^2L_n[\sqrt{2}]))。

计算 (x^k \in G) 通常称为算法的第一阶段。在第二阶段，对所有介于 ([B_1, B_2]) 之间的素数 (q) 计算 (x^{kq} \in G)，这大约需要额外的 (\pi(B_2) - \pi(B_1)) 次群运算，其中 (\pi(b)) 表示小于等于 (b) 的素数个数。选择 (B_1) 和 (B_2) 使得两个阶段的时间大致相同似乎是接近最优的选择，此时 (B_2) 比 (B_1) 大一个数量级。

椭圆曲线方法的最优参数选择取决于 (n) 的未知因子 (p)，因此事先无法确定。常见的做法是先进行几次 (B_1) 较小的试验，若失败则增加 (B_1) 再进行几次试验，直到分解出 (n) 或放弃尝试。

我们在 16K MasPar 上并行运行了 16K 次试验，每个 PE 独立处理自己的椭圆曲线，使用单个随机种子唯一生成。为此，我们将 [4] 中使用的多精度整数运算适配到了 MasPar 上，使得每个 PE 能独立处理自己的扩展精度整数，与其他 PE （在活动集中）同时进行操作。这样，我们可以并行执行 16K 次扩展精度整数的基本运算（加、减、乘、商和余数），无需处理器间通信。

由于模 (n) 求逆运算速度较慢，特别是在 SIMD 机器上，其成本由需要最多迭代次数的处理器决定，因此我们记录群元素的分子和分母（模 (n)），而不进行求逆运算。但由于这些求逆运算可能导致 (n) 的分解，因此无法完全避免。我们定期计算所有分母的模 (n) 乘积（在 PE 阵列上使用 11 次模 (n) 乘法），并在前端（速度快得多）计算所得乘积与 (n) 的最大公约数 (g)。若 (g > 1)，则检查 PE 以确定有多少个 PE 找到了因子，若 (g) 不是素数，则进一步细化分解结果。通常，只有在存在多个非常小的因子（最多 10 位）时才会出现后一种情况；较大的因子通常仅由一个 PE 找到。

模 (n) 加法和减法在 SIMD 模式下容易实现高效运算，但模 (n) 乘法更具挑战性，这是因为模 (n) 余数计算的指令流更依赖于涉及的值。我们在整个程序中使用了 Montgomery 表示法，它允许进行与数据无关的模乘法，而不影响加法或减法运算。

尽管这些改进显著提高了我们的第一个实现的性能，但由于 0.2 MIPS 的 PE 处理器速度较慢，对于 100 位的 (n)，完成 16K 条曲线（(B_1 = 50000)，(B_2 = 10^6)）大约需要一天时间。使用这些参数几乎可以保证找到所有最多 28 位的因子，但这并非最优参数选择。要找到 28 位的因子，使用更少的曲线和更大的 (B_1) 会更好，而对应 16K 次试验的最优 (B_1) 超过 (10^7)。若要寻找 45 位的因子，这是一个不错的选择，但在该机器上需要过长的时间。

因此，这个并行椭圆曲线程序并不完全适合当前的 MasPar 机器。然而，对于未来的 SIMD 机器，如果 PE 的运行速度与当前工作站相当，那么 16K 次并行试验（匹配的界限）最多只需几天即可完成。

鉴于当前 PE 的速度，要在 16K MasPar 上获得更好的并行椭圆曲线程序，似乎只能将每条曲线的工作分配给 (r) 个 PE，这样可以减少试验次数，同时有望将每条曲线的速度提高相同的倍数，从而允许选择更接近 30 位因子最优解的参数。我们通过设计一种与上述方法完全不同的多精度整数运算实现了这一目标，下面将详细介绍。

4. 分块模运算

本节介绍一种适用于 SIMD 机器的替代整数运算方法，该方法本质上可将机器重新配置为处理器数量较少但速度更快的机器，至少对于模固定整数 (n) 的运算如此。

选择一个小整数 (b)，使得 PE 能够高效执行 (b) 位整数的算术运算。由于 PE 执行 32 位整数乘法并得到 64 位结果的效率最高，我们选择 (b = 30)，这样在加法时也不会出现溢出问题。设 (r) 是满足 (2^{b \cdot r} > n) 的最小整数，其中 (n) 是待分解的奇数。将 128 行 128 个 PE 中的每一行划分为 (u = \lfloor 128 / (r + 1) \rfloor) 个不相交的 (r + 1) 个连续 PE 的块，以及 (128 - (r + 1) \cdot u) 个空闲 PE。因此，共有 (128 \cdot u) 个块。活动集由块内的 PE 组成。

假设块内连续的 PE 编号从 0 到 (r)。若第 (i) 个 PE 包含一个 (b) 位整数 (v_i \geq 0)，则 (v_0, v_1, \cdots, v_r) 共同表示数 (v = \sum_{i = 0}^{r} v_i 2^{b \cdot i})。由于我们使用 ({0, 1, \cdots, n - 1}) 中的整数表示模 (n) 的剩余类，通常有 (v < n)，因此 (v_r = 0)。这个额外的 PE 用于模 (n) 乘法。

设 (v = \sum_{i = 0}^{r} v_i 2^{b \cdot i}) 和 (w = \sum_{i = 0}^{r} w_i 2^{b \cdot i}) 是两个模 (n) 的整数。为计算 (s = v + w)，(r + 1) 个 PE 计算 (s_i = v_i + w_i) 和 (c_i = \lfloor s_i / 2^b \rfloor)，然后第 (i) 个 PE 将 (c_i) 发送给相邻的第 (i + 1) 个 PE，最后计算 (s_i = s_i - c_i \cdot 2^b + c_{i - 1})。在极少数情况下，若某个 (s_i \geq 2^b)，则重复进位传播过程，直到所有 (s_i < 2^b)。这里需要注意的是，在 MasPar 上只需一条快速指令即可检查是否还有进位需要传播。为完成模 (n) 加法，我们使用类似的技术从 (s) 中减去 (n)；若 (s - n) 为负数，则 (s) 为最终结果，否则结果为 (s - n)。

虽然容易构造需要 (r) 次进位传播步骤的示例，使用深度为 (O(\log r)) 的进位传播树可以获得更好的最坏情况性能，但我们的简单方法在平均情况下更快，因为第二次进位传播步骤很少发生。

块内的模 (n) 乘法更为复杂。与第一个椭圆曲线实现一样，我们使用 Montgomery 表示法避免除法。设 (R = 2^{b \cdot r}) 是大于 (n) 的最小 (2^b) 幂。整数 (x) 模 (n) 的 Montgomery 表示法 (\tilde{x}) 是整数 (x \cdot R \bmod n \in {0, 1, \cdots, n - 1})。Montgomery 表示法下的数的加法和减法与普通模 (n) 加法和减法相同，按上述方法执行。然而，乘法比普通模 (n) 乘法简单得多。设 (z = x \cdot y \bmod n)，则 (\tilde{z} = x \cdot y / R \bmod n)。我们可以通过以下步骤高效计算 (\tilde{z})：
1. 计算 (t = x \cdot y)。
2. 设 (v = \sum_{i = 0}^{r - 1} v_i 2^{b \cdot i})，并找到 (d) 使得 (d \cdot n \equiv -1 \bmod 2^b)（由于 (n) 是奇数，(d) 是定义明确的）。
3. 对于 (i = 0, 1, \cdots, r - 1)，依次将 (v) 替换为 (v + 2^{b \cdot i} \cdot n \cdot (v_i \cdot d \bmod 2^b))，其中 (v_i)（(i > 0)）是上一次迭代计算得到的 (v) 的 (2^b) 进制数字。
4. 注意，经过第 (j) 次迭代后，新的 (v_j) 为零，且新的 (v) 与旧的 (v) 模 (n) 同余。因此，最终的 (v_0) 到 (v_{r - 1}) 都为零，通过将结果 (v) 右移即可完成除以 (R) 的操作。结果可能大于等于 (n)，此时只需减去 (n) 一次即可使其小于 (n)。

Montgomery 乘法可以使用 (2r^2 + r) 次 (b) 位整数乘法完成。

将 (x) 和 (y) 的（普通）乘法迭代与模 (n) 除以 (R) 的操作合并后，在 PE 块中可以很容易地实现 Montgomery 表示法下的数的乘法。直接应用上述算法会得到一种块式模乘法，需要在 (r) 个 PE 上并行执行 3r 次乘法，而不是我们期望的 ((2r^2 + 1) / r = 2r + 1) 次乘法。如 [2] 所示，可将 3r 次乘法改进为 (2r + 2) 次，接近最优比例。使用这种方法，我们获得了可接受的模乘法速度：对于 95 位的 (n)，一次模乘法大约需要 0.003 秒，由于 (r = 11) 意味着每行有 (\lfloor 128 / (11 + 1) \rfloor = 10) 个块，因此可以同时执行 1280 次这样的乘法。

块式运算的另一个优点是所有相关值可以保存在寄存器中，从而避免了代价高昂的内存读取和存储操作。详细描述可参考 [2]。

graph TD;
    A[开始] --> B[选择b和r];
    B --> C[划分PE为块];
    C --> D[执行加法];
    D --> E{是否有进位};
    E -- 是 --> D;
    E -- 否 --> F[完成加法模n];
    F --> G[执行乘法];
    G --> H[使用Montgomery表示法];
    H --> I[计算结果];
    I --> J[结束];

5. 第二个实现及结果

将上一节的模运算融入我们的第一个 SIMD 椭圆曲线实现中，得到了程序的第二个版本。该版本每个块运行一条曲线，而非每个处理单元运行一条曲线，这意味着试验次数取决于 (n) 的大小。对于 95 位的 (n)，有 1280 次试验；对于 80 位的 (n)，(r) 变为 9，试验次数增加到 1536 次。

为降低计算曲线上点 (x) 的 (k) 次幂所需的群运算次数，我们采用了以下基本方法。根据第 3 节中 (k) 的定义，(k = \prod_{i = 1}^{l} q_i)，其中 ({q_1, q_2, \cdots, q_l}) 是小于 (B_1) 的素数幂集合。通常计算 (x^k) 的方法是先将 (x) 提升到 (q_1) 次幂，然后将结果提升到 (q_2) 次幂，依此类推，直到处理完所有 (q_i)。

对于某个整数 (m)，定义其权重 (w(m)) 为 (m) 的二进制表示中 1 的个数。若使用普通的重复平方和乘法进行幂运算，计算 (x^k \in G) 的代价为 (\sum_{i = 1}^{l} \lceil \log_2 q_i \rceil) 次群 (G) 中的平方运算和 (\sum_{i = 1}^{l} (w(q_i) - 1)) 次群 (G) 中的乘法运算。

设 (S = S_1 \cup S_2 \cup \cdots \cup S_s) 是 ({1, 2, \cdots, l}) 的一个划分，且 (Q_j = \prod_{i \in S_j} q_i)。显然，也可以通过先将 (x) 提升到 (Q_1) 次幂，然后将结果提升到 (Q_2) 次幂，依此类推，直到 (Q_s)，来计算 (x^k)。这种计算所需的群 (G) 中的平方运算次数与上述方法大致相同，但通过选择一个使得 (\sum_{j = 1}^{s} w(Q_j)) 较小的划分 (S)，可以显著减少群 (G) 中的乘法运算次数。然而，找到关于此度量的最佳划分通常是一个难题，在实践中，我们只能在合理的时间内找到可行的划分。

使用简单的贪心算法（(B_1 = 100000)），我们找到了一个子集基数至多为 2 的划分，其权重约为原始权重的一半。Bill Cook 在相同限制下得到了一个最优划分，但结果权重并没有显著降低。接着，我们考虑子集基数至多为 3 的情况，这使得权重约为原始权重的三分之一。鉴于寻找这些三元组所花费的时间，这可能是我们所能达到的最佳结果。

这些三元组是通过在 16K MasPar 上使用贪心算法找到的。为使运行时间可接受，我们按区间处理素数。具体而言，对于每个 (i \in {1, \cdots, 20})，我们确定区间 ([(i - 1) \cdot 10^5, i \cdot 10^5]) 内素数的小于 (2 \cdot 10^6) 的素数幂的三元组划分。这种按区间划分的方法也有助于选择界限。使用此技术，程序的总运行时间减少了 18%。

例如，素数 1028107、1030639 和 1097101 的权重分别为 10、16 和 11，但它们的乘积的权重为 8。

对于 (r = 11) 的 (n)，完成 1280 次椭圆曲线试验（(B_1 = 10^6)，(B_2 = 2 \cdot 10^7)）大约需要 34 小时。对于 (r = 12) 的 (n)，完成 1152 次试验（相同界限）大约需要 (34 \cdot 12 / 11) 小时，其他情况的运行时间也可类似推导。因此，这个第二个椭圆曲线程序允许在搜索 30 位因子时更平衡地选择参数。

我们使用该程序对多个数进行了分解，包括一些来自特定列表的合数和 RSA 数据安全分解挑战的“分区列表”中的数。截至目前，我们找到的最大因子有 40 位，这创造了新的椭圆曲线分解记录：
(p(11279) = 2658418735626949973617503123207968956766268614820186399554424770378507734924917342278622201969372653526213641483293)

被分解的数是最后两个因子的 89 位乘积，(p(11279)) 表示第 11279 个分区数。这个因子是在 1408 次试验中的一次，在第二阶段 (q = 1208269) 时找到的，这意味着我们相当幸运，因为此时找到该因子的概率仅为 0.7%。

此外，我们的 SIMD 椭圆曲线程序还有一个意外应用。在相关研究中，引入了一类新的伪随机生成器，这类生成器基于整数 (b)、(r) 和 (s)，使得 (m = b^r \pm b^s \pm 1) 为素数，其周期等于 (b) 模 (m) 的阶。确定 (b) 的阶需要对 (m - 1) 进行分解，这通常是困难的。

我们发现，例如 (b = 2^{53} - 1052) 在 511 位素数 (m = b^{31} + b^{16} - 1) 下的阶为 ((m - 1) / 2)，其中 (m - 1) 的因子 (b^8 + 1) 的分解是使用我们的 SIMD 椭圆曲线实现完成的：
(b^8 + 1 = 257 \cdot 16673 \cdot 275422758002613571762817294646047967240555733940016572674871759755285633142985718355214139747587789390113133069907813629489)

这是少数几个前端全局求逆产生合数因子的例子之一：一个块找到了 24 位因子，同时另一个块找到了 26 位因子。这个 (b) 和 (m) 为仅使用加法生成双精度浮点伪随机数提供了一种有吸引力的方法。

数据位数	试验次数	(B_1)	(B_2)	运行时间
95 位	1280	(10^6)	(2 \cdot 10^7)	约 34 小时
80 位	1536	-	-	-
(r = 12) 的 (n)	1152	(10^6)	(2 \cdot 10^7)	约 (34 \cdot 12 / 11) 小时

graph TD;
    A[计算x的k次幂] --> B[划分素数幂集合];
    B --> C[选择子集基数];
    C --> D{基数为2?};
    D -- 是 --> E[贪心算法找划分];
    D -- 否 --> F{基数为3?};
    F -- 是 --> G[贪心算法找三元组];
    F -- 否 --> H[其他情况];
    E --> I[计算权重];
    G --> I;
    I --> J[优化运算次数];
    J --> K[完成计算];

综上所述，我们通过不断优化椭圆曲线分解算法在 SIMD 机器上的实现，显著提高了分解性能，并在实际应用中取得了良好的成果。未来，随着硬件技术的不断发展，该算法有望在更大规模的分解任务中发挥更重要的作用。同时，其在伪随机数生成等领域的意外应用也为相关研究提供了新的思路。