Coppersmith

2026-05-16/更新于 2026-07-15

Coppersmith 算法是由美国数学家 Don Coppersmith 于 1996 年提出的一类 基于格的多项式时间算法，核心用于求解 模大合数 的低次多项式小根与整数域上多元多项式小根问题。

理论推导

假设我们有一个首一的多项式 $F (x)$ ，其度数（最高次幂）为 $d$ 。我们想要解同余方程

F (x) \equiv 0 (\mod n)

如果模数 $n$ 可以被分解为 $p_{1}^{e_{1}}, p_{2}^{e_{2}}, \dots p_{k}^{e_{k}}$ ，那么根据模运算的性质，原方程可以转化为

\begin{matrix} F (x) \equiv 0 (\mod p_{1}^{e_{1}}) \\ F (x) \equiv 0 (\mod p_{2}^{e_{2}}) \\ ⋮ \\ F (x) \equiv 0 (\mod p_{k}^{e_{k}}) \end{matrix}

这样 $k$ 个方程的组合。对于每一个素数 $p_{i}$ ，我们先求解 $F (x) \equiv 0 (\mod p_{i})$ 。在计算机中，求解 $G F (p)$ 上的多项式方程是成熟的。等到解出这些基本的同余方程之后，使用 Hensel 引理就能得到原方程模 $p^{2}, p^{3}, \dots$ 的根。所以一旦 $n$ 被分解，求根就不再是难事。对于已知分解的模数 $n$ ，使用 sage 求根的脚本如下：

# sage
from itertools import product
def solve_mod(coeffs, factors):
    # coeffs 是按照降序排列的系数列表（最高次幂在前），如 x^2 - 1 是 [1, 0, -1]
    # factors 是 n 的质因数分解，如 [(p1, k1), (p2, k2), ...]
    # sage 底层传参要求传升序排列的列表，反转一下
    coeffs_rev = coeffs[::-1]
    # 存放每一个模为 pi^ki 的同余方程的解，如 [[1, 2], [1, 4]]
    all_roots = []
    # 存放 pi^ki
    mods = []
    
    for p, k in factors:
        mod = p**k
        mods.append(mod)
        # 构建子环及其多项式环
        R = Zmod(mod)
        P.<x> = PolynomialRing(R)
        # 传入系数列表，生成多项式
        f = P(coeffs_rev)
        # 调用求根函数。sage 返回的是元组(解，重数)，取 r [0] 即可
        roots = [r[0] for r in f.roots()]
        all_roots.append(roots)
    
    final_roots = []
    # product 函数用于求笛卡尔积，也即列表内所有参数的全部排列组合，组合的数目即 all_roots 内的列表数目
    for rems in product(*all_roots):
        # 类型转换
        rems_int = [Integer(r) for r in rems]
        # 用 crt 解最小的同余方程组
        combined_root = crt(rems_int, mods)
        final_roots.append(combined_root)
    return sorted(final_roots)

例如求解 $x^{2} - 1 \equiv 0 (\mod 15)$ ，传入 solve_mod([1, 0, -1], [(3, 1), (5, 1)])，就能得到解集 [1, 4, 11, 14]。

但如果 $n$ 难以分解，求根则进化为 NP 难的问题，此时引入 Coppersmith 方法。

对于在模 $n$ 下度为 $d$ 的首一多项式 $F (x) = x^{d} + a_{d - 1} x^{d - 1} + \dots + a_{1} x + a_{0}$ ，使用 Coppersmith 方法能在多项式时间内找到满足 $F (x_{0}) \equiv 0 (\mod n)$ 且 $| x_{0} | < n^{\frac{1}{d} - ϵ}$ 的小解 $x_{0}$ ，其中 $ϵ$ 随着格维度的增大而趋近于 $0$ ，在 sage 中可以显式调整（epsilon = ），默认为 $0.125$ 。如果初始多项式不是首一的，只需要调用一下 f = f.monic() 就好。

接下来讲一讲算法思路。先引入 Howgrave-Graham 定理：

设 $h (x) \in Z [x]$ 是一个最高次为 $d$ 的整系数多项式， $n$ 是一个正整数。假设存在一个未知的整数 $x_{0}$ 满足：

