原文链接：https://tover.xyz/p/HSSP-note/
看看能不能抢到2024年的最后一篇帖子
以下正文：

现在是个人都会用正交格攻击了，搞得我不学好像就落后似的，所以抽空学习了一下

其中，和正交格相关的最出名的就是HSSP问题了，于是下面就把HSSP问题怼一遍

HSSP问题

HSSP（Hidden Subset Sum Problem）问题大概如下

令M为大整数，整数\alpha_1, \cdots, \alpha_n \in \mathbb{Z}_M，向量\pmb{x_1}, \cdots, \pmb{x_n} \in \mathbb{Z}^m为m维向量，且\pmb{x_i}的元素落在\{0, 1\}中，令
\pmb{h} = (h_1, \cdots, h_m) \equiv \sum_{i=1}^n \alpha_i \pmb{x_i} \pmod M
现知道M和\pmb{h}，求\pmb{\alpha} = (\alpha_1, \cdots, \alpha_n)和\pmb{x_1}, \cdots, \pmb{x_n}

PS：根据Coron的测试数据，M的大小大概是M = O(2 \iota n^2 + n \cdot log\ n)，其中\iota = 0.035，（也有m>n）

和经典背包问题（SSP）的区别是，在HSSP中隐藏（Hidden）了\alpha_1, \cdots, \alpha_n，所以无法直接通过构造格来解

给一个生成问题的样例代码

def genHssp(m, n, M):
  R.<z> = PolynomialRing(Zmod(M))
  x = [R([randint(0, 1) for mi in range(m)]) for ni in range(n)]
  a = [randint(0, M-1) for ni in range(n)]
  h = sum([a[i] * x[i] for i in range(n)])
  return (a, [xi.list() for xi in x]), (M, h.list())

m = 200
n = 100
M = 199999999999997
(a, X), (M, h) = genHssp(m, n, M)

线性代数知识

在讲正交格前先来补充一点线性代数的前置知识，更详细的内容可以参考

Strang, Gilbert. Introduction to linear algebra. Wellesley-Cambridge Press, 2022.

讲到线性代数，第一个想到的应该都是像
A \pmb{x} = \pmb{b}
这样的矩阵和向量的运算

如果不把矩阵和向量看成单纯的数字，而是看成是空间和空间中的点的话，就可以得到著名的四大基本子空间：行空间、列空间、零空间（核空间）和左零空间

给定一个m \times n（m行n列）的矩阵A，可以把A的每一列看成空间的基，用这个行向量的基张成（Span，即进行线性组合）的空间就叫列空间，数学方式表示大概是（如果是格的话就是\forall \pmb{x} \in \mathbb{Z}^n）
\{b \in \mathbb{R}^m\ |\ A \pmb{x} = b,\ \forall \pmb{x} \in \mathbb{R}^n\}
如果把A的每一行看作空间的基，那么张成的空间就叫行空间，数学表示大概是
\{b \in \mathbb{R}^n\ |\ \pmb{x} A = b,\ \forall \pmb{x} \in \mathbb{R}^m\}
如果只关注A \pmb{x} = \pmb{0}这个方程，那么方程的所有解落在零空间（又叫核空间）中
\{\pmb{x} \in \mathbb{R}^n\ |\ A \pmb{x} = \pmb{0}\}
如果把\pmb{x}放在A的左边，得到的空间又叫左零空间
\{\pmb{x} \in \mathbb{R}^m\ |\ \pmb{x} A = \pmb{0}\}
令r为矩阵A的秩（Rank），那么列空间和行空间的维度都是r，零空间的维度是n-r，左零空间的维度是m-r

直观上看，A消元后非零的列和行的数量都是r，零的列是n-r，零的行是m-r，详细的证明可以看书或者网上找找

在四大基本子空间中有一个重要的结论是，行空间与零空间相互垂直，列空间与左零空间相互垂直

