应用LLL算法解决背包问题

应用LLL算法解决背包问题

$Instruction$

背包问题，从实际生活上来描述就是，假定一个背包容重$W$，现有$n$个物品，其重量分别为$a_1,a_2,\cdots ,a_n$，需要知道装哪些物品能刚好将背包装满，用数学符号描述即是
$$a_1{x}_1+a_2{x}_2+\cdots +a_n{x}_n=W$$
问题需要求得的即是$x_i$的值；而根据问题的条件限制，比如 x i ∈ { 0 , 1 } x_i\in\{0,1\} xi​∈{0,1}时，该背包问题即为0-1背包问题；而 x i ∈ { 0 , 1 , ⋯   , n } x_i\in\{0,1,\cdots,n\} xi​∈{0,1,⋯,n}时，该背包问题为有界背包问题；当$x_i\ge 0$时，该背包问题为无界背包问题；关于背包问题的分类还有很多种，之后的例子中我们会将矩阵乘法问题转换为有界背包问题来求解

而LLL算法实际上是求解格上的SVP问题（最短向量），那么要将LLL算法应用于解决背包问题，我们需要通过基向量构造一个合适的格，并且验证我们需要求得的$x_i$在构造的向量在该格上并且是该格的最短向量

那么我们先将背包问题以矩阵的形式进行表示：
$$c={xx}_{1\times n}\cdot {aa}_{n\times 1}$$
其中$a_{n\times 1}=\left(\begin{array}{c}{a}_1\\ a_2\\ a_3\\ ⋮\\ a_n\end{array}\right)$，$x_{1\times n}=\left(\begin{array}{c}{x}_1{x}_2{x}_3\cdots x_n\end{array}\right)$，结果为整数$c$；当然，也可以转置过来构造为$c={aa}_{1\times n}\cdot {xx}_{n\times 1}$，但是为了之后的例子方便书写，我们用前者的构造方式

那么为了将背包问题与格相联系起来，我们构造一个矩阵乘法等式
$$KKAA=YY$$
其中$AA$是一个$(n+1)\times (n+1)$维的矩阵，为
$$A=\left(\begin{array}{c}{II}_{n\times n}{aa}_{n\times 1}\\ 0_{1\times n}-c\end{array}\right)$$
另外，$YY$赋值为
$$YY=\left(\begin{array}{c}{x}_1{x}_2{x}_3\cdots x_{n}0\end{array}\right)$$
$KK$为
$$KK=\left(\begin{array}{c}{x}_1{x}_2{x}_3\cdots x_{n}1\end{array}\right)$$
代入$KKAA=YY$中进行矩阵乘法后可以发现满足该等式

到这里我们成功将带有$x_i$的行向量用矩阵乘法的方式表示出来，可以这么说$AA$作为一个格的基向量组合，通过线性变换（左乘$KK$）得到行向量$YY$；那么行向量$YY$是在由基向量$AA$组成的格上

成功构造出一个让有$x_i$的向量在已知的格上，那么这个向量还必须满足是该格的最短向量，也即是SVP问题的描述

格的最短向量长度根据高斯启发式，满足这样的约等式
$$\sqrt[]{\frac{n}{2\pi e}}(|L|{)}^{\frac{1}{n}}$$
其中$n$为构造格的由基向量组成的矩阵的维度；在实际应用中我们会用到这个式子来验证我们所构造的$YY$是否为最短向量

通过LLL算法求得的矩阵的第一行即是$YY$；而为什么是结果矩阵的第一行可能是LLL算法内部的原因

Python实现

def solve(suq_a, c, n):
    A = Matrix(ZZ, n + 1, n + 1) # 构造一个(n+1)x(n+1)维的矩阵
    for i in range(n):
        A[i, i] = 1 
    for i in range(n):
        A[i, n] = suq_a[i] 
    A[n, n] = -c
    res = A.LLL()[0] # 取结果矩阵的第一行数据
    return res

$Example$

$Problem$

来源于2022 羊城杯CTF_Crypto linearAlgebra

from secret import flag
import libnum
import random
import os

assert flag[:7] == b'DASCTF{' and flag[-1] == ord('}')

n = 32
asize = 512

def pad(m, n):
  assert len(m) < n*n
  return m + b'\x00' + os.urandom(n*n - len(m) - 1)

def bytes2Matrix(m, n):
  assert len(m) == n*n
  return matrix(ZZ, [[m[r*n+c] for c in range(n)] for r in range(n)])

def randMatrix(s, n):
  while True:
    Mt = matrix(ZZ, [[random.randint(0, 2^s) for c in range(n)] for r in range(n)])
    if Mt.det() != 0:
      return Mt

def corrupt(Mt, s, n):
  Mt_ = matrix(ZZ, Mt)
  for i in range(n):
    r_ = random.randint(0, n-1)
    c_ = random.randint(0, n-1)
    Mt_[r_, c_] = random.randint(0, 2^s)
  return Mt_

m = pad(flag[7:-1], n)
M = bytes2Matrix(m, n)
A = randMatrix(asize, n)
C = A*M
Ac = corrupt(A, asize, n)

print('C:')
print(C)
print('Ac:')
print(Ac)

$Analysis$