同余条件： $h (x_{0}) \equiv 0 (\mod n)$
小根边界： $| x_{0} | < X$ （即 $X$ 是 $x$ 的上界， $X$ 的数值推导后面会说）

为什么要特意提到上界 $X$ ？Howgrave-Graham 定理的核心思想就是把 $G F (n)$ 上的方程转化为整数环上的方程。如果一个多项式满足 $h (x_{0}) \equiv 0 (\mod n)$ ，且 $| h (x_{0}) | < n$ ，那么 $h (x_{0}) = 0$ 。这就实现了从模方程到普通方程的转化。但我们该怎么保证 $| h (x_{0}) | < n$ ？就只能用 $| h (X) |$ 来估计。

我们定义多项式 $h (X x)$ 。假设 $h (x) = \sum_{i = 0}^{d} a_{i} x^{i}$ ，那么 $h (X x) = \sum_{i = 0}^{d} a_{i} X^{i} x^{i}$ 。定义该函数的 欧几里得范数 为所有系数构成的向量的长度，即：

| | h (X x) | | = \sqrt{\sum_{i = 0}^{d} (a_{i} X^{i})^{2}}

从这里可以看出， $X$ 的参与实现了加权的 LLL 算法，起到了约束高次幂的作用。Howgrave-Graham 定理指出，如果长度 $| | h (X x) | |$ 满足：

| | h (X x) | | < \frac{n}{\sqrt{d + 1}}

其中 $d + 1$ 可以记为 $ω$ ，它的意义是多项式的项数。此时 $x_{0}$ 不仅是模方程的根，同时也是多项式在整数环上的根，即 $h (x_{0}) = 0$ 。下面证明一下。

设新系数为 $A_{i} = a_{i} X^{i}$ ，既然 $| x_{0} | < X$ ，那么 $| \frac{x_{0}}{X} | < 1$ ，有下面的等式：

| h (x_{0}) | = | \sum a_{i} x_{0}^{i} | = | \sum (a_{i} X^{i}) (\frac{x_{0}}{X})^{i} | = | \sum A_{i} (\frac{x_{0}}{X})^{i} |

因为 $| \frac{x_{0}}{X} | < 1$ ，所以

| h (x_{0}) | \leq \sum | A_{i} |

对任意向量 $u, v$ ，有 $| u \cdot v | \leq | | u | | \cdot | | v | |$ 。这个不等式被称作柯西 - 施瓦兹不等式。

对于 $d + 1$ 维向量 $A = (A_{0}, A_{1}, \dots, A_{d})$ ，取 $u = A, v = (1, 1, \dots, 1)$ ，就有：

\sum_{i = 0}^{d} | A_{i} | \leq \sqrt{d + 1} \cdot \sqrt{\sum_{i = 0}^{d} A_{i}^{2}}

代入后得到 $| h (x_{0}) | \leq \sqrt{d + 1} \cdot | | h (X x) | |$ 。只要找到的 $| | h (X x) | | < \frac{n}{\sqrt{d + 1}}$ ，再次代入，就绝对能保证 $| h (x_{0}) | < n$ 。在实际使用 Coppersmith 方法时，由于我们会构造模 $n^{m}$ 的多项式，Howgrave-Graham 定理的判定条件也会相应地放宽到 $| | h (X x) | | < \frac{n^{m}}{\sqrt{ω}}$ 。

Howgrave-Graham 定理告诉我们，只要能找到一个系数向量的范数 足够小 的多项式，我们就可以将困难的模方程问题转化为简单的整数方程问题。面对这样一个求范数最小值的问题，LLL 算法可以高效解决。但原始方程只有一个，LLL 无从下手。此时 Coppersmith 提出了构造辅助多项式的方法。

选择一个整数 $m$ （决定了格的维度， $m \approx \frac{1}{d \cdot ϵ}$ ），定义以下的辅助多项式集合：

g_{i, j} (x) = n^{m - j} \cdot x^{i} \cdot f (x)^{j}

其中 $0 \leq j < m$ 且 $0 \leq i < d$ 。这个辅助多项式的三部分都很重要。