可以简单证明一下，令\pmb{b} = \pmb{x}_b A为行空间的一个向量，令\pmb{x}为零空间的一个向量，那么两个向量相乘
\pmb{b} \cdot \pmb{x} = \pmb{x}_b A \pmb{x}
根据零空间的性质，A \pmb{x} = \pmb{0}，所以
\pmb{b} \cdot \pmb{x} = \pmb{x}_b \pmb{0} = 0
也就是任意一个行空间的向量与任意一个零空间的向量都相互垂直，即行空间与零空间相互垂直

列空间与左零空间的证明类似

另一个重要的结论是，行空间与零空间可以张成整个\mathbb{R}^n空间，列空间与左零空间可以张成整个\mathbb{R}^m空间

直观上看，行空间与零空间相互垂直就是不相关，然后两个空间的维度加起来刚好是n

列空间与左零空间的也类似

由这两个结论可得，如果要求一个格的正交格（就是相互垂直的），那么只要求他的零空间（行看作基）或者左零空间（列看作基）就好

Flatter

Flatter是一个比LLL更快的格规约算法

和LLL不同的是，目前SageMath没有原生集成Flatter，所以需要装一个

安装方法可以直接看Github，大概就是

git clone https://github.com/keeganryan/flatter.git
sudo apt install libgmp-dev libmpfr-dev fplll-tools libfplll-dev libeigen3-dev
mkdir build && cd ./build
cmake ..
make -j8
# 软链接路径改成自己的PATH
ln -s `pwd`/bin/flatter [imath:0]HOME/.local/bin
flatter -h

然后我直接抄了@Neobeo的做法，通过子进程调用Flatter二进制应用

# https://github.com/Neobeo/HackTM2023/blob/main/solve420.sage
# faster LLL reduction to replace `M.LLL()` wiith `flatter(M)`
def flatter(M):
    from subprocess import check_output
    from re import findall
    M = matrix(ZZ,M)
    # compile https://github.com/keeganryan/flatter and put it in [/imath:0]PATH
    z = '[[' + ']\n['.join(' '.join(map(str,row)) for row in M) + ']]'
    ret = check_output(["flatter"], input=z.encode())
    return matrix(M.nrows(), M.ncols(), map(int,findall(b'-?\\d+', ret)))

在SageMath中用的时候直接把正常用的M.LLL()换成flatter(M)即可

HSSP问题的格通常都比较大，所以用Flatter会比LLL节约不少时间

Nguyen-Stern算法

接下来就看看这个HSSP到底要怎么解，以下内容我参考了

Coron, Jean-Sébastien, and Agnese Gini. "A polynomial-time algorithm for solving the hidden subset sum problem." Annual International Cryptology Conference. Cham: Springer International Publishing, 2020.

还有文章对应的代码，和@tl2的文章，有很多被我忽略掉的内容都可以在这篇文章中看到

Nguyen和Stern的做法是，给定\pmb{h} \pmod M，首先找与\pmb{h}垂直的向量\pmb{u}，那么就有
\pmb{u} \cdot \pmb{h} \equiv \sum_{i=1}^n \alpha_i (\pmb{u} \cdot \pmb{x_i}) \equiv 0 \pmod M
令向量
\pmb{p_u} = ((\pmb{u} \cdot \pmb{x_1}), (\pmb{u} \cdot \pmb{x_2}), \cdots, (\pmb{u} \cdot \pmb{x_n}))
那么问题就可以转化为
\pmb{p_u} \cdot \pmb{\alpha} \equiv 0 \pmod M
也就是\pmb{p_u}和\pmb{\alpha}在模M的情况下相互垂直

然后如果\pmb{u}是短向量的话，那么\pmb{p_u}也会是短向量（因为\pmb{x_i}的元素落在\{0, 1\}中），如果\pmb{p_u}比所有与\pmb{\alpha}垂直的非零向量都短的话，那么就只能是\pmb{p_u} = \pmb{0}