简单来说，题目是将flag以单个字符的形式填入$32\times 32$维的矩阵，从第一行第一列的元素开始从左到右从上到下依次填入，剩余的部分以随机字节的形式也以相同方式依次填入矩阵，该矩阵作为$MM$；同时生成，生成随机数（大小满足$(0,2^{512})$）填入同样维数为$32\times 32$的一个新的矩阵$AA$；使用矩阵乘法计算出$CC=AAMM$；再对$AA$中随机32个元素进行随机替换，生成$AcAc$

到此为止，问题的计算过程就结束了；已知量有$CC$和$AcAc$，求解$MM$的部分元素（第一行以及第二行的部分元素，由于flag的长度有限，很可能只占有第一行以及第二行的部分元素，剩余元素均是随机生成的字节）

由于矩阵乘法的特殊性，每行元素依次乘以每列元素的和，其实可以将其看做是一个背包问题；也就是结果矩阵的每一个元素的生成过程实际上都是背包问题的数学表达，所以我们可以将这里的矩阵乘法问题以元素为单位转换为背包问题

由于已知的是$AA$经过数据扰乱后的$AcAc$，但是数据扰乱是随机选取位置进行的，而该矩阵总的元素量有32 * 32个；那么很有可能矩阵中很多行的元素没有变化，这里以行为单位看元素没有变化是因为，在矩阵乘法$CC=AAMM$中，矩阵$AA$是以每行的元素依次乘以矩阵$MM$的每列的元素

那么如果我们在矩阵$AcAc$中找到整行里的元素都没有被改变的第$k$行，该行在矩阵乘法$CC=AAMM$中会乘以$MM$的每列元素，生成结果矩阵$CC$中的第$k$行，而我们需要的是矩阵$MM$的第一行以及第二行的部分元素（其实就是$MM$的每列元素的第一个以及第二个）；将需要知道的$MM$的每个元素的生成过程都看作是一个背包问题，那么已知相当于背包问题的$a_i$序列（也就是$AcAc$的第$k$行的所有元素）以及背包问题的结果$c$；通过以上条件我们就可以应用LLL算法来求解该问题，得到结果$MM$的每一列元素的大小

那么关键又落在了如何知道矩阵$AcAc$的哪一行元素和$AA$是一样的，其实也可以跳过此步骤但是这里为了进一步体现LLL算法的用法，我们选择先判断哪一行元素没有变化；而判断某行元素有没有变化，实际上也可以说是求解$MM$的本身的一部分

我们取矩阵$CC$的第一列元素（也可以是其他列的元素），该列元素的生成过程是，矩阵$AA$的每行元素乘以矩阵$MM$的第一列元素得到的结果；那么我们将矩阵$CC$的第一列按元素为单位来看，那么就是32个背包问题（行数是32），使用矩阵$AcAc$应用相同的方法求解得到每个背包问题的结果；正确结果应该是矩阵$MM$的第一列元素，并且大小都在0~255之间（因为是每个元素实际上是单个字节转换为整数），那么如果出现不同的结果（这里是所有元素都远大于255），即是矩阵$AcAc$中的该行元素遭到替换，与原始矩阵$AA$的不同，所以导致经过LLL算法后的结果也不同；而出现正确结果，则说明矩阵$Ac$的该行元素与原始矩阵$AA$一致

$Solution$

# type: ignore
C = matrix(ZZ, 32, 32, c)
Ac = matrix(ZZ, 32, 32, Ac)

def solve(seq_a, s, n):
    A = Matrix(ZZ, n + 1, n + 1)
    for i in range(n):
        A[i, i] = 1
    for i in range(n):
        A[i, n] = seq_a[i] 
    A[n, n] = -s
    res = A.LLL()[0]
    return res

for k in range(32):
    seq_a = Ac[k]
    c_i = C[k][0]
    n = 32
    nbit = n
    res = solve(seq_a, c_i, n)
    if 0 < res[0] < 255:
        seq_a = Ac[k]
        for i in range(32):
            c_i = C[k][i]
            flag_i = solve(seq_a, c_i, n) # 得到的是M矩阵的某一列的元素
            print(chr(flag_i[0]),end="") # 取某一列的第一个元素
        for i in range(10):
            c_i = C[k][i]
            flag_i = solve(seq_a, c_i, n)
            print(chr(flag_i[1]),end="") # 取某一列的第二个元素
        break

---

---

那么我们接下来来验证之前提到的高斯启发式

import gmpy2
n = 32
seq_a = Ac[2] # 对比发现矩阵Ac与矩阵A第三行相同
s = C[2][0]
A = Matrix(ZZ, n + 1, n + 1)
for i in range(n):
    A[i, i] = 1
for i in range(n):
    A[i, n] = seq_a[i] 
A[n, n] = -s
res = A.LLL()
# 高斯启发式
print(float(sqrt((n + 1) / (2 * pi * e))) * gmpy2.iroot(res.det(), n + 1)[0])
print((int(sqrt((n + 1) / (2 * pi * e))) * gmpy2.iroot(res.det(), n + 1)[0]).bit_length())
# 求向量长度
tmp = 0
for i in res[0]:
    tmp += i ^ 2
print(gmpy2.iroot(tmp,2)[0])
print((gmpy2.iroot(tmp,2)[0]).bit_length())
'''
81803.830770217814
16
104181
17
'''

发现几乎相同