对于 $f (x)^{j}$ 来说，它的作用是针对性地放大模数。如果 $f (x)$ 有根 $x_{0}$ ，即 $f (x_{0}) = k n$ 。那么两边同时幂次，就有 $f (x)^{j} = k^{j} n^{j}$ ，同余的模数从 $n$ 抬高到了 $n^{j}$ 。
$n^{m - j}$ 的作用是抹平差异，把所有构造出来的多项式拉平到模 $n^{m}$ 上。 $f (x_{0})^{j} \cdot n^{m - j} = k^{j} n^{m}$ ，现在所有的多项式都统一模 $n^{m}$ 为 $0$ 了。
$x^{i}$ 主要用来填充矩阵。如果不加这一项，仅靠 $f (x)^{j} \cdot n^{m - j}$ 只能构造出 $m$ 个多项式，而格基矩阵的维度需要是 $m \times d$ 。由于 $0 \leq i < d$ ，对于每一个 $j$ 阶段生成的多项式， $x^{i}$ 的参与（相当于将多项式系数整体向左平移 $i$ 位）把数量扩充到了原来的 $d$ 倍，这样就能凑齐 $m \times d$ 个多项式，保证格基矩阵是一个 满秩的下三角矩阵。从代数角度考虑的话，我们已知 $x_{0}$ 满足 $f (x_{0})^{j} \cdot n^{m - j} \equiv 0 (\mod n^{m})$ 。当我们构造新多项式，并代入 $x_{0}$ 时，由于 $x_{0}^{i}$ 是一个具体的整数，乘以 $0$ 依然得 $0$ 。因此， $x_{0}^{i} \cdot f (x_{0})^{j} \cdot n^{m - j} \equiv 0 (\mod n^{m})$ 始终成立。到此我们就证明了引入 $x^{i}$ 既能扩容格基矩阵，也不会改变原方程的根。

解析了辅助多项式之后，我们来看它形成的格基。假设目标多项式是 $f (x) = x^{2} + a x + b$ ，我们选择维度参数 $m = 2$ 。那么根据辅助多项式的定义，需要构造以下四个多项式：

$g_{0, 0} (x) = n^{2} \cdot f (x)^{0} = n^{2}$
$g_{1, 0} (x) = n^{2} \cdot x \cdot f (x)^{0} = n^{2} x$
$g_{0, 1} (x) = n^{1} \cdot f (x)^{1} = n (x^{2} + a x + b)$
$g_{1, 1} (x) = n^{1} \cdot x \cdot f (x)^{1} = n (x^{3} + a x^{2} + b x)$

观察到辅助多项式出现了高于 $x^{2}$ 的 $x^{3}$ 项，这说明辅助多项式的最高次幂与原始多项式无关，只与维度参数 $m$ 有关。如果选择 $m = 3$ ，最终会出现 $x^{5}$ 。但无论次数有多高，转化到整数环上求解都很简单，影响运行时间的是与格维度息息相关的 LLL 算法。接着为了满足 Howgrave-Graham 定理的范数要求，把 $x$ 替换为 $X x$ ：

$g_{0, 0} (X x) = n^{2}$
$g_{1, 0} (X x) = n^{2} X x$
$g_{0, 1} (X x) = n b + n a (X x) + n (X x)^{2}$
$g_{1, 1} (X x) = 0 + n b (X x) + n a (X x)^{2} + n (X x)^{3}$

最后提取各多项式的系数，行顺序就是上面的编号顺序，列按照 $x^{0}, x^{1}, x^{2}, x^{3}$ 排列，就得到了格基矩阵 $L$ ：

L = [\begin{matrix} n^{2} & 0 & 0 & 0 \\ 0 & n^{2} X & 0 & 0 \\ n b & n a X & n X^{2} & 0 \\ 0 & n b X & n a X^{2} & n X^{3} \end{matrix}]

这是一个下三角矩阵。事实上，无论格多么复杂，最终得到的一定是一个下三角矩阵，这是由辅助多项式的代数性质决定的。考虑普遍的情况，设格的维度为 $ω = m \cdot d$ 。下三角矩阵的行列式为

det (L) = n^{\frac{d m (m + 1)}{2}} X^{\frac{ω (ω - 1)}{2}}

其中 $n$ 是模数， $d$ 是最高幂次， $m$ 是格基的维度参数， $X$ 是辅助因子。 $ω = m d$ ，是格基矩阵的维度。在这样的格中，LLL 算法能在多项式时间内找到一个相对最短的系数向量，其对应着一个系数极小的多项式 $h (x)$ 。LLL 算法保证了这个最短向量的长度上界：