而如果\pmb{p_u} = \pmb{0}的话，就是\pmb{u} \cdot \pmb{x_i} = 0，令L_x是以\pmb{x_1}, \cdots, \pmb{x_n}为基的格，就可以得到\pmb{u} \in L_x，即\pmb{u}是L_x的正交格L_x^\bot中的向量

L_x^\bot的维度是m-n：因为L_x^\bot的基的秩是r=n（n \le m），然后我这里是看成行向量为基的空间（行空间），且L_x的基是n行m列的，所以根据前面的线性代数知识，与行空间垂直的零空间的维度就是m-r = m-n

所以，如果我们可以找到m-n个满足条件的向量\pmb{u}的话，就相当于找到了L_x的正交格L_x^\bot，进而使用L_x^\bot找到L_x，最后由L_x恢复基\pmb{x_1}, \cdots, \pmb{x_n}

于是，最后得到的攻击路线就是

用\pmb{h}构造格基，LLL找到m-n个短向量\pmb{u}_i
用\pmb{u}_i构造格L_x^\bot，用L_x^\bot找L_x的正交补\bar{L_x}（可以看作是和L_x同一个空间，但基不是\pmb{x}_i）
对\bar{L_x}使用BKZ恢复\pmb{x}_i

Part.1 找短向量u

这里我直接用Coron论文中的方法造格

首先拆开
\pmb{u} \cdot \pmb{h} \equiv 0 \pmod M
得到
\sum_{i=1}^m u_i h_i \pmod 0 \pmod M
然后提出其中的u_1
u_1 h_1 + \sum_{i=2}^m u_i h_i \pmod 0 \pmod M
两边乘h_1^{-1} \pmod M
u_1 + \sum_{i=2}^m u_i (h_i h_1^{-1}) \pmod 0 \pmod M
最后拆开模数M，并换一下位置
kM + \sum_{i=2}^m u_i (-h_i h_1^{-1}) = u_1
根据这个关系就可以构造格基
B_1 = \begin{bmatrix} M & \\ -h_2 h_1^{-1} & 1 \\ -h_3 h_1^{-1} & & 1 & \\ \vdots & & & \ddots \\ -h_m h_1^{-1} & & & & 1 \end{bmatrix}_{m \times m}
令
\begin{aligned} \pmb{v}_1 &= (k, u_2, u_3, \cdots, u_m) \\ \pmb{w}_1 &= (u_1, u_2, u_3, \cdots, u_m) \end{aligned}
那么就是
\pmb{w}_1 \cdot B_1 = \pmb{w}_1
根据Coron文章第三章的分析，可以保证对B_1规约后的前m-n行是满足条件的向量\pmb{u}，这个，可以自己看论文...

B = matrix(ZZ, m)
B[0, 0] = M
h0i = Integer(h[0]).inverse_mod(M)
for i in range(1, m):
  B[i, 0] = - h[i] * h0i
  B[i, i] = 1
L = flatter(B)

vh = vector(Zmod(M), h)
print([vector(Zmod(M), list(l)) * vh  for l in L])

另外，还可以构造另一种更直观的格基
B_2 = \begin{bmatrix} M & \\ h_1 & 1 \\ h_2 & & 1 & \\ \vdots & & & \ddots \\ h_m & & & & 1 \end{bmatrix}_{(m+1) \times (m+1)}
令
\begin{aligned} \pmb{v}_2 &= (-k, u_1, u_2, u_3, \cdots, u_m) \\ \pmb{w}_2 &= (0, u_1, u_2, u_3, \cdots, u_m) \end{aligned}
那么就是
\pmb{w}_2 \cdot B_2 = \pmb{w}_2
这个格基在Coron的文章和@tl2的文章都有类似的，可以去参考一下

Part.2 恢复格Lx

这一步就比较简单

首先根据上面分析，用L的前m-n就可以构造L_x^\bot

然后只需要求L_x^\bot的零空间就可以得到L_x的正交补\bar{L_x}

这里我直接用SageMath的right_kernel求令空间，亲测把algorithm指定为pari的话会快一点

