很早之前就了解利用Coppersmith来解决CTF中的问题，但是只是懂了皮毛，运用起来一直是一知半解，现在好好研究一下它的原理。

原理解释

度为$d$的首一多项式：$F(x)=x^d+a_{d-1}x^{d-1}+…+a_1x^1+a_0(\bmod M)$

我们需要寻找$x_0$，使得$F(x_0)\equiv 0(\bmod M)$，这里的$|x_0|<M^{\frac 1d}$。

构造$G(x_0)=0$来简化$F(x_0)$（这里取的是等号并非同余），转换为等式之后就可以用求根公式或者牛顿迭代法等方法进行快速求根。而$G(x)$的构造往往与格的LLL算法有关。

已知$X$为$|x_0|$的取值上界，可以用行向量$b_F=(a_0,a_1X,…,a_dX^d)$表示$F(x)$。

原理举例

取模数$M=323=17×19$，多项式$F(x)\equiv x^2+33x+215(\bmod M)$。

构造$G(x)=9F(x)-M(x+6)=9x^2-26x+3$，得到$G(3)=0$，也就是$F(3)\equiv 0(\bmod M)$，$x_0=3$。

Howgrave-Graham定理

如果$|x_0|＜X$，$F(x_0)\equiv 0(\bmod M)$，$|b_F|<\frac M{\sqrt{d+1}}$，那么$F(x_0)=0$。

证明：
$$
|F(x_0)|=|\sum^d_{i=0}a_ix_0^i|≤\sum^d_{i=0}|a_i||x_0|^i≤\sum^d_{i=0}|a_i|X^i≤\sqrt{d+1}|b_F|＜\frac{M\sqrt{d+1}}{\sqrt{d+1}}=M
$$
其中，第三个不等式用到了柯西不等式：
$$
(\sum^n_{i=1}x_iy_i)^2≤(\sum^n_{i=1}x_i^2)(\sum^n_{i=1}y_i^2)
$$
令$x_i≥0,y_i=1$，有：
$$
\sum^n_{i=1}x_i≤\sqrt{n(\sum^n_{i=1}x_i^2)}=\sqrt{d+1}|b_F|
$$
这样就能得到$-M＜F(x_0)＜M$，又有$F(x_0)\equiv 0(\bmod M)$，所以$F(x_0)=0$。

寻找$G(x)$

首先要把$F(x)=\sum^d_{i=0}a_ix^i$这个多项式变为首一多项式，如果$F(x)$不是首一多项式但是$gcd(a_d,M)=1$，那么可以对式子同时乘$a_d^{-1}$；如果$gcd(a_d,M)＞1$，那么可以把$M$拆开变为两个式子进行计算。

对于$d+1$个多项式$G_i(x)=Mx^i(0\leq i＜d)$，根据$F(x)$构造$d+1$维方阵$B$：
$$
B=\begin{pmatrix}M & 0 & … &0 &0 \0 & MX & … &0 &0 \… & & … & &… \0 & 0 & … &MX^{d-1} &0 \a_0 & a_1X & … &a_{d-1}X^{d-1} &X^d \\end{pmatrix}
$$
这里有$det(B)=M^dX^{\frac {d(d+1)}2}$。

LLL格基规约后得到的首行行向量$b’$根据LLL结论有$|b’|\leq 2^{\frac{n-1}4}det(B)^{\frac 1n}=2^{\frac d4}M^{\frac d{d+1}}X^{\frac d2}$。

根据上面的Howgrave-Graham定理，知道$2^{\frac d4}M^{\frac d{d+1}}X^{\frac d2}＜\frac M{\sqrt{d+1}}$，可以根据这个估计$X$上界的估计值。

例如我们取$M=10001$，$F(x)=x^3+10x^2+5000x-222$，$X=10$，求根，可以构造：
$$
B=\begin{pmatrix}M & 0 & 0 &0 \0 & MX &0 &0 \0 & 0 &MX^{2} &0 \-222 & 5000 X &10X^{2} &X^3 \\end{pmatrix}
$$
LLL后得到$b’=(444,10,-2000,-2000)$，构造$G(x)=444+x-20x^2-2x^3$，牛顿迭代得到$x_0=4$。

Coppersmith方法

通过上文的推导得到了$2^{\frac d4}M^{\frac d{d+1}}X^{\frac d2}＜\frac M{\sqrt{d+1}}$，Coppersmith方法通过改变了矩阵的构造对该结果进行了优化，这个式子可以化为等价形式：
$$
2^{\frac {n-1}4}M^{\frac d{n}}X^{\frac {d(d+1)}{2n}}＜\frac M{\sqrt{n}}
$$
所以可以通过增加格的维度$n$和模数$M$来扩大上界$X$，以便找到满足条件的解：

方法一可以通过插入多项式$xF(x),x^2F(x),…,x^kF(x)$，称作$x-shift$。

