原文链接：https://tover.xyz/p/LLL-attack-equation/
上一次写格攻击的教程已经是两年前的Wiener's v.s Lattices，两年来新学了不少的东西，也纠正了以前不少的错误，所以决定新开一篇文章，系统性地讲一下格攻击之解小未知数方程。
另一个比较重要的原因是，今年在山石夏令营中刚好讲了相关的内容，由于版权问题当时的PPT我就不公开了，反正文章应该会比当时说的更详细。
想要看懂以下内容，你可能需要：
引言
首先格攻击准确地讲应该是格规约攻击，格规约是一个很强大的工具，它可以解决很多种问题，这里只说最常被用来解的一种问题，解小未知数方程
所以先看看什么是小未知数方程，顾名思义就是他的未知数是“小”的，但什么是“小”的呢，大小其实是需要比较的，准确来说对于具体问题需要算一个界，小于这个界的未知数才是小的
所以不妨先来看一个栗子
bits = 512
while True:
q = random_prime(2^bits, lbound=2^(bits - 1))
R = GF(q)
f = R(random_prime(2^(bits//2 - 1)))
g = R(random_prime(2^(bits//2 - 1), lbound=2^(bits//4 - 1)))
if gcd(f, q*g) == 1:
h = f^(-1) * g
break
print('q = %d' % q)
print('h = %d' % h)
# q = 11771107769919685639025680019132176714690799312096361547135769537209645333910623545026199750211710876592602157283472114417112848369348753900020624817689903
# h = 10244740934987488963223294332625298216569051945560244475215923413678863448556753932134017729009652572232483956393290621358806004252483727672873796530687963
# Aim: find f and g
这个例子，很多比赛都曾出现过和这个栗子类似的题目，其实他是教科书[HPS14]中1维NTRU的密钥生成部分
整理一下这个栗子，现在有条件f, g < \sqrt{q/2}、\text{gcd}(f, qg) = 1和h \equiv f^{-1} g \pmod q，已知(q, h)，求(f, g)
首先撇开大小关系不看，目前只有一个同余方程，但这个方程里面有两个未知数（其实还有一个隐藏的k，因为还需要展开同余式），根据我们中学所学的知识，一个方程两个未知数是不能解或者说是多解的，代入这个栗子就是不能找到生成密钥时所用的f和g
但是栗子中所给的大小关系改变了上述情况，虽然符合条件h \equiv f^{-1} g \pmod q的(f, g)有多解，但如果现在我多加一个条件，只需要其中f和g乘起来绝对值最小的一组解（当然只是一个栗子，实际的界会复杂一点），那么同时符合这两个条件的解的数量就会大幅减少，只要遍历里面少数的解就可解出生成密钥所用的(f, g)
剧透一下，上面的条件f, g < \sqrt{q/2}就是用来限定他是一组最小的解，或者说这是一个用来判断未知数(f, g)是否“小”的一个界
现在问题是，即使我知道(f, g)是一个最小的解，我如何根据h \equiv f^{-1} g \pmod q把这个最小解找出来呢，方法就是这里说的格规约，比如我可以用以下代码找出来：
M = matrix(ZZ, [
[1, h],
[0, q]
])
L = M.LLL()
f, g = L[0]
f = abs(f)
g = abs(g)
R = GF(q)
assert R(h) == R(f)^(-1) * R(g)
print('f = %d' % f)
print('g = %d' % g)
PS：以上代码需要在SageMath中运行，如果你实在不想用SageMath也可以用Python的fpylll替换上面代码，但会比较麻烦
前置知识
虽然我也想马上就解释上面代码的原理，但这不是一两句话就可以说清楚的，所以还是从前置知识开始，以下的内容基本是教科书的内容，如果你都知道，可以直接跳到下一章
线性代数知识
向量
向量可以理解为是空间上的一个点，或者从空间原点指向这个点的一个箭头，通常用空间上这个点的坐标表示这个向量，比如(1, 2)就是第一个轴的坐标是1，第二个轴的坐标是2的向量
在运算时通常不会把向量的坐标全写出来，而会附一个变量，一般表示向量的变量会用粗体字母或上加箭头的字母表示，如\mathbf{v}或\vec{v}
向量还有横向量和列向量的区分，可以简单理解为横向量就是在纸上横着写的向量（？），列向量就是在纸上竖着写的向量（？？），如果你是零基础的话，这个留到后面说到矩阵的时候会好理解一点。
横向量在代码上会好用一点，而列向量会在数学中好用一点，由于后面我会有大量的代码和数学表达混用，所以需要区分一下横向量和列向量
目前没想到特别好的表达方法，所以后面我会用粗体字母表示横向量，如\mathbf{v}，用上加箭头字母表示列向量，如\vec{v}，同一个字母表示同一个向量，只是不同表示，即\mathbf{v}^\text{T} = \vec{v}，\vec{v}^\text{T} = \mathbf{v}
PS：我想我应该会比较多用横向量，然后\mathbf{v}^\text{T}表示列向量
关于向量的运算，比较简单的如加法运算，加法运算需要两个向量有相同数量的元素，比如
\mathbf{a} = (a_0, a_1, \cdots, a_{n-1}) \\
\mathbf{b} = (b_0, b_1, \cdots, b_{n-1})
则
\mathbf{a} + \mathbf{b} =
(a_0 + b_0, a_1 + b_1, \cdots, a_{n-1} + b_{n-1})
另一个复杂一点的运算叫点乘，点乘其实就是把两个向量对应位置的元素相乘，然后把所有相乘结果加起来，即（在一些文章中也有用尖括号像
<\mathbf{a}, \mathbf{b}>表示点乘，看到的时候认识就好）
\mathbf{a} \cdot \mathbf{b} = \sum_{i=0}^{n-1} a_i \cdot b_i
如果把向量看成是矩阵的话，点乘也可以看成是前面向量的横向量乘上后面向量的列向量
\mathbf{a} \cdot \mathbf{b} =
\mathbf{a} \mathbf{b}^\text{T} =
\begin{bmatrix} a_0 & \cdots & a_{n-1} \end{bmatrix}
\begin{bmatrix} b_0 \\ \vdots \\ b_{n-1} \end{bmatrix} =
\sum_{i=0}^{n-1} a_i \cdot b_i
以后看到这样的矩阵相乘也想到对应位置相乘再相加就好了，方面后面向矩阵的学习过度
向量还有个长度，前面说了向量可以看成是空间上的一个点，长度就是这个点到空间原点的距离，这个长度在不同的范数（Norm）中的值和大小关系会不同，常见的范数如
在这里我用常用的L2-Norm，然后省略下标2简写为
\|\mathbf{v}\| = \sqrt{\mathbf{v} \cdot \mathbf{v}}
矩阵
上面说了，向量可以看成矩阵，而矩阵其实也可以看成同元素数量向量的堆叠，通常用大写字母表示一个矩阵，如
M = \begin{bmatrix}
1 & 2 & 3 \\
4 & 5 & 6 \\
\end{bmatrix}_{2 \times 3}
可以看成是两个元素数量为
3的横向量叠在一起，也可以看成是三个元素数量为
2的列向量叠在一起，通常我会习惯在右下角写下这个数量的行数和列数，比如
2 \times 3就是矩阵
M有两行三列
矩阵的运算有加法和乘法，加法需要运算的两个矩阵有相同的规模，即行数和列数都相同，然后把对应位置的元素加起来就好了
乘法会复杂一点，如果用通用一点的式子表达的话，就是，令
\begin{array}{cc}
A = \begin{bmatrix}
a_{0, 0} & \cdots & a_{0, q-1} \\
\vdots & \ddots & \vdots \\
a_{p-1, 0} & \cdots & a_{p-1, q-1} \\
\end{bmatrix}_{p \times q}
&
B = \begin{bmatrix}
b_{0, 0} & \cdots & b_{0, r-1} \\
\vdots & \ddots & \vdots \\
b_{q-1, 0} & \cdots & b_{q-1, r-1} \\
\end{bmatrix}_{q \times r}
\end{array}
则
AB = \begin{bmatrix}
\sum_{i=0}^{q-1} a_{0, i}b_{i, 0} & \cdots & \sum_{i=0}^{q-1} a_{0, i}b_{i, r-1} \\
\vdots & \ddots & \vdots \\
\sum_{i=0}^{q-1} a_{p-1, i}b_{i, 0} & \cdots & \sum_{i=0}^{q-1} a_{p-1, i}b_{i, r-1} \\
\end{bmatrix}_{p \times r}
实际上就是
A对应行的行向量和
B对应列的列向量做点乘，生成
AB中对应行列的元素，可以看上面通用表达自行体会
需要注意的是，做矩阵乘法运算的两个矩阵，需要前面矩阵的列数等于后面矩阵的行数，然后乘法结果的矩阵的行数等于前面矩阵的行数，列数等于后面矩阵的列数
理所当然地，一般情况下矩阵乘法是不可逆的，即AB \ne BA，否则上面行列关系就不会符合
线性组合
在格攻击中，最重要的一步是构造线性方程然后构造格基，在讲线性方程前线看看线性组合
从定义上来看，给定一堆向量\mathcal{B} = \{\mathbf{b_0}, \mathbf{b_1}, \cdots, \mathbf{b_{n-1}}\}（当然可以看成一个矩阵，但严谨一点我用集合表示），如果存在一堆系数a_0, a_1, \cdots, a_{n-1}，使得\mathbf{v} = \sum_{i=0}^{n-1} a_i \mathbf{b_i}，则称\mathbf{v}是\mathcal{B}的线性组合，抽象点看，就是向量\mathbf{v}可以仅使用\mathcal{B}中的向量，通过加法和数乘组合得到
若\mathcal{B}中某个向量是其他向量的线性组合，则称\mathcal{B}线性相关，否则称\mathcal{B}线性无关，只是两个定义，一般来说我们会更想要构造出线性无关
接下来看看线性方程，假设有下面一组线性方程
a_{0, 0} x_0 + a_{0, 1} x_1 + \cdots + a_{0, n-1} x_{n-1} = y_0 \\
a_{1, 0} x_0 + a_{1, 1} x_1 + \cdots + a_{1, n-1} x_{n-1} = y_1 \\
\vdots \\
a_{m-1, 0} x_0 + a_{m-1, 1} x_1 + \cdots + a_{m-1, n-1} x_{n-1} = y_{m-1} \\
如果是普通的解方程的话，接下来就是消元，但是上面说了，一个方程多个未知数是不能通过消元解的，所以需要把问题逐渐过渡到格的问题，或者至少是线性代数的问题
首先提取出每组相同的x_i，再把a_{i, j}和y_j组合成（列）向量
\begin{bmatrix} a_{0, 0} \\ a_{1, 0} \\ \vdots \\ a_{m-1, 0} \end{bmatrix} x_0
+
\begin{bmatrix} a_{0, 1} \\ a_{1, 1} \\ \vdots \\ a_{m-1, 1} \end{bmatrix} x_1
+ \cdots +
\begin{bmatrix} a_{0, n-1} \\ a_{1, n-1} \\ \vdots \\ a_{m-1, n-1} \end{bmatrix} x_{n-1}
=
\begin{bmatrix} y_0 \\ y_1 \\ \vdots \\ y_{m-1} \end{bmatrix}
最后把
x_i也组合成（列）向量
\begin{bmatrix}
a_{0, 0} & a_{0, 1} & \cdots & a_{0, n-1} \\
a_{1, 0} & a_{1, 1} & \cdots & a_{1, n-1} \\
\vdots & \vdots & \ddots & \vdots \\
a_{m-1, 0} & a_{m-1, 1} & \cdots & a_{m-1, n-1} \\
\end{bmatrix}_{m \times n}
\begin{bmatrix} x_0 \\ x_1 \\ \vdots \\ x_{n-1} \end{bmatrix}_{n \times 1}
=
\begin{bmatrix} y_0 \\ y_1 \\ \vdots \\ y_{m-1} \end{bmatrix}_{m \times 1}
这样就可以得到用矩阵和向量表示的线性方程组，如果把这个方程表示成
A \vec{x} = \vec{y}的话，向量
\vec{y}就是矩阵
A中所有列向量的线性组合，组合的系数是
\vec{x}中的元素
也可以把这个线性方程写成横向量的表达，但记得矩阵需要转置
\begin{array}{l}
\begin{bmatrix}
x_0 & x_1 & \cdots & x_{n-1}
\end{bmatrix}
\begin{bmatrix}
a_{0, 0} & a_{1, 0} & \cdots & a_{m-1, 0} \\
a_{0, 1} & a_{1, 1} & \cdots & a_{m-1, 1} \\
\vdots & \vdots & \ddots & \vdots \\
a_{0, n-1} & a_{1, n-1} & \cdots & a_{m-1, n-1} \\
\end{bmatrix}_{n \times m} \\
=
\begin{bmatrix}
y_0 & y_1 & \cdots & y_{m-1}
\end{bmatrix}
\end{array}
格知识
从定义上讲，令\mathbf{b_0}, \mathbf{b_1}, \cdots, \mathbf{b_{n-1}}为空间\mathbb{R}^m上一组线性无关的向量（m \ge n），则以下集合被称作一个格：
\mathcal{L}(\mathbf{b_0}, \mathbf{b_1}, \cdots, \mathbf{b_{n-1}}) :=
\{\sum_{i=0}^{n-1} x_i \mathbf{b_i}\ :\ x_i \in \mathbb{Z}\}
其中矩阵
B = \begin{bmatrix} \mathbf{b_0} & \mathbf{b_1} & \cdots & \mathbf{b_{n-1}} \end{bmatrix} 被称作这个格的格基，所以又可表示为
\mathcal{L}(B) :=
\{\mathbf{B}\mathbf{x}\ :\ \mathbf{x} \in \mathbb{Z}^n\}
简单来讲就是对格基
B的整数线性组合的集合，看个图可能好理解一点
在格的世界中有一系列的问题，而这里解方程只会涉及到一个，叫（Search‑）SVP问题（Search Shortest Vector Problem），即寻找一个格中的最短非零向量
很不幸的是，高维格的SVP是个困难问题，所以给我们留下的只有两个选择：只解低维格的SVP，或者解apprSVP问题（Approximate SVP）
关于apprSVP，只需要知道可以使用LLL或类似的格规约算法找到格中近似最短的非零向量即可，建议使用SageMath或者fpylll的LLL算法
PS：在后面格攻击中，如果把LLL算法的界代进去的话会发现可以攻击的界非常糟糕，但不知道是LLL的界太松还是SageMath有优化的问题，实际应用中发现效果还是非常好嘀，所以关于LLL，不用管界直接怼就是了
附SageMath中的参考代码：
M = matrix(ZZ, [[1, h], [0, q]])
L = M.LLL()
# or
L = M.BKZ(block_size=2)
shortestV = L[0]
一般用BKZ效果会比LLL好，其中的block_size设置越大效果越好，最大是格的维度
关于LLL算法和BKZ算法，这里我不会讲，给以后留个坑（逃
格攻击理论
这里的格攻击其实是一个启发式，所谓启发式可以大概理解为，你不能证明它是对的，但实际用起来有很好的效果的一些算法
其中用到的最重要的一个东西是高斯启发式
当然不是直接用的，下面一步步细说，首先介绍几个工具
第一个是基本域，设B = [\mathbf{b_0}, \mathbf{b_1}, \cdots, \mathbf{b_{n-1}}]为一组格基，则这组格基的基本域为
\mathcal{F}(\mathbf{b_0}, \cdots, \mathbf{b_{n-1}}) :=
\{ t_0 \mathbf{b_0} + \cdots + t_{n-1} \mathbf{b_{n-1}}\ :\ 0 \le t_i < 1 \}
一般在全文只有一组格基的时候会简写为
\mathcal{F}，即假设你知道
\mathcal{F}对应什么格基
从图像上看，\mathcal{F}就是这个格里原点最近的其中“一格”，看图，比如二维中的栗子
后面会用到的是基本域的容积，从结果上看，基本域的容积等于格基的行列式（的绝对值），证明略，可以参考[HPS14]
\text{Vol}(\mathcal{F}) = \text{det}(\mathcal{L}(B))
第二个工具是
超球体，其实就是高维空间的球体，同样也需要用到它的容积
设\mathbb{B}_R(\mathbf{a})为一个以点\mathbf{a}为圆心，半径为R的n维球体，那么\mathbb{B}_R(\mathbf{a})的容积（或者说体积）为
这个式子有点复杂，所以一般除了维度较小的情况，不大会用这个式子，而是用它的估算值，在维度n足够大时，可以作以下估算
高斯说，一个以零点为圆心的超球体中格点的数量可以用格的基本域的容积估算
把上面式子中向量数量定为1，即可估算只含一个格点的超球半径，即上面高斯启发式的图中的\sigma(L)
PS：“只含一个格点”其实只是一个估算，实际上因为\|\mathbf{v}\| = \|-\mathbf{v}\|，所以随便圈一个球里面应该都有偶数个非零格点
逆向思维一下就是说，长度短于\sigma(L)的格向量（约）只有一个，所以如果我们通过计算得到一个向量长度（约）小于\sigma(L)，那么大概率可以确定它是一个最短向量，在维度不高的情况下，可以用LLL把他求出来
最后需要注意的是，这个方法只是一个启发式方法，在实际操作中能不能求到这个最短向量还需要考虑几个东西
- 格的维度不能太高，不然LLL的准确度不能保证，而且LLL会用很长时间，然后是上面的容积估算也不能保证
- 格需要“足够随机”，上面图也说了，这只是一个概率估算，所以需要足够随机，后面会说到一些不太随机的格基，比如背包的格和NTRU的格，在实际操作中，这两个问题的格基维度其实都不能太大
最后，虽然它是一个启发式，但如果是比赛赛题的话，就是一定能解的，不妨先怼一遍，不然再想能不能降维打击或者用其他方法
格攻击应用
有了上面一大堆知识就可以讲应用了，先看回上面引言中的栗子：
有条件f, g < \sqrt{q/2}、\text{gcd}(f, qg) = 1和h \equiv f^{-1} g \pmod q，已知(q, h)，求(f, g)
直接使用同余式构造格基好像有点困难，因为LLL规约使用的格基里面的元素都需要在\mathbb{Z}上，所以不妨先把他展开，展开前先整理一下，因为现在我只知道f是小未知数，不能保证f^{-1}也是小未知数，所以先消去逆，两边乘一下f即可，得到
fh \equiv g \pmod q
然后可以插入一个未知的
k，把同余式展开，习惯上我会把一个小的未知数放一边（如这里的
g），方便后续构造格基
hf - qk = g
构造格基的方法其实就是构造出线性方程组，如果你是像我这样把小未知数放等号右边的话，那么一般是等号左边有多少未知数就需要构造多少组方程
现在等号左边有(f, k)两个未知数，所以还需要多一个方程，在没有更多的条件可用的情况下，可以插入一些恒等式，比如f = f，然后就可以得到
PS：加入恒等式时，尽量把格基构造成上三角或下三角矩阵，方便求行列式
1 \cdot f - 0 \cdot k = f \\
h \cdot f - q \cdot k = g
把这组线性方程用矩阵-向量的形式表达出来，这里我习惯用横向量
\begin{bmatrix} f & -k \end{bmatrix}
\begin{bmatrix}
1 & h \\
0 & q \\
\end{bmatrix}
=
\begin{bmatrix} f & g \end{bmatrix}
为了方便后面表述，不妨令
\mathbf{v} = (f, -k)、
B = \begin{bmatrix} 1 & h \\ 0 & q \end{bmatrix}、
\mathbf{w} = (f, g)，然后上面方程就可以简写为
\mathbf{v} B = \mathbf{w}
首先很明显，\mathbf{v}中的f和-k都是整数，根据格的定义，\mathbf{w}是格基B中向量的整数线性组合，所以\mathbf{w}是以B为格基的格\mathcal{L}(B)中的一个格点
接下来我会期望\mathbf{w}是\mathcal{L}(B)的最短向量，这就要用到上面说的估算方法，首先计算一下\mathbf{w}的长度\|\mathbf{w}\|和高斯的期望值\sigma(\mathcal{L}(B))
\|\mathbf{w}\| = \sqrt{f^2 + g^2} < \sqrt{\frac{q}{2} + \frac{q}{2}} = \sqrt{q} \\
\sigma(\mathcal{L}(B)) \approx \text{det}(\mathcal{L}(B))^{1/2} = \sqrt{q} \\
于是就有
\|\mathbf{w}\| \lessapprox \sigma(\mathcal{L}(B))，根据前面分析，
\mathbf{w}大概率为
\mathcal{L}(B)的最短向量，LLL解之
先来说说LLL，LLL是一个格规约算法，即把原来格基规约成一个向量更短的新格基，注意这两个格基都表达同一个格，即同一个格可以有多个不同的格基，具体可参考[HPS14]的7.4节，里面会介绍如何用一个行列式绝对值为1的U矩阵去换基，以及证明
以下是LLL规约后的格基的一些属性
从结果上看，LLL规约后的格基的第一个向量即我们想要的近似最短向量，可以不妨先拿出来当作\mathbf{w} = (f, g)，或者也可能是-\mathbf{w} = (-f, -g)，因为前面也说了，两者长度相等，但我们知道(f, g)都是正数，所以取出来后取个绝对值就好了
然后如无意外到这里就可以从|\mathbf{w}|中提取(f, g)，最后根据具体的题目条件验证(f, g)是否正确解就好
于是根据上面思路就得到引言中的EXP代码
配平
如果仅使用上面应用中的知识的话，在实际应用中会发现构造出来的格基很多都不会符合条件，其实还需要结合一个配平的技巧
还是先看一个栗子
bits = 512
while True:
q = random_prime(2^bits, lbound=2^(bits - 1))
R = GF(q)
f = R(random_prime(2^(3*bits//4 - 1)))
g = R(random_prime(2^(bits//4 - 1)))
if gcd(f, q*g) == 1:
h = f^(-1) * g
break
print('q = %d' % q)
print('h = %d' % h)
# q = 8114289192306045529468503816747459961749065015867824745830651826013976103158819451256549284151983553594151248484475469795886232064196980904146067875445697
# h = 838300814476084941348853046958667696551309880234911679429950542076545172518536842358788741627265361068496523434290150119253124253856647527559973827885148
# Aim: find f and g
和上一个栗子类似，只是大小关系发生了变化，不妨令bits为\ell，粗略估算一下，首先是f和g的大小，因为是随机选取，所以大概率是
f \approx 2^{\frac{3}{4} \ell - 1} \\
g \approx 2^{\frac{1}{4} \ell - 1}
所使用的线性方程还是上面的
\mathbf{v} B = \mathbf{w}，接下来估算
\|\mathbf{w}\|和
\sigma(\mathcal{L}(B))的关系
\|\mathbf{w}\| \gtrapprox 2^{\frac{3}{4} \ell - 1}
\sqrt{\frac{q}{2}}
= \text{det}(\mathcal{L}(B))
\approx \sigma(\mathcal{L}(B)) \\
大概率没有
\|\mathbf{w}\| \le \sigma(\mathcal{L}(B))，即大概率
\mathbf{w}不是
\mathcal{L}(B)的短向量
观察一下原因，这个情况是由\mathbf{w}的“不平衡”造成的，即\mathbf{w} = (f, g)中f比g（数量级）大，从而f把\mathbf{w}给“拉长”了
解决的方法是把线性方程中的g扩大到数量级和f相同，首先令D = 2^{\frac{\ell}{2}}，这个D可以使得f \approx D \cdot g
接下来修改线性方程，把与g相乘的一列都乘上D
\begin{bmatrix} f & -k \end{bmatrix}
\begin{bmatrix}
1 & D h \\
0 & D q \\
\end{bmatrix}
=
\begin{bmatrix} f & D g \end{bmatrix}
不妨令这个新方程为
\mathbf{v} B' = \mathbf{w}'，可以自行验证这个方程的正确性
粗略计算一下这个改动会发生什么
\sigma(\mathcal{L}(B')) \approx \sqrt{Dq} > \sqrt{q} \approx \sigma(\mathcal{L}(B)) \\
\|\mathbf{w}'\| \approx 2^{\frac{3}{4} \ell} \approx \|\mathbf{w}\|
可以发现
\sigma(\mathcal{L}(B'))增大了很多，但
\|\mathbf{w}'\|只增大了一点，几乎没有影响
接下来重新进行估算
\|\mathbf{w}'\| = \sqrt{f^2 + (Dg)^2} \le 2^{\frac{3}{4}\ell - \frac{1}{2}} \\
\sigma(\mathcal{L}(B')) \approx \text{det}(\mathcal{L}(B'))^{1/2} = \sqrt{Dq} \approx 2^{\frac{3}{4}\ell} \\
配平后就有
\|\mathbf{w}'\| \lessapprox \sigma(\mathcal{L}(B'))，LLL解之即可
注意最后解出的是\mathbf{w}' = (f, Dg)，要获得g的话记得除回D
参考代码：
D = 2^(512//2)
M = matrix(ZZ, [[1, h*D], [0, q*D]])
L = M.LLL()
f, g = L[0]
f = abs(f)
g = abs(g//D) # div D
R = GF(q) # check
assert R(h) == R(f)^(-1) * R(g)
print('f = %d' % f)
print('g = %d' % g)
总结
最后总结一下使用格解小未知数方程的套路：
- 找小未知数，构造线性方程、短向量和格基
- 配平，直到\|\mathbf{w}\| \lessapprox \sigma(\mathcal{L}(B))
- 如果配不到说明方程找错了/这题不能用格规约解
- LLL解短向量，注：维度不能过高
- 短向量中提取未知数
- 代回题目验证
栗子
栗子部分我打算举一些常见的可以用格攻击解决的经典问题，其中应该会插一些上面没讲但在栗子中会用到的东西
Wiener攻击
首先在两年前的Wiener's v.s Lattices中我写过Wiener攻击和某一种格攻击是等价的，但这篇文章写的有点乱（毕竟是两年前写的），下面在这里复述一下用格做Wiener的方法
Legendre’s theorem
首先Wiener攻击主要用到的是Legendre’s theorem
所以原理部分只用证可以用格攻击解Legendre’s theorem的问题就好了，把问题抽象一下即给定a， b，和条件
|\frac{a}{b}-\frac{c}{d}| \le \frac{1}{2d^2}
解出
c和
d（只是Legendre’s theorem换了一种表述）
为了拆开绝对值，不妨假设a/b \ge c/d，左右同乘bd，为了方便计算，令其左边为s
s = ad-bc \le \frac{b}{2d}
于是就可以构造格基（后面我为了方便都用圆括号表示横向量了）
(d, c)
\begin{bmatrix} 1 & a \\ 0 & -b \\ \end{bmatrix}
= (d, s)
明显
(d, s)不平衡，所以还需要配平，需要一个
D \approx \frac{d}{s}，具体怎么算这里就略了，反正这个约等的界挺松的
然后配平
(d, c)
\begin{bmatrix} 1 & D a \\ 0 & - D b \\ \end{bmatrix}
= (d, D s)
不妨令方程为
\mathbf{v} B = \mathbf{w}，粗略估算
\|\mathbf{w}\| \approx \sqrt{2d^2} \\
\sigma(\mathcal{L}(B)) \approx \sqrt{Db} \approx \sqrt{\frac{db}{s}} > \sqrt{\frac{db}{\frac{b}{2d}}} = \sqrt{2d^2}
就有
\|\mathbf{w}\| \lessapprox \sigma(\mathcal{L}(B))，LLL解之
PS：其实这里假设了s \le d，实际还有可能s > d，需要另一种配平方式，可以参考原文
PPS：在二维格基的情况下，LLL是可以解SVP问题的，所以LLL解出来的向量可以看作是最短非零向量
代入Wiener
下面把上面用格攻击做Legendre’s theorem的方法代进Wiener攻击中再推一下
首先在RSA密钥生成中有
ed - k\varphi = 1
而
\varphi = (p-1) (q-1) = N - (p + q - 1)
令
s = (p + q - 1)即
\varphi = N - s
把上面所有东西整合在一起，即
ed - kN = ks + 1
然后就构造格基
(k, d)
\begin{bmatrix} 1 & -N \\ 0 & e \\ \end{bmatrix}
= (k, ks+1)
然后配平，其实只用配个
D \approx s即可，一般的RSA中，
p \approx q，所以一般配
D = N^{1/2}，如果是
p, q不平衡的RSA就要另外考虑了
(kD, d)
\begin{bmatrix} D & -N \\ 0 & e \\ \end{bmatrix}
= (kD, ks+1)
不妨令方程为
\mathbf{v} B = \mathbf{w}，估算
\|\mathbf{w}\| \approx \sqrt{2k^2D^2}\\
\sigma(\mathcal{L}(B)) \approx \sqrt{D e}
现在需要
\|\mathbf{w}\| \lessapprox \sigma(\mathcal{L}(B))，也即需要
2 k^2 D^2 \lessapprox De \\
2 k^2 N^{1/2} \lessapprox e \\
k^2 \lessapprox \frac{e}{2 N^{1/2}}
这里有一个比较麻烦的
k需要先消去，翻回上面，
k是在展开同余式的时候产生的
ed \equiv 1 \pmod {\varphi} \\
ed - k\varphi = 1
整理一下
k = \frac{ed - 1}{\varphi} > \frac{ed}{N}
代回去就是
\frac{e^2 d^2}{N^2} < k^2 \lessapprox \frac{e}{2 N^{1/2}}
整理得
e d^2 \lessapprox \frac{1}{2} N^{3/2}
时可以实施攻击
注意看，和Wiener攻击的条件d < \frac{1}{3} N^{1/4}相比，这是一个更广义的条件，因为实际上即使用连分数做攻击，e的大小也会有影响，而Wiener当时的文章为了省事直接令e \approx N才推出了一个狭义的结果，可以代入e \approx N看看两个界是否一样
PS：由于我上面计算用了些估算，所以前面的系数可能不太一样-
参考攻击代码：
D = 2^(n.nbits()//2)
m = matrix(ZZ, [
[D, n+1],
[0, -e]
])
L = m.LLL()
w = L[0]
v = m.solve_left(w)
k = abs(v[0])
d = abs(v[1])
phi = (e*d-1) // k
p_plus_q = n + 1 - phi
p_min_q = (p_plus_q^2 - 4*n)^(1/2)
p = (p_plus_q + p_min_q) // 2
q = n // p
assert p*q == n
print('p = %s' % p)
print('q = %s' % q)
print('d = %s' % d)
01背包密码
又简称背包密码或者Knapsack，先看栗子代码
from secret import flag
from random import getrandbits
import libnum
BITS = 1024
fb = bin(libnum.s2n(flag))[2:].zfill(8*len(flag))
fb = [int(i) for i in fb]
M = [getrandbits(BITS) for i in range(len(fb))]
S = sum([fb[i]*M[i] for i in range(len(fb))])
print('M = %s' % M)
print('S = %s' % S)
简单来说，现在知道S = \sum_{i=0}^{n-1}x_iM_i，知道M_i和S，求x_i，其中x_i \in \{0, 1\}
首先构造格基，过程我就省略了，直接上结果
(x_0, x_1, ..., x_{n-1}, -1)
\begin{pmatrix}
2 & 0 & 0 & \cdots & 0 & M_0 \\
0 & 2 & 0 & \cdots & 0 & M_1 \\
0 & 0 & 2 & \cdots & 0 & M_2 \\
\vdots & \vdots & \vdots & \ddots & \vdots & \vdots \\
0 & 0 & 0 & \cdots & 2 & M_{n-1} \\
1 & 1 & 1 & \cdots & 1 & S \\
\end{pmatrix}_{(n+1)*(n+1)} \\
= (2 x_0-1, 2 x_1-1, ..., 2 x_{n-1}-1, 0)\\
不妨令方程为
\mathbf{v} B = \mathbf{w}，类似地做个估算，不妨设一个
M的下界
M_\text{min}，使得所有
M_i \ge M_\text{min}
然后先计算行列式，这里可以搜箭形/爪形的行列式解法，比如：这
\text{det}(\mathcal{L}(B))
= | S - \frac{1}{2}\sum_{i=0}^{n-1}M_i |
= | \sum_{i=0}^{n-1} (x_i - \frac{1}{2}) M_i |
\ge \frac{1}{2} n M_\text{min}
然后就可以估算长度
\|\mathbf{w}\| = \sqrt[n+1]{\sum_{i=0}^{n-1} (2x_i-1)^2}
= \sqrt[n+1]{\sum_{i=0}^{n-1} 1}
= n ^ \frac{1}{n+1} \\
\sigma(\mathcal{L}(B))
\ge \sqrt{n} \text{det}(\mathcal{L}(B))^{\frac{1}{n}}
\ge \sqrt{n} (\frac{1}{2} n M_\text{min})^\frac{1}{n}
\gtrapprox M_\text{min}^\frac{1}{n}
所以当
n ^ \frac{1}{n+1} \le M_\text{min}^\frac{1}{n}时
\|\mathbf{w}\| \le \sigma(\mathcal{L}(B))，问题可解，不过这个条件是挺容易实现的，就不用每次都算一遍了
不过如上文所说的，背包的格不是随机的，所以虽然条件符合，但维度一大就会解不出来
PS：可以对格基最后一列配个大系数，但实测作用不大，原理未明
参考代码：
with open('./out', 'r') as f:
exec(f.read())
n = len(M)
B = zero_matrix(ZZ, n+1)
for i in range(n):
B[i, i] = 2
B[-1, i] = 1
B[i, -1] = M[i]
B[-1, -1] = S
L =B.LLL()
print(L[0]) # 理论上只含(-1, 0, 1)且只有最后一位是0，否则是没解出
import libnum
fb = ''.join([str((li+1)//2) for li in L[0][:-1]])
flag = libnum.n2s(int(fb, 2))
print(flag) # 若解出 2 xi - 1
fb = ''.join([str((li+1)//2) for li in -L[0][:-1]])
flag = libnum.n2s(int(fb, 2))
print(flag) # 若解出 1 - 2 xi
Learning with errors
简称LWE问题，栗子就不写了，可以直接看我很久之前写的Embedding Technique里的栗子
当时写Embedding Technique的时候没有写证明，所以这里主要写条件证明部分
简单来说，LWE的问题是有条件A_{m \times n} \vec{x}_{n \times 1} + \vec{e}_{m \times 1} \equiv \vec{b}_{m \times 1}，m < n，知道A、\vec{b}，求\vec{x}或\vec{e}，注意这里都是列向量
由Embedding Technique可以构造下面线性方程，先设一个e的上界e_\text{max}，使得所有e_i \le e_\text{max}，注意下面粗体是横向量，\mathbf{0}^T是全0的列向量
(-\mathbf{x}_{1 \times n}\ |\ 1)
\begin{bmatrix}
(A^T)_{n \times m} & (\mathbf{0}^T)_{n \times 1} \\
\mathbf{b}_{1 \times m} & e_\text{max} \\
\end{bmatrix}_{(n+1) \times (m+1)} \\
= (\mathbf{e}_{1 \times m}\ |\ e_\text{max})
不妨令方程为
\mathbf{v} B = \mathbf{w}，同样地先计算行列式，由于这里的格基不一定是方阵，所以算行列式会有点麻烦
首先非方阵格基的行列式公式
代公式，然后还要利用行列式性质和Cofactor Formula（可以参考[Str22]）
\begin{aligned}
\text{det}(B^TB)
&=
\begin{vmatrix}
\begin{bmatrix}
A & \mathbf{0} \\
\mathbf{b}^T & e_\text{max} \\
\end{bmatrix}
\begin{bmatrix}
A^T & \mathbf{0}^T \\
\mathbf{b} & e_\text{max} \\
\end{bmatrix}
\end{vmatrix} \\
&=
\begin{vmatrix}
A A^T & e_\text{max} \mathbf{b}^T \\
e_\text{max} \mathbf{b} & e_\text{max}^2 \\
\end{vmatrix} \\
&= e_\text{max}^2
\begin{vmatrix}
A A^T & \mathbf{b}^T \\
\mathbf{b} & 1 \\
\end{vmatrix} \\
&=
e_\text{max}^2 (
\begin{vmatrix}
A A^T & \mathbf{0}^T \\
\mathbf{b} & 1 \\
\end{vmatrix}
+
\begin{vmatrix}
A A^T & \mathbf{b} \\
\mathbf{b}^T & 0 \\
\end{vmatrix} ) \\
&\ge e_\text{max}^2
\begin{vmatrix}
A A^T & \mathbf{0}^T \\
\mathbf{b} & 1 \\
\end{vmatrix} \\
&=
e_\text{max}^2 \cdot \text{det}(A A^T)
\end{aligned}
所以
\text{det}(\mathcal{L}(B))
= \sqrt{\text{det}(B^TB)}
\ge \sqrt{e_\text{max}^2 \cdot \text{det}(A A^T)}
= e_\text{max} \cdot \sqrt{\text{det}(A A^T)}
然后估算长度
\|\mathbf{w}\| = \sqrt[m+1]{\sum_{i=0}^{m-1} e_i^2 + e_\text{max}^2}
\le \sqrt[m+1]{(m+1) e_\text{max}^2}
\le e_\text{max}^{\frac{2}{m+1}} \\
\sigma(\mathcal{L}(B))
\ge \sqrt{n} \text{det}(\mathcal{L}(B))^{\frac{1}{n}}
\ge \sqrt{n} (e_\text{max} \cdot \sqrt{\text{det}(A A^T)})^\frac{1}{n}
\gtrapprox (e_\text{max} \cdot \sqrt{\text{det}(A A^T)})^\frac{1}{n}
所以当
e_\text{max}^{\frac{2}{m+1}} \le (e_\text{max} \cdot \sqrt{\text{det}(A A^T)})^\frac{1}{n}时
\|\mathbf{w}\| \le \sigma(\mathcal{L}(B))时可解，因为
\mathbf{e}是一个误差向量，所以里面元素都很小，所以这也是一个很容易符合的条件，而且如果不是特别构造，
A算是随机矩阵，解的效果会比背包好很多
PS：但是算吐了-
参考代码：参考Embedding Technique里的代码
NTRU
全称：N-th degree Truncated polynomial Ring Units
建议直接看之前的NTRU学习笔记
Ring-LWE
或者叫RLWE，其实就是环上的LWE问题
建议直接看之前的2023D3CTF的d3bdd中RLWE的部分
Hidden Number Problem
隐藏数问题，又叫HNP
建议直接看之前的HNP学习笔记
参考
- [HPS14] Hoffstein J, Pipher J, Silverman J H, et al. An Introduction to Cryptography[M]. Springer New York, 2014.
- [Str22] Strang G. Introduction to linear algebra[M]. Wellesley-Cambridge Press, 2022.
- [MG02] Micciancio D, Goldwasser S. Complexity of lattice problems: a cryptographic perspective[M]. Springer Science \& Business Media, 2002.