Lxo = matrix(ZZ, L[:m-n])
Lxc = Lxo.right_kernel(algorithm='pari').matrix() # faster
print('right_kernel done.')

Lx_real = matrix(ZZ, [xi + [0] * (m - len(xi)) for xi in X])
rsc = Lxc.row_space()
print([xi in rsc for xi in Lx_real])

Part.3 恢复xi

理论上直接对Lxc求个LLL或者BKZ就可以恢复\pmb{x_1}, \cdots, \pmb{x_n}，但实际上并没有

细看一下，\pmb{x_1}, \cdots, \pmb{x_n}的元素在\{0, 1\}中，这在01背包问题中也遇到过类似的问题，所以可以利用类似的解决方法，即把\pmb{x_1}, \cdots, \pmb{x_n}转化为
2\pmb{x_1}-\pmb{1}, 2\pmb{x_2}-\pmb{1}, \cdots, 2\pmb{x_n}-\pmb{1}
就可以把元素转化到(-1, 1)中

虽然这对向量长度影响不大，但乘上去的系数2会增大格基的行列式，就更容易筛掉无关的变量

于是就可以构造这样一个格基（其中E是元素全为1、大小和\bar{Lx}一样的矩阵）
B_3 = \begin{bmatrix} -E \\ \hline 2 \bar{Lx} \end{bmatrix}_{2n \times m}
令（U是n \times n，看作一种映射就好）
U \bar{L_x} = L_x
可以得到关系
[I_{n}, U]_{n \times 2n} \cdot B_3 = [2 U \bar{Lx} - E]_{n \times m} = [2L_x - E] = \begin{bmatrix} 2\pmb{x_1}-\pmb{1} \\ \vdots \\ 2\pmb{x_n}-\pmb{1} \end{bmatrix}
所以对B_3归约后就可能得到2\pmb{x_1}-\pmb{1}, 2\pmb{x_2}-\pmb{1}, \cdots, 2\pmb{x_n}-\pmb{1}

进一步观察发现其实B_3中的E的每一行都是相关的（甚至相同的），实际作用的就一行，对B_3规约后也发现有n-1行全为0

所以不妨令\pmb{e}为全为1的行向量，就可以把格简化为
B_4 = \begin{bmatrix} -\pmb{e} \\ \hline 2 \bar{Lx} \end{bmatrix}_{n+1 \times m}
参考代码

 def checkMatrix(M, wl=[-1, 1]):
  M = [list(_) for _ in list(M)]
  ml = list(set(flatten(M)))
  logging.debug(ml)
  return sorted(ml) == sorted(wl)
    
e = matrix(ZZ, [1] * m)
B = block_matrix([[-e], [2*Lxc]])
Lx = B.BKZ()
assert checkMatrix(Lx)
assert len(set(Lx[0])) == 1

最后恢复一下\pmb{x}_i和\alpha_i

Lx = Lx[1:]
E = matrix(ZZ, [[1 for c in range(Lxc.ncols())] for r in range(Lxc.nrows())])
Lx = (Lx + E) / 2

Lx2 = []
e = vector(ZZ, [1] * m)
rsc = Lxc.row_space()
for lx in Lx:
  if lx in rsc:
    Lx2 += [lx]
    continue
  lx = e - lx
  if lx in rsc:
    Lx2 += [lx]
    continue
  print('Something wrong?')
Lx = matrix(Zmod(M), Lx2)
vh = vector(Zmod(M), h)
va = Lx.solve_left(vh)

PS：其实用\begin{bmatrix} 2 \bar{Lx} \\ \hline -\pmb{e} \end{bmatrix}做格也可以，但是干扰的那一行就不会放在第一行，还要另外写代码找出来（就是全为1或者全为-1的行）

模板/参考代码

最后把上面所有的代码整合一下

import logging
logging.basicConfig(
    level=logging.DEBUG,
    format="[%(levelname)s] %(message)s"
)