方法二可以利用$F(x)$的幂增加模数，比如$F(x)^k\equiv 0(\bmod M^k)$。

一元Coppersmith

定理：设$0＜\epsilon＜min{0.18,d}$，如果$|x_0|<\frac 12M^{\frac 1d-\epsilon}$，那么可以在多项式时间内找到小根。

一个$e$阶的多项式，给定$\beta$，可以快速求得模$n$或者模$n$的因数$p$下多项式在$n^\frac {β^2}e$以内的小根，这里有$p≥n^β(0＜\beta≤1)$。

在实际运用时通常使用sage的small_roots()函数，根据上面所引入的一系列理论，大致进行一下Coppersmith攻击：

生成一个以$x$为符号的一元多项式环，并将其赋值给名为$PR$的变量，这样可以直接使用多项式环中定义的未知数$x$来创建多项式：PR.<x> = PolynomialRing(Zmod(n))

定义求解的函数f，就是构造一个在模$n$或模$n$的因数为$0$条件下的方程。

在已知$p$高位的攻击中：f = x + p_high

规定$p,q≥n^β$，可知$β≠0.5$，应该稍小一些，$p,q$二进制位数相同时一般取$0.4$即可。

small_roots的$β$参数实际计算只用到的是第一位的小数，如果为两位也只使用一位小数的求解范围，所以一般形如p = getPrime(bits),q = geyPrime(bits)的$p$高位攻击通常只取$0.4$。

x0 = f.small_roots(X=2^kbits, beta=0.4)，这里内部使用的算法是LLL算法的变体。

这里引入一个问题，实际问题中我们需要知道$p$的多少位数才能进行攻击呢？

这里我做了一个实验，测试了一下究竟需要多少位已知高位才可以进行攻击：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
from Crypto.Util.number import getPrime