| | h (X x) | | \leq 2^{\frac{ω - 1}{4}} det (L)^{\frac{1}{ω}}

重新回到 Howgrave-Graham 定理。由于我们之前用辅助多项式把模数提升到了 $n^{m}$ ，所以有

| | h (X x) | | < n^{m} \cdot \frac{1}{\sqrt{ω}}

忽略 $2^{\frac{ω - 1}{4}}$ 和 $\frac{1}{\sqrt{ω}}$ 两个常数项，代入后有关系式

det (L)^{\frac{1}{ω}} < n^{m}

代入之前行列式的值：

{(n^{\frac{d m^{2}}{2}} \cdot X^{\frac{m^{2} d^{2}}{2}})}^{\frac{1}{m d}} < n^{m}

初步化简：

n^{\frac{m}{2}} \cdot X^{\frac{m d}{2}} < n^{m}

两边同除以 $n^{\frac{m}{2}}$ ：

X^{\frac{m d}{2}} < n^{\frac{m}{2}}

开方，得到最终小根 $x_{0}$ 的理论上界 $X$ ：

X < n^{\frac{1}{d}}

但是这只是理论上界。实际上，如果想尽量逼近这个上界，就需要原始关系式 $X < n^{\frac{1}{d} - ϵ}$ 中的 $ϵ$ 尽量小。一旦 $ϵ$ 开始小，就直接影响到矩阵维度 $ω = \frac{1}{ϵ}$ ，耗时也会大幅增加。

实际应用

m 泄露

假设题目给出了 $n, e, m_{p a r t}$ ，结合 flag 的头尾（可选），在小公钥指数（ $e = 3$ ）的情况下，可以直接还原出 $m_{f u l l}$ 。这里举 $m$ 高位泄露的情况：

c \equiv m^{e} \equiv (m_{h i g h} + x)^{e} (\mod n)

构造 $f (x) = (m_{h i g h} + x)^{e} - c$ ，则有 $f (x) \equiv 0 (\mod n)$ 。此时的小根上界 $X = n^{1 / e}$ 。也就是说只要 $m$ 的未知位数小于 $\frac{1}{3} nbits$ ，就能还原完整的 $m$ 。

from Crypto.Util.number import *

n = 
e = 3
c = 
m_high = 
k = # 低位丢失了多少位信息

R.<x> = PolynomialRing(Zmod(n))
f = (m_high + x)^e - c
roots = f.small_roots(X = 2^k)

m = m_high + int(roots[0])
print(long_to_bytes(m).decode())

值得注意的是，虽然有理论上界，但真正在跑内置的 small_roots 方法时，一定要填实际情况下小根的上界（比如截断丢失了 300 位， $m_{l o w}$ （也就是 $x$ ）的上界就是 $2^{300}$ ）

p 泄露

对于 $p$ 高位泄露的问题，假设已知的高位是 $p_{h i g h}$ （低位已全部补 0），低位未知。那么就有

p = p_{h i g h} + x_{0}

构造一个简单的多项式 $f (x) = p_{h i g h} + x$ ，当代入未知低位 $x_{0}$ 时，有：

f (x_{0}) = p_{h i g h} + x_{0} = p \equiv 0 (\mod p)

同余式里是 $(\mod p)$ ，而不是 $(\mod n)$ 。我们不知道 $p$ ，就需要调用 sage 中带有 beta 参数的 small_roots，它可以解决这种因子未知而其倍数（ $n = p \cdot q$ ）已知的情况。模方程 beta 参数决定了 $p$ 相对于 $n$ 的 大小下界，即

p \geq n^{β}

如果不显式声明的话， $β$ 默认为 $1$ ，也即求解的模方程中模数就是 $n$ 本身。但在 $p$ 高位已知的情况下， $p \approx \sqrt{n} = n^{0.5}$ 。事实上，如果设定 beta = 0.5，当 $p, q$ 的比特数相同时，会有 50%的概率找不出 p。最优的数值是 beta = 0.499，一定能保证找到 $p$ ，且能保证小根上界 $X$ 不过度坍缩。下面是解释：

已知：

1 2	`p = getPrime(1024) q = getPrime(1024)`

