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

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

HSSP问题·

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

令MM为大整数，整数α1,⋯ ,αn∈ZM\alpha_1, \cdots, \alpha_n \in \mathbb{Z}_M，向量x1,⋯ ,xn∈Zm\pmb{x_1}, \cdots, \pmb{x_n} \in \mathbb{Z}^m为mm维向量，且xi\pmb{x_i}的元素落在{0,1}\{0, 1\}中，令

h=(h1,⋯ ,hm)≡∑i=1nαixi(modM)\pmb{h} = (h_1, \cdots, h_m) \equiv \sum_{i=1}^n \alpha_i \pmb{x_i} \pmod M

现知道MM和h\pmb{h}，求α=(α1,⋯ ,αn)\pmb{\alpha} = (\alpha_1, \cdots, \alpha_n)和x1,⋯ ,xn\pmb{x_1}, \cdots, \pmb{x_n}

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

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

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

1
2
3
4
5
6
7
8
9
10
11
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.

讲到线性代数，第一个想到的应该都是像

Ax=bA \pmb{x} = \pmb{b}

这样的矩阵和向量的运算

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

给定一个m×nm \times n（mm行nn列）的矩阵AA，可以把AA的每一列看成空间的基，用这个行向量的基张成（Span，即进行线性组合）的空间就叫列空间，数学方式表示大概是（如果是格的话就是∀x∈Zn\forall \pmb{x} \in \mathbb{Z}^n）

{b∈Rm ∣ Ax=b, ∀x∈Rn}\{b \in \mathbb{R}^m\ |\ A \pmb{x} = b,\ \forall \pmb{x} \in \mathbb{R}^n\}

如果把AA的每一行看作空间的基，那么张成的空间就叫行空间，数学表示大概是

{b∈Rn ∣ xA=b, ∀x∈Rm}\{b \in \mathbb{R}^n\ |\ \pmb{x} A = b,\ \forall \pmb{x} \in \mathbb{R}^m\}

如果只关注Ax=0A \pmb{x} = \pmb{0}这个方程，那么方程的所有解落在零空间（又叫核空间）中

{x∈Rn ∣ Ax=0}\{\pmb{x} \in \mathbb{R}^n\ |\ A \pmb{x} = \pmb{0}\}

如果把x\pmb{x}放在AA的左边，得到的空间又叫左零空间

{x∈Rm ∣ xA=0}\{\pmb{x} \in \mathbb{R}^m\ |\ \pmb{x} A = \pmb{0}\}

令rr为矩阵AA的秩（Rank），那么列空间和行空间的维度都是rr，零空间的维度是n−rn-r，左零空间的维度是m−rm-r

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

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

可以简单证明一下，令b=xbA\pmb{b} = \pmb{x}_b A为行空间的一个向量，令x\pmb{x}为零空间的一个向量，那么两个向量相乘

b⋅x=xbAx\pmb{b} \cdot \pmb{x} = \pmb{x}_b A \pmb{x}

根据零空间的性质，Ax=0A \pmb{x} = \pmb{0}，所以

b⋅x=xb0=0\pmb{b} \cdot \pmb{x} = \pmb{x}_b \pmb{0} = 0

也就是任意一个行空间的向量与任意一个零空间的向量都相互垂直，即行空间与零空间相互垂直

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

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

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

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

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

Flatter·

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

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

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

1
2
3
4
5
6
7
8
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 $HOME/.local/bin
flatter -h

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

1
2
3
4
5
6
7
8
9
10
# 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 $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的做法是，给定h(modM)\pmb{h} \pmod M，首先找与h\pmb{h}垂直的向量u\pmb{u}，那么就有

u⋅h≡∑i=1nαi(u⋅xi)≡0(modM)\pmb{u} \cdot \pmb{h} \equiv \sum_{i=1}^n \alpha_i (\pmb{u} \cdot \pmb{x_i}) \equiv 0 \pmod M

令向量

pu=((u⋅x1),(u⋅x2),⋯ ,(u⋅xn))\pmb{p_u} = ((\pmb{u} \cdot \pmb{x_1}), (\pmb{u} \cdot \pmb{x_2}), \cdots, (\pmb{u} \cdot \pmb{x_n}))

那么问题就可以转化为

pu⋅α≡0(modM)\pmb{p_u} \cdot \pmb{\alpha} \equiv 0 \pmod M

也就是pu\pmb{p_u}和α\pmb{\alpha}在模MM的情况下相互垂直

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

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

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

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

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

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