for bits in range(512,1025):
    p = getPrime(bits)
    q = getPrime(bits)
    n = p * q
    for i in range(bits//2,0,-1):
        test_bits = i
        p_h = p >> test_bits
        R.<x> = PolynomialRing(Zmod(n))
        f = (p_h << test_bits) + x
        x0=f.small_roots(2**test_bits,0.4)
        if x0:
            print(i,n.bit_length(),float(i/bits))
            break

实际得到的结果可以发现，$\frac {p_{unknown}}{p_{bits}}$的数值大约在$0.442$到$0.445$之间，通过这个规律，可以大致推断出大概需要多少位的已知位数。

f.small_roots()还可以通过修改$\epsilon$调整精度改变根的上限$X$的影响，默认值为$0.05$，通过测试如果调为$0.01$那么$\frac {p_{unknown}}{p_{bits}}$的数值大约在$0.482$到$0.485$之间，改为$0.1$那么范围会调整为$0.425$到$0.429$，随着数值的减小，我们所需的已知位数也在减小，但是计算时间也会逐渐增加。

二元Coppersmith

$Z^2$表示整数的二维笛卡尔坐标集合，包括所有由整数组成的有序对。

我们把含有两元的多项式$F(x,y)\equiv 0(\bmod M)$的度记为$d$（$x^iy^j,i+j≤d$）,同样有$x,y$的上界$X,Y$，如果$XY＜M^{\frac 1d-\epsilon}$，利用一元引入$G(x)$的思想，可以找到并计算出多项式$F_1(x,y),F_2(x,y)$，使得所有$(x_0,y_0)\in Z^2,|x_0|＜X,|y_0|＜Y$而且$F(x_0,y_0)\equiv 0(\bmod M)$，使得$F_1(x_0,y_0)=F_2(x_0,y_0)=0$。

定理：我们取$d$使其满足$deg_x(F(x,y)),deg_y(F(x,y))≤d$，定义：
$$
F(x,y)=\sum_{0≤i,j≤d}F_{i,j}x^iy^j
$$
取$W$使其满足：
$$
W=max_{0≤i,j≤d}{|F_{i,j}|X^iY^j}
$$
如果$XY＜W^{\frac 2{3d}}$，那么就存在一种算法能在多项式时间里求出所有的$(x_0,y_0)\in Z^2$，满足$F(x_0,y_0)=0$，此时有$|x_0|≤X,|y_0|≤Y$。

简单讲一下Coron对这个定理的解释：

用到了多项式移位构造矩阵的思想，选择了一个较大的$k\in N$满足$0\leq i,j＜ d+k$并考虑$k^2$多项式：
$$
s_{a,b}(x,y)=x^ay^bF(x,y)\ \ \mathrm{for}\ \ 0\leq a,b＜k
$$
紧接着构造一个$k^2+(d+k)^2$行$(d+k)^2$列的矩阵。

例如这里取$F(x,y)=axy+bx+cy+d=127xy-1207x-1461y+21$，已知$X=30,Y=20$，使得模数$M=127^4$。

我们取$k=2$，所以可以构造一个$13×9$的矩阵：

构造这个矩阵由上下两部分组成，上部分根据$k=2$和$s_{a,b}(x,y)$可以得到$a,b\in{0,1}$，得到四个相应的多项式：
$$
\begin{cases}XY\cdot F(x,y)=aX^2Y^2+bX^2Y+cXY^2+dXY\X\cdot F(x,y)=aX^2Y+bX^2+cXY+dX\Y\cdot F(x,y)=aXY^2+bXY+cY^2+dY\F(x,y)=aXY+bX+cY+d\end{cases}
$$
根据这四个多项式可以构成一个$9×4$的矩阵，每一行根据多项式每一项的幂来填充，总共$9$种情况，为了能让最终得到的矩阵能规约出较小的向量，再构造一个$9×9$的对角矩阵拼接在多项式矩阵下面。拼接后的矩阵可以提供更多的线性相关性，使得规约过程中可以对向量进行更多的操作：
$$
A=(X^2Y^2,X^2Y,XY^2,XY,X^2,Y^2,X,Y,1)
$$

计算它的Hermite标准型矩阵（空位全部为$0$，星号为无关数值）：

取$L$矩阵由$B’$的第五到第九行、第五到第九列这个$5×5$上三角矩阵构成，对其进行LLL格基规约，可以得到第一行相应转化为多项式：
$$
(−16129X^2, −16129Y^2, 1048258X, 983742Y, −28446222)
$$

$$
G(x,y)=−16129x^2−16129y^2+ 1048258x+983742y−28446222
$$

计算得到$(x,y)=(21,21),(23,19)$。

这里如果我们取$k$为其他值，那么相应的矩阵也会不同，对应的移位多项式也会越多，格基规约的精度越高，但是计算速度也会变慢。例如取$k=3$，则会得到$16×25$矩阵，移位多项式总共有$7$个，转化后计算即可。

这里提供一个根据上面的推导的Coppersmith代码。理论上根据同样的思想可以求解$n$元的Coppersmith：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
import itertools

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

    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 []

使用时注意f是模$N$的同余式，bounds是未知量（可以是一个或者多个）的上界，m是一个参数使得$f(x)^m\equiv0(\bmod N^m)$，可以增加精度，d跟上文推导的移位用的$k$是一回事，实际使用时需要调试m和d的值。

[CISCN 2021华南]small

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
import random, hashlib
from Crypto.Util.number import getPrime
from secret import x, y, flag

BITS = 70

assert(2**BITS < x < 2**(BITS+1))
assert(2**BITS < y < 2**(BITS+1))

m = str(x) + str(y)
h = hashlib.sha256()
h.update(m.encode())
assert(flag == "flag{" + h.hexdigest() + "}")


p = getPrime(512)
a = getPrime(510)
b = getPrime(510)
c = (1 + a * x * y ** 2 + b * x ** 2 * y) % p
print("p = " + str(p))
print("a = " + str(a))
print("b = " + str(b))
print("c = " + str(c))

'''
p = 8813834626918693034209829623386418111935369643440896703895290043343199520112218432639643684400534953548489779045914955504743423765099014797611981422650409
a = 2817275225516767613658440250260394873529274896419346861054126128919212362519165468003171950475070788320195398302803745633617864408366174315471102773073469
b = 1763620527779958060718182646420541623477856799630691559360944374374235694750950917040727594731391703184965719358552775151767735359739899063298735788999711
c = 2298790980294663527827702586525963981886518365072523836572440106026473419042192180086308154346777239817235315513418426401278994450805667292449334757693881
'''

观察到核心式子：
$$
c\equiv 1+axy^2+bx^2y\pmod p
$$
这里的$x,y$都比较小，构造二元coppersmith，直接求解就能得到flag：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
from Crypto.Util.number import *
import itertools
import hashlib

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

    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 []

p = 8813834626918693034209829623386418111935369643440896703895290043343199520112218432639643684400534953548489779045914955504743423765099014797611981422650409
a = 2817275225516767613658440250260394873529274896419346861054126128919212362519165468003171950475070788320195398302803745633617864408366174315471102773073469
b = 1763620527779958060718182646420541623477856799630691559360944374374235694750950917040727594731391703184965719358552775151767735359739899063298735788999711
c = 2298790980294663527827702586525963981886518365072523836572440106026473419042192180086308154346777239817235315513418426401278994450805667292449334757693881

p.<x,y>=Zmod(p)[]
f=(1+a*x*y**2+b*x**2*y)-c
x,y=small_roots(f,bounds=(2^70,2^71))[0]
m = str(x) + str(y)
h = hashlib.sha256(m.encode()).hexdigest()
flag='NSSCTF{'+h+'}'
print(flag)