由 $β$ 的定义， $β = \frac{\log_{2} (p)}{\log_{2} (n)}$ 。 $l o g_{2} (n) = \log_{2} (p) + \log_{2} (q)$ ，有

β = \frac{\log_{2} (p)}{\log_{2} (p) + \log_{2} (q)} = \frac{1}{1 + \frac{\log_{2} (q)}{\log_{2} (p)}}

假设最坏的情况，即 $β$ 最小， $\frac{\log_{2} (q)}{\log_{2} (p)}$ 最大。根据 getPrime() 函数的范围， $q$ 的上界趋近于 $2^{1024}$ ， $p$ 的下界趋近于 $2^{1023}$ 。此时得到的 $β$ 最小为 $\frac{1023}{2047} = 0.4997$ 。设定的 beta = 0.4990 比下界更小。同理，对于 RSA-1024 和 RSA-4096，它们的 $β$ 下界是 $0.4995$ 和 $0.4999$ 。因此，设定 beta = 0.499 是完全合适的。相比于 beta = 0.5，只会让小根上界从 $n^{0.25}$ 下降到 $n^{0.249}$ ，但却保证了一定能找到 $p$ 。

根据 Coppersmith 在未知因子上寻找小根的定理，实际能够求出的根的上限 $X$ 由以下公式确定：

X = n^{\frac{β^{2}}{d} - ϵ}

代入 $d = 1, β = 0.5$ ：

X = n^{0.25 - ϵ}

也就是未知低位的真实大小 $x_{0}$ 必须小于 $n^{0.25 - ϵ}$ 。转换为比特数，就有：

k_{u n k n o w n} = \frac{L}{4} - ϵ \cdot L

对于 RSA-1024，能够求解 $p$ 的未知位数最多为 $L / 4$ ，即 256 位。现实情况中对于不同的 $ϵ$ ，未知位数不能超过以下数值：

$ϵ = 0.05$ ， $k_{u n k n o w n} = 204 bits$ ， $k_{k n o w n} = 308 bits$
$ϵ = 0.02$ ， $k_{u n k n o w n} = 235 bits$ ， $k_{k n o w n} = 277 bits$
$ϵ = 0.01$ ， $k_{u n k n o w n} = 245 bits$ ， $k_{k n o w n} = 267 bits$

随着 $ϵ$ 的下降，能求解的未知位数会变多，但求解花费的时间也会大幅增长。

需要注意的是，下面的脚本需要 $p, q$ 位数相同时才能使用，即标准 RSA。检测脚本：

beta = 
epsilon = 
nbits = 
d = 1

maxRootBits = nbits * ((beta * beta) / d - epsilon)

print(f"小根上界为{maxRootBits:.0f}bits")

一旦满足条件，coppersmith 就能发挥作用。高位泄露常见形式是 p >> k 或 p >> k << k，还可以写成掩码形式 p & (((1 << k) - 1) << (pbits - k))。脚本如下：

from Crypto.Util.number import *

n = 
e = 
c = 
p_high = 
k = # 低位丢失了多少位信息

pbits = isqrt(n).nbits()
# 检查是 p >> k 还是 p >> k << k
if pbits != p_high.nbits():
    p_high = p_high << k
R.<x> = PolynomialRing(Zmod(n))
f = p_high + x
roots = f.small_roots(X = 2^k, beta = 0.499)

for r in roots:
    # 提取模环元素要加 int()类型转换
    p = p_high + int(r)
    if n % p == 0:
        q = n // p
        phi = (p - 1) * (q - 1)
        d = inverse_mod(e, phi)
        m = pow(c, d, n)
        print(long_to_bytes(m).decode())

低位泄露常见形式是 p % (2**k)，或 p & ((1 << k) - 1)。脚本：

from Crypto.Util.number import *

n = 
e = 
c = 
p_low = 
k = # 低位保留了多少位信息，不用 p_low.nbits()的原因是避免丢失前导 0

R.<x> = PolynomialRing(Zmod(n))
f = 2^k * x + p_low
f = f.monic()
pbits = isqrt(n).nbits()
roots = f.small_roots(X = 2^(pbits - k), beta = 0.499)

for r in roots:
    p = p_low + int(r) * 2^k
    if n % p == 0:
        q = n // p
        phi = (p - 1) * (q - 1)
        d = inverse_mod(e, phi)
        m = pow(c, d, n)
        print(long_to_bytes(m).decode())