Part.1 找短向量u·

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

首先拆开

u⋅h≡0(modM)\pmb{u} \cdot \pmb{h} \equiv 0 \pmod M

得到

∑i=1muihi(mod0)(modM)\sum_{i=1}^m u_i h_i \pmod 0 \pmod M

然后提出其中的u1u_1

u1h1+∑i=2muihi(mod0)(modM)u_1 h_1 + \sum_{i=2}^m u_i h_i \pmod 0 \pmod M

两边乘h1−1(modM)h_1^{-1} \pmod M

u1+∑i=2mui(hih1−1)(mod0)(modM)u_1 + \sum_{i=2}^m u_i (h_i h_1^{-1}) \pmod 0 \pmod M

最后拆开模数MM，并换一下位置

kM+∑i=2mui(−hih1−1)=u1kM + \sum_{i=2}^m u_i (-h_i h_1^{-1}) = u_1

根据这个关系就可以构造格基

B1=[M−h2h1−11−h3h1−11⋮⋱−hmh1−11]m×mB_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}

令

v1=(k,u2,u3,⋯ ,um)w1=(u1,u2,u3,⋯ ,um)\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}

那么就是

w1⋅B1=w1\pmb{w}_1 \cdot B_1 = \pmb{w}_1

根据Coron文章第三章的分析，可以保证对B1B_1规约后的前m−nm-n行是满足条件的向量u\pmb{u}，这个，可以自己看论文…

1
2
3
4
5
6
7
8
9
10
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])

另外，还可以构造另一种更直观的格基

B2=[Mh11h21⋮⋱hm1](m+1)×(m+1)B_2 = \begin{bmatrix} M & \\ h_1 & 1 \\ h_2 & & 1 & \\ \vdots & & & \ddots \\ h_m & & & & 1 \end{bmatrix}_{(m+1) \times (m+1)}

令

v2=(−k,u1,u2,u3,⋯ ,um)w2=(0,u1,u2,u3,⋯ ,um)\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}

那么就是

w2⋅B2=w2\pmb{w}_2 \cdot B_2 = \pmb{w}_2

这个格基在Coron的文章和@tl2的文章都有类似的，可以去参考一下

Part.2 恢复格Lx·

这一步就比较简单

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

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

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

1
2
3
4
5
6
7
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就可以恢复x1,⋯ ,xn\pmb{x_1}, \cdots, \pmb{x_n}，但实际上并没有

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

2x1−1,2x2−1,⋯ ,2xn−12\pmb{x_1}-\pmb{1}, 2\pmb{x_2}-\pmb{1}, \cdots, 2\pmb{x_n}-\pmb{1}

就可以把元素转化到(−1,1)(-1, 1)中

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

于是就可以构造这样一个格基（其中EE是元素全为11、大小和Lxˉ\bar{Lx}一样的矩阵）

B3=[−E2Lxˉ]2n×mB_3 = \begin{bmatrix} -E \\ \hline 2 \bar{Lx} \end{bmatrix}_{2n \times m}

令（U是n×nn \times n，看作一种映射就好）

ULxˉ=LxU \bar{L_x} = L_x

可以得到关系

[In,U]n×2n⋅B3=[2ULxˉ−E]n×m=[2Lx−E]=[2x1−1⋮2xn−1][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}

所以对B3B_3归约后就可能得到2x1−1,2x2−1,⋯ ,2xn−12\pmb{x_1}-\pmb{1}, 2\pmb{x_2}-\pmb{1}, \cdots, 2\pmb{x_n}-\pmb{1}

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

所以不妨令e\pmb{e}为全为11的行向量，就可以把格简化为

B4=[−e2Lxˉ]n+1×mB_4 = \begin{bmatrix} -\pmb{e} \\ \hline 2 \bar{Lx} \end{bmatrix}_{n+1 \times m}

参考代码

1
2
3
4
5
6
7
8
9
10
11
 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

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

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
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：其实用[2Lxˉ−e]\begin{bmatrix} 2 \bar{Lx} \\ \hline -\pmb{e} \end{bmatrix}做格也可以，但是干扰的那一行就不会放在第一行，还要另外写代码找出来（就是全为1或者全为-1的行）

模板/参考代码·

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

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
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
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 $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ιn2+n⋅log n)M = O(2 \iota n^2 + n \cdot log\ n)，偷懒了一下直接偷@tl2的derive_M