大整数分解算法与实践

2023-08-02

大整数分解

1977 年，基于大整数分解难题的密码体系 RSA 诞生，命名自 Ron Rivest、Adi Shamir 和 Leonard Adleman 三位提出者姓氏的首字母。RSA 算法在公钥密码学上作出重要贡献并对全世界产生重要影响，三位提出者于 2002 年获该年度图灵奖。整数分解问题是当今几乎所有公钥密码安全性领域中最重要的课题之一。

RSA算法的安全性依赖于大整数分解难题，即

将一个大整数分解成若干大素因子的乘积。

有意思的是，虽然这是实践证明的一个「公认」难题，但现在仍无法证明其为 $P$问题 还是 $NP$问题。

针对小因子的几种典型的大整数分解算法有：

• Pollard's rho 算法

• Pollard’s p-1 算法

• ECM算法（Lenstra elliptic curve factorization/elliptic-curve factorization method )

此次我们将介绍 Pollard’s p-1 算法、ECM 算法与 ECM 算法的一个具体实践 —— GMP-ECM。

Pollard’s p-1 算法

Powersmooth：如果一个数 $m$ 的所有的质数幂$p^v$都是B-smooth的，即$p^v \le B$，则称$m!$为B-smooth。

Pollard’s p-1 算法是 John Pollard 在 1974 年设计的一种利用费马小定理的整数分解算法。根据费马小定理，$a^{K(p-1)} \equiv 1 \mod p$，假设待分解的 $n$ 存在一个因子为 $p$，我们可以表示为 $p|n$。对于 $p$ 存在$d = a^{p-1} - 1 \equiv 0 \mod p$，因此 $d$ 是 $p$ 的倍数，所以可以得到

$p = gcd(d, n) = gcd(a^{p-1} - 1, n)$

对于数据的选取，Pollard p-1 算法采用随机数的策略。如果随机不出结果，那么就从小于 $B$（B-powersmooth）的质数列表里去生成 $d$，直到质数列表被用尽。

步骤

输入：$n$（合数）

输出：$n$ 或失败的非平凡因子

1. 选择一个平滑边界 $B$
2. 定义：

   $M = \prod_{q \le B} {q^{\left \lfloor log_q{B} \right \rfloor }}$，其中 $q$ 是素数
3. 随机选择一个和 $n$ 互质的数 $a$
4. 计算：

   $g = gcd(a^M - 1, n)$
5. • 如果 $1 < g < n$，返回 $g$
   • 如果 $g=1$，选择一个更大的 $B$ 并转到步骤 2 或返回失败
   • 如果 $g=n$，选择一个更小的 $B$ 并转到步骤 2 或返回失败

ECM 算法

Pollard’s p-1 方法要求 $n$ 存在因子 $p$，并且 $p-1$ 拥有较小的因子，因此可以利用 $p=gcd(a^{p-1} - 1, n)$ 找到 $n$ 的因子。但是当 $p-1$ 的因子较大时，Pollard’s p-1 大概率不能得到 $n$ 的分解。

对于这个方法的改进，ECM 方法在 Pollard’s p-1 的基础上引入了随机的椭圆曲线，将原本的乘法群转换为椭圆曲线的点构成的群。椭圆曲线方法是一种基于数论的算法，利用椭圆曲线上的点的性质和运算规则，寻找满足一定条件的因子。

算法概述

假设 $p$ 是 $n$ 的一个质因数，考虑椭圆曲线 $E_{a,b} \mod p$，基于 $E_{a,b} \mod p$，根据 Hasse 定理，可以知道椭圆曲线的阶 $g$ 满足 $|g-(p+1)|<2 \sqrt{p}$。当椭圆曲线的 $a$ 和 $b$ 发生变化时，椭圆曲线的阶会落在范围 $[ p + 1 - 2\sqrt{p}, p + 1+ 2 \sqrt{p}]$ 中。

利用 ECM 算法寻找 $n$ 的因子分为三步：

1. 首先选择 $B$，并计算 $m = lcm(1,2,\dots ,B)$
2. 在 $Z/NZ$ 上选择一条随机的曲线与曲线上的随机一点 $P$
3. 计算 $mP$。在计算过程中，如果不能计算出对应的椭圆曲线中的点的加法，则找到了 $n$ 的因子，否则可以更换一条曲线进行计算，或返回失败

ECM 算法能够找到 $n$ 因子的原理在于，当随机选取的椭圆曲线阶 $g$ 为 B-smooth 时，在计算 $mP$ 的过程中，如果遇到 $gcd(v,p) = p$ 或者 $gcd(v,q) = q$，会计算得到 $gcd(v,n)$ 的非平凡因子。

一种 ECM 实现 —— GMP-ECM

GMP-ECM 的实现基于原有的 ECM 方法进行了优化。

首先是对整体的算法过程进行的优化。在 GMP-ECM 中计算因子分为两个步骤，每一个步骤选取两个$（B1，B2）$。

算法分为两个阶段，分别使用 $B1$ 与 $B2$ 进行运算。在 Stage 1 的计算与 ECM 类似。

Stage 1

首先计算在椭圆曲线 $E_{a,b}$ 上的点 $P$ 的积：

$$Q = \Pi_{\pi \le B_2} \pi^{\left \lfloor (log B_1)/(log \pi) \right \rfloor} P_0$$