d 泄露

相比于之前的 $m, p$ 泄露，通过 $d$ 泄露实施攻击之所以能够成立，是因为 RSA 的公私钥之间存在极强的代数约束。下面先进行数学推导：

d 低位泄露

根据 $e d \equiv 1 (\mod ϕ (n))$ ，简单变形后有 $e d = k (n - p - q + 1) + 1$ 。在两边同时对 $2^{t}$ 取模：

e d_{0} \equiv k (n_{0} - p_{0} - q_{0} + 1) + 1 (\mod 2^{t})

这里的 $d_{0}$ 就是 $d$ 的低位， $t$ 是泄露的位数。代入 $q_{0} = n_{0} / p_{0}$ ，并在方程两边同乘 $p_{0}$ ，稍加整理：

k p_{0}^{2} + (e d_{0} - k n_{0} - k - 1) p_{0} + k n_{0} \equiv 0 (\mod 2^{t})

在之前 RSA 的文章里，我们已经证明了 $k \in (0, e)$ 。因此只需要枚举 $k$ ，对于每一个猜想的 $k$ ，解模方程（对于有明显分解的模数 $2^{t}$ ，可以用 Hensel 提升法从模 $2$ 一直提升到模 $2^{t}$ ）。解出的 $p_{0}$ 就是 $p$ 的低 $t$ 位，接下来只需要按之前的方法恢复 $p$ 就可以了。脚本：

from Crypto.Util.number import *

n = 
e = 
c = 
d_low = 
t = # 低位保留了多少位信息

n_low = n % (2^t)
R.<x> = PolynomialRing(Zmod(2^t))
candidates = []
for k in range(1, e):
    f = k*x^2 + (e*d_low - k*n_low - k - 1)*x + k*n_low
    # 这里虽然调用了普通的.roots()，但 sage 不能在复合数模环上求重数，需要显式把这个功能关闭
    roots = f.roots(multiplicities=False)
    for r in roots:
        # 这里不再是元组了
        p_low = int(r)
        candidates.append(p_low)

# 退化回 p 的低位泄露
R_prime.<x> = PolynomialRing(Zmod(n))
for p_low in candidates:
    f = 2^t * x + p_low
    f = f.monic()
    pbits = isqrt(n).nbits()
    roots = f.small_roots(X = 2^(pbits - t), beta = 0.499)
    for r in roots:
    	p = p_low + int(r) * 2^t
    	if n % p == 0:
        	q = n // p
        	phi = (p - 1) * (q - 1)
        	d = inverse_mod(e, phi)
        	m = pow(c, d, n)
        	print(long_to_bytes(m).decode())

d 高位泄露

二元 copper

Sagemath 没有内置的二元 coppersmith 实现，只能手搓。这里先贴一下板子：

import itertools

def small_roots(f, bounds, m=1, d=None):
	if not d:
		d = f.degree()

	if isinstance(f, Polynomial):
		x, = polygens(f.base_ring(), f.variable_name(), 1)
		f = f(x)

	R = f.base_ring()
	N = R.cardinality()
	
	f /= f.coefficients().pop(0)
	f = f.change_ring(ZZ)

	G = Sequence([], f.parent())
	for i in range(m+1):
		base = N^(m-i) * f^i
		for shifts in itertools.product(range(d), repeat=f.nvariables()):
			g = base * prod(map(power, f.variables(), shifts))
			G.append(g)

	B, monomials = G.coefficient_matrix()
	monomials = vector(monomials)

	factors = [monomial(*bounds) for monomial in monomials]
	for i, factor in enumerate(factors):
		B.rescale_col(i, factor)

	B = B.dense_matrix().LLL()

	B = B.change_ring(QQ)
	for i, factor in enumerate(factors):
		B.rescale_col(i, 1/factor)

	H = Sequence([], f.parent().change_ring(QQ))
	for h in filter(None, B*monomials):
		H.append(h)
		I = H.ideal()
		if I.dimension() == -1:
			H.pop()
		elif I.dimension() == 0:
			roots = []
			for root in I.variety(ring=ZZ):
				root = tuple(R(root[var]) for var in f.variables())
				roots.append(root)
			return roots

	return []

其中 m 决定格基复杂度。