# https://github.com/Neobeo/HackTM2023/blob/main/solve420.sage
# faster LLL reduction to replace `M.LLL()` wiith `flatter(M)`
def flatter(M, **kwds):
    from subprocess import check_output
    from re import findall
    M = matrix(ZZ,M)
    # compile https://github.com/keeganryan/flatter and put it in [imath:0]PATH
    z = '[[' + ']\n['.join(' '.join(map(str,row)) for row in M) + ']]'
    ret = check_output(["flatter"], input=z.encode())
    return matrix(M.nrows(), M.ncols(), map(int,findall(b'-?\\d+', ret)))

def genHssp(m, n, M):
  R.<z> = PolynomialRing(Zmod(M))
  x = [R([randint(0, 1) for mi in range(m)]) for ni in range(n)]
  a = [randint(0, M-1) for ni in range(n)]
  h = sum([a[i] * x[i] for i in range(n)])
  return (a, [xi.list() for xi in x]), (M, h.list())

def checkMatrix(M, wl=[-1, 1]):
  M = [list(_) for _ in list(M)]
  ml = list(set(flatten(M)))
  logging.debug(ml)
  return sorted(ml) == sorted(wl)

def Nguyen_Stern(h, m, n, M):
  B = matrix(ZZ, m)
  B[0, 0] = M
  h0i = Integer(h[0]).inverse_mod(M)
  for i in range(1, m):
    B[i, 0] = - h[i] * h0i
    B[i, i] = 1
  #L = B.BKZ()	# slooooooow
  L = flatter(B)
  logging.info('flatter done.')

  '''
  vh = vector(Zmod(M), h)
  logging.debug([vector(Zmod(M), list(l)) * vh  for l in L])
  '''

  Lxo = matrix(ZZ, L[:m-n])
  Lxc = Lxo.right_kernel(algorithm='pari').matrix() # faster
  logging.info('right_kernel done.')

  '''
  try:
    Lx_real = matrix(ZZ, [xi + [0] * (m - len(xi)) for xi in X])
    rsc = Lxc.row_space()
    logging.debug([xi in rsc for xi in Lx_real])
  except:
    pass
  '''

  e = matrix(ZZ, [1] * m)
  B = block_matrix([[-e], [2*Lxc]])
  Lx = B.BKZ()
  logging.info('BKZ done.')
  assert checkMatrix(Lx)
  assert len(set(Lx[0])) == 1

  Lx = Lx[1:]
  E = matrix(ZZ, [[1 for c in range(Lxc.ncols())] for r in range(Lxc.nrows())])
  Lx = (Lx + E) / 2

  Lx2 = []
  e = vector(ZZ, [1] * m)
  rsc = Lxc.row_space()
  for lx in Lx:
    if lx in rsc:
      Lx2 += [lx]
      continue
    lx = e - lx
    if lx in rsc:
      Lx2 += [lx]
      continue
    logging.warning('Something wrong?')
  Lx = matrix(Zmod(M), Lx2)

  vh = vector(Zmod(M), h)
  va = Lx.solve_left(vh)
  return Lx, va

# stolen from https://github.com/tl2cents/Implementation-of-Cryptographic-Attacks/blob/main/MultivariateHSSP/A%20Polynomial-Time%20Algorithm%20for%20Solving%20the%20Hidden%20Subset%20Sum%20Problem.ipynb
def derive_M(n):
    iota=0.035
    Mbits=int(2 * iota * n^2 + n * log(n,2))
    M = random_prime(2^Mbits, proof = False, lbound = 2^(Mbits - 1))
    return Integer(M)

m = 200
n = 100
M = derive_M(n)
(a, X), (M, h) = genHssp(m, n, M)
logging.debug('m: %d | n: %d' % (m, n))
logging.debug('%s, %s' % (M, M.nbits()))

Lx, va = Nguyen_Stern(h, m, n, M)
print(sorted(va) == sorted(a))

M的生成采用了Coron的 M = O(2 \iota n^2 + n \cdot log\ n)，偷懒了一下直接偷@tl2的derive_M