在 Stage 1 的计算中，GMP-ECM 优化了椭圆曲线的选取与计算。

随机的椭圆曲线生成

在计算因子的过程中，当随机生成的椭圆曲线的阶在 $[ p + 1 - 2\sqrt{p}, p + 1+ 2 \sqrt{p}]$ 范围内变化，且随机的椭圆曲线是 B1-smooth 时，可以分解 $n$。在 GMP-ECM 中对于椭圆曲线的生成运用一定方法，可以使生成的椭圆曲线阶满足某些特征，例如：

• Suyama's form：保证椭圆曲线的阶 $g$ 能被 12 整除
• Montgomery's form：保证椭圆曲线的阶 $g$ 能被 3 整除

椭圆曲线运算

在 GMP-ECM 中使用 Montgomery's 形式的椭圆曲线，在求点的乘法或加法过程遵循如下规则：

• $P(x_P :: z_P), Q(x_Q :: z_Q) \rightarrow P + Q$
• 计算: $x_{P+Q} = 4z_{P-Q} * (x_P x_Q - z_P z_Q)^2$
• 计算: $z_{P+Q} = 4x_{P-Q} * (x_P z_Q - z_P x_Q)^2$
• $P(x_P :: z_P) \rightarrow 2P$
• 计算: $x_{2P} = (x_p^2 - z_P^2)^2$
• 计算: $z_{2P} = (4x_P z_P)[(x_p - z_P)^2 + d(4x_P z_P)]$, $d = (a+2) / 4$，其中 $a$ 为椭圆曲线 $b y^2 z = x^3 + ax^2z + xz^2$中的参数

Stage 2

GMP-ECM 的第一阶段主要计算椭圆曲线 $E$ 上的点 $Q$。如果失败，即 $gcd(n, zQ) = 1$，则在第二阶段将查找范围扩大至 $[B1，B2]$，并且新范围内存在一个质数 $\pi$，使得 $\pi Q = O_E \mod p$。

Stage 2 的主要思想为中间相遇（meet-in-the-middle）策略：

• 在 Stage 2 计算一个 $\sigma Q = ( x_\sigma : y_\sigma)$ 和 $\tau Q = ( x_\tau : y_\tau)$.
• 若 $\sigma Q + \tau Q = O_E \mod p$ 表明 $x_\sigma = x_\tau \mod p$。因此计算 $gcd(x_\sigma - x_\tau, n)$ 即可获得因子 $p$。

算法复杂度

ECM 方法找到一个数 $n$ 的因子 $p$ 预计需要的时间为 $O(L(p)^{\sqrt{2} + o(1)} M(\log n))$，这其中 $L(p) = e^{\sqrt{\log p \log \log p}}$。

$M(\log n)$ 代表了乘法模 $n$ 的复杂度。

在实际的运算过程中，由于算法存在随机性，找到数 $n$ 因子 $p$ 的时间可能存在差异。另外，算法能否分解与分解所需的时间也与 $B1$、$B2$ 的选取有关（可以参考 GMP-ECM 文档推荐的 $B1$、$B2$ 值进行计算）。

示例

以大整数 $N$（2048 bit）为例：

N = 0xA2862FB145AC7CE580E0BFF3CEFF42646050F8C611D36A3026E6EFB433B5CDD9EEFBA893E3F25C23A4951BA20992680162934CBDB6876CC791738C140C6EA6EC82938F7E18C5F0760367CF20883DCAB1D6885E222BBC7A1B5D434FD9FC148E387FA9AD69CE708FDDDE3E963F9AB2AF2046A37D7DBA21BC92E6F2A631C3E7770C77C95FAD6F1DC60D9F0645AEF2CC5D03D2151E35443FE0F90F1E2EAE58489C4450C8281EDCC58DDB409C306797C616954ACABCE4881E40ABDD2689A9F7596F29CD5CB42C752DA9306E7DB87CAC8957E3CA165421CF9E1B7220759A10588B386E33FD8E762E92C9F79D50150EF5FCA5F411DE23B8DFFE47A95D48ADDCF4797565

调用解析工具使用 GMP-ECM 方法分解：

ecm -c 2  25e4 1.3e8 -mpzmod threads: 2 mod: 1
A2862FB145AC7CE580E0BFF3CEFF42646050F8C611D36A3026E6EFB433B5CDD9EEFBA893E3F25C23A4951BA20992680162934CBDB6876CC791738C140C6EA6EC82938F7E18C5F0760367CF20883DCAB1D6885E222BBC7A1B5D434FD9FC148E387FA9AD69CE708FDDDE3E963F9AB2AF2046A37D7DBA21BC92E6F2A631C3E7770C77C95FAD6F1DC60D9F0645AEF2CC5D03D2151E35443FE0F90F1E2EAE58489C4450C8281EDCC58DDB409C306797C616954ACABCE4881E40ABDD2689A9F7596F29CD5CB42C752DA9306E7DB87CAC8957E3CA165421CF9E1B7220759A10588B386E33FD8E762E92C9F79D50150EF5FCA5F411DE23B8DFFE47A95D48ADDCF4797565
********** Factor found in step 2: 1021791499165844943393503 21 52761ms

成功找到素因子 1021791499165844943393503。