SekaiCTF 2025 Writeup #2 - unfairy-ring

Writeup for unfairy-ring from SekaiCTF 2025.

TLDR: Following the approach proposed by Miura¹ and Cheng², with some careful variable management, we can transform the equations into

\left\{ \begin{align*} x_1^2 + L_1(x_2, \cdots, x_n)x_1 + Q_1(x_2, \cdots, x_n) &= t_1 \\ x_2^2 + L_2(x_3, \cdots, x_n)x_2 + Q_2(x_3, \cdots, x_n) &= t_2 \\ &\vdots \\ x_m^2 + L_m(x_{m+1}, \cdots, x_n)x_m + Q_m(x_{m+1}, \cdots, x_n) &= t_m \\ \end{align*} \right.

Then randomly set $x_{m+1}, \cdots, x_n$ and solve the equations from $x_m$ to $x_1$ .

Introduction#

1
from functools import reduce
2
from uov import uov_1p_pkc as uov # https://github.com/mjosaarinen/uov-py/blob/main/uov.py
3
import os
4
FLAG = os.getenv("FLAG", "SEKAI{TESTFLAG}")
5

6
def xor(a, b):
7
    assert len(a) == len(b), "XOR inputs must be of the same length"
8
    return bytes(x ^ y for x, y in zip(a, b))
9

10
names = ['Miku', 'Ichika', 'Minori', 'Kohane', 'Tsukasa', 'Kanade']
11
pks = [uov.expand_pk(uov.shake256(name.encode(), 43576)) for name in names]
12
msg = b'SEKAI'
13

14
sig = bytes.fromhex(input('Ring signature (hex): '))
15
assert len(sig) == 112 * len(names), 'Incorrect signature length'
16

17
t = reduce(xor, [uov.pubmap(sig[i*112:(i+1)*112], pks[i]) for i in range(len(names))])
18
assert t == uov.shake256(msg, 44), 'Invalid signature'
19

20
print(FLAG)

This challenge is adapted from defund’s fairy-ring, which is a 6 key ring signature scheme based on UOV. The original solution exploits the key reuse vulnerability and solves the linear equation $f(x)-f(x+\delta)=y$ for an arbitrary quadratic function $f$ .

During DiceCTF, I was working on nil-circ and didn’t have chance to check the fairy-ring code. There was also some discussion about solving it without key reuse after the CTF ended, perhaps Neobeo can comment on it.

So it is quite impressive when Neobeo said that he can solve 7 key case without key reuse and will make a revenge in SekaiCTF. The only hints I got from Neobeo are these two screenshots:

hint1 hint2

The actual solution is far more complicated than the hints and includes many tricks. So in the writeup I’ll try to start from the very beginning and explain how we finally get to the 6 keys solution step by step.

Preliminaries#

Some mathematical tools we’ll use later.

Quadratic Forms over Characteristic 2 Field#

Def. A quadratic form is given by $f_{A}(\mathbf{x}) = \sum_{i, j}a_{ij}x_ix_j = \mathbf{x}^\top A\mathbf{x}$ , where $A=(a_{ij})$ . In a characteristic 2 field, the matrix $M$ is not required to be symmetric.

Here, we only consider the homogeneous quadratic forms because the signature scheme doesn’t involve any linear terms.

Def. The polar form of a quadratic form is $f'_{A}(\mathbf{x}, \mathbf{y}) = f_{A}(\mathbf{x}+\mathbf{y}) - f_{A}(\mathbf{x}) - f_{A}(\mathbf{y})$ .

Lemma. The polar form is bilinear.

Expanding the polar form, we get

\begin{align*} f'_{A}(\mathbf{x}, \mathbf{y}) &= f_{A}(\mathbf{x}+\mathbf{y}) - f_{A}(\mathbf{x}) - f_{A}(\mathbf{y}) \\ &= (\mathbf{x}+\mathbf{y})^\top A (\mathbf{x}+\mathbf{y}) - \mathbf{x}^\top A \mathbf{x} - \mathbf{y}^\top A \mathbf{y} \\ &= \mathbf{x}^\top A \mathbf{y} + \mathbf{y}^\top A \mathbf{x} \\ &= \mathbf{x}^\top (A + A^\top) \mathbf{y} \end{align*}

You can see that the polar form is bilinear. The matrix $A + A^\top$ is symmetric and its diagonal is $0$ when over characteristic 2 field.

Prop. If the dimension of $A$ is odd, $A + A^\top$ will be singular over characteristic 2 field.

Proof omitted.

Def. If we can split a quadratic form into two disjoint parts, i.e. $f_A(\mathbf{x})=f_B(\mathbf{x}_B)+f_C(\mathbf{x}_C)$ , these two parts are orthogonal.

Def. A linear transformation of a quadratic form is given by $\mathbf{x} = T\mathbf{y}$ , where $T$ is an arbitrary matrix. After the transformation, the quadratic form becomes $f_A(\mathbf{x} = T\mathbf{y})=(T\mathbf{y})^\top A(T\mathbf{y})=\mathbf{y}^\top T^\top AT\mathbf{y}$ . The new matrix $A'=T^\top AT$ is congruent to the original matrix $A$ .

We don’t ask $T$ to be invertible or square, so the transformation is only injective from $\mathbf{y}$ to $\mathbf{x}$ . As a result, the new quadratic form $f_{A'}$ is a “subset” of the original quadratic form $f_A$ .

The last property of the characteristic 2 field is Frobenius endomorphism $F(x)=x^{2^r}$ . It is a linear map over characteristic 2 field and will be used in the Thomae Wolf’s approach later.

Ring UOV Signature Scheme#

Def. Function $\mathcal{F}\in\mathbb{F}^n\to\mathbb{F}^m$ is a multivariate quadratic map if every component of $\mathcal{F}$ is a quadratic form, i.e., $\mathcal{F}(\mathbf{x})=(f_{A_1}(\mathbf{x}), \cdots, f_{A_m}(\mathbf{x}))$ .

A detailed explanation of the UOV signature scheme can be found in defund’s post. Basically, UOV is a multivariate quadratic map with trapdoor that can solve equation $\mathcal{F}(\mathbf{x})=\mathbf{t}$ efficiently.

We can rewrite UOV into equations of the form

\renewcommand{\vec}[1]{\mathbf{#1}} f_{A_i}(\vec{x})=\vec{x}^\top A_i\vec{x}=t_i, 1\leq i \leq m

Def. $MQ^{m}(n)$ denotes a multivariate quadratic problem with $m$ equations and $n$ variables.

Then we can extend the definition to a ring signature scheme.

Def. $MQ^{m}(n_1, \cdots, n_k)$ denotes a ring signature scheme involving $k$ quadratic maps $\mathcal{F}_i\in\mathbb{F}^{n_i}\to\mathbb{F}^m$ , s.t. $\sum_{i=1}^k \mathcal{F}_i(\mathbf{x}_i) = \mathbf{t}$ , where $\mathbf{x}_i\in\mathbb{F}^{n_i}$ .

This challenge can be formulated as an instance of $MQ^{m}(n, \cdots, n)$ with $k=6, n=112$ and $m=44$ , and no trapdoor is present.

Thomae Wolf’s Method#

In the first hint, Neobeo gives two papers Furue et al.³ and Hashimoto⁴ and both papers explains Thomae Wolf’s idea. I recommend reading this slide from PQCrypto 2021.

It’s not relevant to the actual solution, but I’ve finished this section before the drama happened, so I’ll keep it here.

Theorem. (Thomae Wolf) For an underdetermined $MQ^m(n)$ problem, there’s a polynomial time algorithm to reduce it to $MQ^{n-\alpha}(n-\alpha)$ , where $\alpha=\lfloor\frac{n}{m}\rfloor-1$ .

The algorithm consists of three steps:

Find a linear transformation $T$ such that in $f_{M_i'}, 1\leq i\leq \alpha$ , the variables $x_1, x_2, \cdots, x_m$ only appear in squared terms.
Choose $x_{m+1},\cdots, x_n$ properly so that the equations become linear in $x_1^2, x_2^2, \cdots, x_m^2$ .
Apply Frobenius endomorphism to transform these equations into linear equations in $x_1, x_2, \cdots, x_m$ .

Here I’ll formulate the steps in polar form.

Step 1#

Firstly we’ll find a set of linear independent vectors $\mathbf{v}_1,\mathbf{v}_2,\cdots,\mathbf{v}_m$ such that

f_{M_l}'(\mathbf{v}_i, \mathbf{v}_j)=0, \forall 1\leq i<j\leq m, 1\leq l\leq\alpha

We can start with $\mathbf{v}_1=\mathbf{e}_1$ . For each $2\leq j\leq m$ , solve $f_{M_l}'(v_i, v_j)=0, 1\leq i\leq\alpha, l<j$ and get $\mathbf{v}_j$ .

T=(\mathbf{v}_1, \mathbf{v}_2, \cdots, \mathbf{v}_m)^\top

Extend $T$ to a complete basis and apply the transformation to the equations. As a result, the first $m$ variables in the first $\alpha$ equations will be linear in squared terms. i.e. it could be written as

f_{M_l}=\sum_{i=1}^ma^{(l)}_{ii}x_i^2+\sum_{i=1}^mL_i^{(l)}(x_{m+1},\cdots, x_n)x_i+Q^{(l)}(x_{m+1},\cdots, x_n)=t_l, 1\leq l\leq \alpha

where $L_i^{(l)}$ is linear and $Q^{(l)}$ is quadratic.

To ensure that $\mathbf{v}_m$ satisfies the polar form conditions and is linearly independent of the first $m-1$ vectors, $\alpha$ can’t be too large. Therefore, the condition $n\geq\alpha (m-1)+m$ is required.

Step 2#

In the second step, we want to eliminate the linear terms $L_i^{(l)}(x_{m+1},\cdots, x_n)$ . In our definition, quadratic forms don’t contain linear terms, hence we can directly set $x_j=0$ for $j>m$ . However, in the Thomae-Wolf algorithm, the linear terms are non-zero and it’s necessary to solve the equations $L_i^{(l)}(x_{m+1},\cdots, x_n)=0$ and eliminate them.

To ensure the existence of a solution, $n-m\geq\alpha m$ must hold. Finally, we get

\sum_{i=1}^ma^{(l)}_{ii}x_i=t_l', 1\leq l\leq \alpha

Step 3#

In $\mathbb{F}_{2^r}$ , we can raise it to $2^{r-1}$ power and obtain

\sum_{i=1}^m\left(a^{(l)}_{ii}\right)^{2^{r-1}}x_i=\left({t_l'}\right)^{2^{r-1}}

After elimination, we get a $MQ^{n-\alpha}(n-\alpha)$ problem. We need $n-m\geq\alpha m$ to ensure the step 1 and 2 has a solution, thus $\alpha\leq\lfloor\frac{n}{m}\rfloor-1$ .

Complexity of Guessing Analysis#

In Furue’s work, they improved Thomae Wolf’s algorithm by adding a guessing step. It is very helpful for improving the complexity of the algorithm.

For example, we can ignore one equation and solve the remaining $m-1$ equations. This would cause a guessing complexity of $\frac{1}{256}$ for each guessed equation. We also need the probability of sampling a low rank matrix, which is estimated by the following lemma.

Prop. For a uniformly distributed matrix $M\in\mathbb{F}_{256}^{m\times n}, n \geq m, \Pr[\text{rank}(M) \leq m-k] \approx \frac{1}{256^{k(n-m+k)}}.$

Although the proof is not rigorous, we can first fix the first $m-k$ rows of $M$ and it is likely to be linear independent, then estimate the probability of the remaining $k$ rows being linear dependent on the first $m-k$ rows. for each row, only $m-k$ entries are free, and the remaining $n-m+k$ entries are determined by the first $m-k$ rows. Therefore, the probability of each row being linear dependent on the first $m-k$ rows is $\frac{1}{256^{n-m+k}}$ . So ignoring the probability of the first $m-k$ rows being linear dependent, the probability of sampling a rank $m-k$ matrix is $\frac{1}{256^{k(n-m+k)}}$ .

In most cases of a quadratic form, $n$ is much larger than $m$ , thus it’s very unlikely to sample a low rank matrix.

Solution#

The solution is similar to the algorithm in Cheng’s paper². However, I wasn’t aware of it until two months later when I was writing this writeup. Luckily, it didn’t become a full implement-paper challenge because the orthogonal structure of ring signature is needed to solve the 6 keys case.

Iterative reduction#

We’ve got two basic operations to reduce the equations,

Replace $f_{A_i} = t_i$ with $f_{A_i}-cf_{A_j} = t_i - ct_j$ , where $c$ is a constant.
Apply a linear transformation $\mathbf{x} = T\mathbf{y}$ . $T$ may be non-square and $\mathbf{y}$ is a subspace of $\mathbf{x}$ .

With these two operations, we can reduce the number of equations iteratively.

Lemma. $\text{MQ}^{m}(n_1, \cdots, n_k)$ can be reduced to $\text{MQ}^{m-1}(n_1-m-1, \cdots, n_k)$ .

The plan is to use linear transformation to make one variable orthogonal. In general, it’s impossible to make all equations orthogonal, so we must cost some variables to achieve that.

For a linear transformation $T$ , let’s say the columns of $T$ are $\mathbf{t}_0, \cdots, \mathbf{t}_k$ . We can express the orthogonal condition by the polar form $f_{A_i}'(\mathbf{t}_0,\mathbf{t}_j)=0, 1\leq j\leq k$ . Hence if we randomly set $\mathbf{t}_0$ , the rest $\mathbf{t}_j$ are all in the kernel of linear equations $f_{A_i}'(\mathbf{t}_0,\mathbf{t}_j)=0$ . Since we’ve got $m$ equations, the dimension of the kernel is at least $n-m$ . Note that $f_{A_i}'(\mathbf{t}_0,\mathbf{t}_0)=0$ always holds, so the largest subspace we can get is $k=n-m-1$ .

Once we have an orthogonal variable, we can do the first operation to eliminate the square terms $x_0^2$ . Finally, we will get

\left\{ \begin{align*} x_0^2 + Q_1(x_1, \cdots) &= t_1 \\ Q_2(x_1, \cdots) &= t_2 \\ &\vdots \\ Q_m(x_1, \cdots) &= t_m \\ \end{align*} \right.

Now we can keep the last $m-1$ equations and keep reducing it. After finding a solution of $x_1, \cdots$ , $x_0$ can be easily solved by $x_0^2 = t_1 - Q_1(x_1, \cdots)$ and the Frobenius endomorphism ensures that $x_1$ always has a solution. The operations we’ve done don’t damage the orthogonal structure and only the first $n_1$ variables are affected during the reduction. Hence, the last $m-1$ equations is a $MQ^{m-1}(n_1-m-1, n_2, \cdots, n_k)$ problem.

Improving the reduction#

In the previous reduction, the main cost comes from the $m$ equations $f_{A_i}'(\mathbf{t}_0,\mathbf{x})=0$ . But these equations aren’t promised to be linear independent. Actually, the rank of the equations may be less than $m$ .

There’s two possible ways to lower the rank of the equations.

Choose good $\mathbf{t}_0$ .
Randomly do many reductions and hope that the rank is less than $m$ .

Choosing good $\mathbf{t}_0$ #

The linear equations can be written as a matrix whose rows are $\mathbf{t}_0^\top(A_i + A_i^\top)$ . If we can find a linear combination of the rows to be zero, the rank of the matrix will be less than $m$ . For odd $n$ , the matrix $(A_i + A_i^\top)$ is always singular, so simply letting $\mathbf{t}_0$ be the kernel will work. For even $n$ , we have to do some combinations to find a singular matrix. Consider $(A_i + A_i^\top)+t(A_j + A_j^\top), t\in \mathbb{F}_{256}$ , Its determinant is a polynomial of $t$ . If we can find a root of the polynomial, the matrix will be singular.

With the rank minimizing trick, each reduction only costs $m$ variables now. But it’s still unlikely to find $\mathbf{t}_0$ that’s in the kernel of two matrices.

Guessing low rank matrix#

What’s the probability of sampling a low rank matrix?

With a good $\mathbf{t}_0$ , it is $m-1$ independent equations. After exculding the trivial solution $\mathbf{t}_0$ , there’s $n-1$ variables in the equations. Hence the shape of the matrix is $(m-1, n-1)$ and we’ve shown that the probability of the rank being less than $m$ is $\frac{1}{256^{n-m+1}}$ . So only if $n-m$ is small, it is possible to sample a low rank matrix.

Lemma. $\text{MQ}^{m}(m, n_2, \cdots, n_k)$ can be reduced to $\text{MQ}^{m-1}(1, n_2, \cdots, n_k)$ .

This is the only case that guessing may get low rank matrix with prob $\frac{1}{256}$ .

Guessing equations#

Before reduction, we can remove some equations and guess them later. Every equation has $\frac{1}{256}$ probability to match the target, so we can at most remove $4$ equations at first, which gives $2^{32}$ complexity.

7 keys#

The whole reduction procedure can be described as follows:

1
n = 112
2
m = 42 # plus two guessed equations
3
k = 7
4

5
t = [n]*k
6
seq = [6,6,6,5,2,5,5,3,2,1,4,3,0,4,4,4,3,3,2,2,1,1,1,1,1,0,0,0,0,0,0,0,0,0,0,0]
7

8
for ind in seq:
9
    print(t, m, ind)
10
    if t[ind] < m:
11
        if t[ind] == 0:
12
            raise ValueError
13
        t[ind] = 0
14
        m -= 1
15
        continue
16
    elif t[ind] <= m+1:
17
        t[ind] = 0
18
        m -= 2
19
        continue
20
    else:
21
        t[ind] -= m
22
        m -= 1

seq is the sequence of key indices that we will reduce. It’s non-trivial to find the optimal sequence, so I just manually solved it. We need to guess 2 equations, which is $2^{16}$ attempts.

Implementation & optimization#

The backward step involves inverse linear transformation and solving $x_0^2 = t_1 - Q_1(x_1, \cdots)$ . It is the most time-consuming part when brute forcing the solutions. Since the inverse transformation is still linear, the only non-linearity is the square root, which happens only $m$ times. Hence the matrix of $Q_1(x_1, \cdots)$ can be compressed into a smaller size.

It is wrapped in a Transform class, which has compress method to compress the matrix when reduction and compress_backward method to do the backward step with compressed matrix.

1
class Transform:
2
    def __init__(self):
3
        pass
4

5
    def forward(self, problem: ProblemData):
6
        """
7
        Reduce the problem to a smaller size.
8
        """
9
        raise NotImplementedError
10

11
    def compress(self, Ts):
12
        """
13
        Compress the matrix for the backward step.
14
        """
15
        raise NotImplementedError
16

17
    def backward(self, x):
18
        """
19
        Basic backward.
20
        """
21
        raise NotImplementedError
22

23
    def compress_backward(self, x_compressed):
24
        """
25
        Do the backward step with compressed matrix.
26
        """
27
        raise NotImplementedError

After compression, the backward step can reach 3000 iter/s with 16 cores in pure sage implementation.

Final code for 7 keys#

1
from tqdm import tqdm
2
from dataclasses import dataclass
3
from functools import lru_cache
4
from concurrent.futures import ProcessPoolExecutor
5
import random, os
6
from uov import uov_1p_pkc as uov
7

8
F = GF(2**8, name='z', modulus=x^8 + x^4 + x^3 + x + 1)
9
z = F.gen()
10
R = PolynomialRing(F, 'xbar')
11
xbar = R.gen()
12

13
elms = [F.from_integer(i) for i in range(2**8)]
14

15
n = 112
16
m = 42
17

18
def gen(n, m):
19
    U = random_matrix(F, n, n)
20
    while U.det() == 0:
21
        U = random_matrix(F, n, n)
22
    pk = []
23
    for _ in range(m):
24
        M = random_matrix(F, n, n)
25
        M[:m, :m] = 0
26
        M = U.transpose() * M * U
27
        pk.append(M)
28
    return pk
29

30
def from_uov_vector(elms):
31
    m = uov.m
32
    bs = elms.to_bytes(m, 'big')
33
    vs = [F.from_integer(i) for i in bs]
34
    return vs
35

36
def from_uov(pk):
37
    v = uov.v
38
    m = uov.m
39
    Ms = [[[0] * n for _ in range(n)] for _ in range(m)]
40
    m1 = uov.unpack_mtri(pk, v)
41
    m2 = uov.unpack_mrect(pk[uov.p1_sz:], v, m)
42
    m3 = uov.unpack_mtri(pk[uov.p1_sz + uov.p2_sz:], m)
43

44
    for i in range(v):
45
        for j in range(i, v):
46
            vec = from_uov_vector(m1[i][j])
47
            for k in range(m):
48
                Ms[k][i][j] = vec[k]
49

50
    for i in range(v):
51
        for j in range(m):
52
            vec = from_uov_vector(m2[i][j])
53
            for k in range(m):
54
                Ms[k][i][j + v] = vec[k]
55

56
    for i in range(m):
57
        for j in range(i, m):
58
            vec = from_uov_vector(m3[i][j])
59
            for k in range(m):
60
                Ms[k][i + v][j + v] = vec[k]
61

62
    Ms = [matrix(F, M) for M in Ms]
63
    return Ms
64

65
def to_uov(sol):
66
    output = []
67
    for s in sol:
68
        output.append([x.to_integer() for x in s.list()])
69
    output = [bytearray(x) for x in output]
70
    return output
71

72
def F_sqrt(x):
73
    return x.square_root()
74

75
class EQ:
76
    def __init__(self, Ms: list, v):
77
        self.Ms = Ms
78
        self.v = v
79

80
    def evaluate(self, sigs):
81
        result = 0
82
        for M, sig in zip(self.Ms, sigs):
83
            result += sig * M * sig
84
        result += self.v
85
        return result
86

87
    def compress(self, Ts):
88
        new_Ms = []
89
        for M, T in zip(self.Ms, Ts):
90
            new_Ms.append(T.transpose() * M * T)
91
        return EQ(new_Ms, self.v)
92

93
    @property
94
    def ns(self):
95
        return [M.nrows() for M in self.Ms]
96

97
    def __mul__(self, other):
98
        other = F(other)
99
        new_Ms = [other * M for M in self.Ms]
100
        new_v = other * self.v
101
        return EQ(new_Ms, new_v)
102

103
    def __add__(self, other):
104
        assert isinstance(other, EQ)
105
        assert len(self.Ms) == len(other.Ms)
106
        new_Ms = [M + other.Ms[i] for i, M in enumerate(self.Ms)]
107
        new_v = self.v + other.v
108
        return EQ(new_Ms, new_v)
109

110
    def  __rmul__(self, other):
111
        other = F(other)
112
        return self.__mul__(other)
113

114
class ProblemData:
115
    def __init__(self, eqs: list[EQ]):
116
        self.eqs = eqs
117
        self.ns = eqs[0].ns
118
        self.k = len(self.ns)
119
        self.m = len(eqs)
120
        for i in range(1, self.m):
121
            assert eqs[i].ns == self.ns
122

123
    def check(self, sigs):
124
        vs = [eq.evaluate(sigs) for eq in self.eqs]
125
        return all(v == 0 for v in vs)
126

127
    def get_M(self, i, index):
128
        return self.eqs[i].Ms[index]
129

130
    def get_key(self, index):
131
        Ms = [self.get_M(i, index) for i in range(self.m)]
132
        return Ms
133

134
    def compress(self, Ts):
135
        new_eqs = [eq.compress(Ts) for eq in self.eqs]
136
        return ProblemData(new_eqs)
137

138
class Transform:
139
    def __init__(self):
140
        pass
141

142
    def forward(self, problem: ProblemData):
143
        raise NotImplementedError
144

145
    def compress(self, Ts):
146
        raise NotImplementedError
147

148
    def backward(self, x):
149
        raise NotImplementedError
150

151
    def compress_backward(self, x_compressed):
152
        raise NotImplementedError
153

154
class Transform0(Transform):
155
    def __init__(self, index):
156
        self.index = index
157
        self.eq = None
158

159
    def forward(self, problem: ProblemData):
160
        Ms = problem.get_key(self.index)
161
        Ms = [M + M.transpose() for M in Ms]
162
        assert all(M[0].is_zero() for M in Ms)
163
        m = len(Ms)
164
        eqs = problem.eqs[:]
165
        for i0 in range(m):
166
            if Ms[i0][0, 0] != 0:
167
                break
168
        eqs[0], eqs[i0] = eqs[i0], eqs[0]
169
        coeffs = [eqs[i].Ms[self.index][0, 0] for i in range(problem.m)]
170
        eqs[0] = eqs[0] * coeffs[0].inverse()
171
        for i in range(1, problem.m):
172
            eqs[i] = eqs[i] + (-coeffs[i]) * eqs[0]
173
        assert eqs[0].Ms[self.index][0, 0] == 1
174
        for i in range(problem.m):
175
            eqs[i].Ms[self.index] = eqs[i].Ms[self.index][1:, 1:]
176
        self.eq = eqs.pop(0)
177
        return ProblemData(eqs)
178

179
    def compress(self, Ts):
180
        self.compressed_eq = self.eq.compress(Ts)
181
        new_Ts = Ts[:]
182
        new_Ts[self.index] = block_matrix(F, [[identity_matrix(1), 0], [0, Ts[self.index]]])
183
        return new_Ts
184

185
    def backward(self, x):
186
        u = self.eq.evaluate(x)
187
        x0 = F_sqrt(u)
188
        v0 = x[self.index]
189
        v1 = vector(F, [x0] + v0.list())
190
        x_new = x[:]
191
        x_new[self.index] = v1
192
        return x_new
193

194
    def compress_backward(self, x_compressed):
195
        u = self.compressed_eq.evaluate(x_compressed)
196
        x0 = F_sqrt(u)
197
        v0 = x_compressed[self.index]
198
        v1 = vector(F, [x0] + v0.list())
199
        x_new = x_compressed[:]
200
        x_new[self.index] = v1
201
        return x_new
202

203
class Transform1(Transform):
204
    def __init__(self, index, T):
205
        self.index = index
206
        self.T = T
207

208
    def forward(self, problem):
209
        return ProblemData([self.apply_T(eq) for eq in problem.eqs])
210

211
    def apply_T(self, eq: EQ):
212
        Ms = eq.Ms[:]
213
        Ms[self.index] = self.T.transpose() * Ms[self.index] * self.T
214
        return EQ(Ms, eq.v)
215

216
    def compress(self, Ts):
217
        new_Ts = Ts[:]
218
        new_Ts[self.index] = self.T * Ts[self.index]
219
        return new_Ts
220

221
    def backward(self, x):
222
        x0 = x[self.index]
223
        x_new = x[:]
224
        x_new[self.index] = self.T * x0
225
        return x_new
226

227
    def compress_backward(self, x_compressed):
228
        return x_compressed
229

230
class Transform2(Transform):
231
    def __init__(self, index):
232
        self.index = index
233
        self.eq = None
234
        self.pad = 0
235

236
    def forward(self, problem: ProblemData):
237
        Ms = problem.get_key(self.index)
238
        m = len(Ms)
239
        self.pad = Ms[0].nrows() - 1
240
        eqs = problem.eqs[:]
241
        for i0 in range(m):
242
            if Ms[i0][0, 0] != 0:
243
                break
244
        eqs[0], eqs[i0] = eqs[i0], eqs[0]
245
        coeffs = [eqs[i].Ms[self.index][0, 0] for i in range(problem.m)]
246
        eqs[0] = eqs[0] * coeffs[0].inverse()
247
        for i in range(1, problem.m):
248
            eqs[i] = eqs[i] + (-coeffs[i]) * eqs[0]
249
        assert eqs[0].Ms[self.index][0, 0] == 1
250
        new_eqs = []
251
        for i in range(problem.m):
252
            eq = eqs[i]
253
            Ms = eq.Ms[:]
254
            Ms.pop(self.index)
255
            new_eqs.append(EQ(Ms, eq.v))
256
        self.eq = new_eqs.pop(0)
257
        return ProblemData(new_eqs)
258

259
    def compress(self, Ts):
260
        self.compressed_eq = self.eq.compress(Ts)
261
        new_T = matrix(F, [[1]] + [[0] for _ in range(self.pad)])
262
        new_Ts = Ts[:self.index] + [new_T] + Ts[self.index:]
263
        return new_Ts
264

265
    def backward(self, x):
266
        u = self.eq.evaluate(x)
267
        x0 = F_sqrt(u)
268
        v1 = vector(F, [x0] + [0] * self.pad)
269
        x_new = x[:self.index] + [v1] + x[self.index:]
270
        return x_new
271

272
    def compress_backward(self, x_compressed):
273
        u = self.compressed_eq.evaluate(x_compressed)
274
        x0 = F_sqrt(u)
275
        v1 = vector(F, [x0])
276
        x_new = x_compressed[:self.index] + [v1] + x_compressed[self.index:]
277
        return x_new
278

279
def solve_T_even(Ms):
280
    Ms2 = [M + M.transpose() for M in Ms]
281
    n = Ms2[0].nrows()
282
    k = len(Ms2)
283
    if k >= n+1:
284
        raise ValueError
285
    for e in elms:
286
        M0 = Ms2[0] + e * Ms2[1]
287
        if M0.det() != 0:
288
            break
289
    M1 = Ms2[1] * M0.inverse()
290
    M2 = Ms2[2] * M0.inverse()
291
    M3 = Ms2[3] * M0.inverse()
292
    for e1 in elms:
293
        for e2 in elms:
294
            MM = M1 + e1 * M2 + e2 * M3
295
            r = MM.charpoly().roots()
296
            if len(r) == 0:
297
                continue
298
            r0 = r[0][0]
299
            b = (MM-r0*identity_matrix(F, n)).left_kernel().basis()[0]
300
            if b[0] == 0:
301
                continue
302
            b = b / b[0]
303
            T = identity_matrix(F, n)
304
            T[0] = b
305
            Ms3 = [T*M*T.transpose() for M in Ms]
306
            C = zero_matrix(F, k, n-1)
307
            for i in range(k):
308
                M_tmp = Ms3[i] + Ms3[i].transpose()
309
                C[i] = M_tmp[0][1:]
310
            if C.rank() >= n-1:
311
                continue
312
            bs = C.right_kernel().matrix()
313
            T2 = block_matrix(F, [[identity_matrix(1),0], [0,bs]])
314
            return (T2*T).transpose()
315

316
def solve_T_small(Ms):
317
    n = Ms[0].nrows()
318
    k = len(Ms)
319
    if k >= n+1:
320
        raise ValueError
321
    C = zero_matrix(F, k, n-1)
322
    for i in range(k):
323
        M_tmp = Ms[i] + Ms[i].transpose()
324
        C[i] = M_tmp[0][1:]
325
    bs = C.right_kernel().matrix()
326
    T2 = block_matrix(F, [[identity_matrix(1),0], [0,bs]])
327
    return T2.transpose()
328

329
def solve_T_odd(Ms):
330
    Ms2 = [M + M.transpose() for M in Ms]
331
    n = Ms2[0].nrows()
332
    k = len(Ms2)
333
    if k >= n+1:
334
        raise ValueError
335
    M0 = Ms2[0]
336
    M1 = Ms2[1]
337
    M2 = Ms2[2]
338
    for e1 in elms:
339
        for e2 in elms:
340
            MM = M0 + e1 * M1 + e2 * M2
341
            b = MM.left_kernel().basis()[0]
342
            if b[0] == 0:
343
                continue
344
            b = b / b[0]
345
            T = identity_matrix(F, n)
346
            T[0] = b
347
            Ms3 = [T*M*T.transpose() for M in Ms]
348
            C = zero_matrix(F, k, n-1)
349
            for i in range(k):
350
                M_tmp = Ms3[i] + Ms3[i].transpose()
351
                C[i] = M_tmp[0][1:]
352
            if C.rank() >= n-1:
353
                continue
354
            bs = C.right_kernel().matrix()
355
            T2 = block_matrix(F, [[identity_matrix(1),0], [0,bs]])
356
            return (T2*T).transpose()
357

358
def make_one_linear(problem: ProblemData, index):
359
    ns = problem.ns
360
    Ms = problem.get_key(index)
361
    if problem.m <= 2:
362
        T = solve_T_small(Ms)
363
    elif ns[index] % 2 == 0:
364
        T = solve_T_even(Ms)
365
    else:
366
        T = solve_T_odd(Ms)
367
    if T is None:
368
        raise ValueError
369
    transform = Transform1(index, T)
370
    problem = transform.forward(problem)
371
    return problem, transform
372

373
def eliminate(problem: ProblemData, index):
374
    transform = Transform0(index)
375
    problem = transform.forward(problem)
376
    return problem, transform
377

378
def remove_one_key(problem: ProblemData, index):
379
    transform = Transform2(index)
380
    problem = transform.forward(problem)
381
    return problem, transform
382

383
names = ['Miku', 'Ichika', 'Minori', 'Kohane', 'Tsukasa', 'Kanade', 'Mai']
384
pks_uov = [uov.expand_pk(uov.shake256(name.encode(), 43576)) for name in names]
385
pks = [from_uov(pk) for pk in pks_uov]
386
t = uov.shake256(b'shrooms', 44)
387
target = [F.from_integer(i) for i in t]
388

389
# pks = [gen(n, 44) for _ in range(7)]
390
# target = random_vector(F, m)
391
eqs = [EQ(list(Ms), t) for Ms, t in zip(list(zip(*pks)), target)]
392

393
problem_orig = ProblemData(eqs)
394

395
# random shuffle
396
for i in range(256):
397
    ind = random.randrange(1, 44)
398
    elm = random.choice(elms)
399
    eqs[ind] = eqs[ind] + elm * eqs[0]
400
    eqs[0], eqs[ind] = eqs[ind], eqs[0]
401
    if i > 128 and eqs[0].v != 0:
402
        break
403

404
problem0 = ProblemData(eqs[:m])
405
rest_eq0, rest_eq1 = eqs[m:]
406

407
def walk(problem, index):
408
    print(problem.ns, problem.k, problem.m, index)
409
    n = problem.ns[index]
410
    m = problem.m
411
    if n == 0:
412
        raise ValueError
413
    if n < m:
414
        return [remove_one_key(problem, index)]
415
    elif n > m+1:
416
        prob1, transform1 = make_one_linear(problem, index)
417
        prob2, transform2 = eliminate(prob1, index)
418
        return [(prob1, transform1), (prob2, transform2)]
419
    elif n == m+1 or n == m:
420
        prob1, transform1 = make_one_linear(problem, index)
421
        prob2, transform2 = eliminate(prob1, index)
422
        prob3, transform3 = remove_one_key(prob2, index)
423
        return [(prob1, transform1), (prob2, transform2), (prob3, transform3)]
424

425
seq = [6,6,6,5,2,5,5,3,2,1,4,3,0,4,4,4,3,3,2,2,1,1,1,1,1,0,0,0,0,0,0,0,0,0,0,0]
426

427
path = []
428
problem = problem0
429

430
for ind in seq:
431
    path.extend(walk(problem, ind))
432
    problem = path[-1][0]
433

434
print(problem.ns, problem.k, problem.m, ind)
435
eq0 = problem.eqs[0]
436
v = eq0.v
437
M = eq0.Ms[0]
438

439
Ts = [identity_matrix(F, 5)]
440

441
for problem, transform in reversed(path):
442
    Ts = transform.compress(Ts)
443

444
problem0_compressed = problem0.compress(Ts)
445
rest_eq0_compressed = rest_eq0.compress(Ts)
446
rest_eq1_compressed = rest_eq1.compress(Ts)
447

448
def try_solve():
449
    x0 = random_vector(F, 5)
450
    val  = (x0 * M * x0)
451
    if val == 0:
452
        return
453
    c = F_sqrt(v / val)
454
    x0 = c * x0
455
    sol = [x0]
456

457
    for problem, transform in reversed(path):
458
        # assert problem.check(sol), (problem.ns, problem.k, problem.m, len(sol), [len(x) for x in sol], type(transform))
459
        # sol = transform.backward(sol)
460
        sol = transform.compress_backward(sol)
461
    # assert problem0.check(sol)
462
    # assert problem0_compressed.check(sol)
463
    if rest_eq0_compressed.evaluate(sol) == 0:
464
        if rest_eq1_compressed.evaluate(sol) == 0:
465
            print("Found solution")
466
            return sol
467

468
pbar = tqdm()
469
def worker_init():
470
    set_random_seed(os.getpid())
471

472
with ProcessPoolExecutor(max_workers=16, initializer=worker_init) as executor:
473
    while True:
474
        tasks = [executor.submit(try_solve) for _ in range(256)]
475
        for task in tasks:
476
            sol = task.result()
477
            if sol is not None:
478
                break
479
            pbar.update(1)
480
        for task in tasks:
481
            task.cancel()
482
        if sol is not None:
483
            break
484
pbar.close()
485

486
real_sol = [T*x for T, x in zip(Ts, sol)]
487
assert problem0.check(real_sol)
488
assert rest_eq0.evaluate(real_sol) == 0
489
assert rest_eq1.evaluate(real_sol) == 0
490
assert problem_orig.check(real_sol)
491
sol = to_uov(real_sol)
492

493
from pwn import xor
494
t = xor(*[uov.pubmap(s, pk) for s, pk in zip(sol, pks_uov)])
495
print(t.hex())
496
print(uov.shake256(b'shrooms', 44).hex())

More reduction tricks#

The last challenge is to make it to 6 keys. If we just apply the same method, we need to guess 5 equations, which is $2^{40}$ complexity. From Neobeo’s experiment which is written in c, it takes $~1$ hour 20 cores to hit $2^{32}$ times, so hitting $2^{40}$ is not feasible for a two day CTF.

Extending to affine variable#

The idea is that we can reduce $a_i(x_0^2+x_0)+Q_i(x_1,x_2,\cdots,x_n)$ . Here the coefficients of $x_0$ are not $0$ . However, it is still proportional to $x_0^2$ , so the elimination still works.

The equations now becomes $f_{A_i}'(\mathbf{t}_0,\mathbf{t}_0)=kf_{A_i}(\mathbf{t}_0)$ , where $k$ is a constant. After eliminating $k$ , we get $m-1$ equations. Still we can choose good $\mathbf{t}_0$ to remove one equation, but the requirement is more strict now. For the corresponding $(A+A^\top)\mathbf{t}_0=0$ , $f_{A}(\mathbf{t}_0)=0$ is also required. The square terms are hard to control, so we can only guess it with probability $\frac{1}{256}$ .

Another problem is that only half of the equations x^2+x=? have solutions. It seems like a big problem because each time we have half of the possibility to fail. But actually, the rest of the half have two solutions, so in average we still have one solution in every backward step. This is also the main idea in Cheng’s paper².

Therefore, we derive the following reduction:

Lemma. $\text{MQ}^{m}(n_1, \cdots, n_k)$ can be reduced to $\text{MQ}^{m-1}(n_1-m+1, \cdots, n_k)$ .

Homogeneous#

After the reduction, we can solve it with 4 guessed equations. The solution done by Neobeo implemented it and can finish in 1 hour with 20 cores. But it’s still not enough for my pure sage implementation. So I came up with this idea.

You may have noticed that the equations of UOV are homogeneous, so it has a nice property $f(ax)=a^2f(x)$ . My idea is to get one free equation from the homogeneous property.

The homogeneous says that we can sign all $a^2\mathbf{t}$ simultaneously. If $\mathbf{t}=0$ , we get 255 different non-zero signs of $\mathbf{t}=0$ , that makes a free equation possible.

By doing $f_{A_i}-cf_{A_j} = t_i - ct_j$ at the beginning, we can let $t_i=0$ except for one equation. Therefore, only $2^{24}$ attempts are needed to guess the rest 3 equations.

6 keys#

Similar to 7 keys, the script shows the variable usage in the reduction procedure.

1
n = 112
2
m = 40 # plus three guessed equations and one free equation
3
k = 6
4

5
t = [n]*k
6

7
seq = [4,5,5,5,4,4,4,2,3,2,3,3,3,2,2,1,1,1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0]
8

9
for ind in seq:
10
    print(t, m, ind)
11
    if t[ind] < m:
12
        if t[ind] == 0:
13
            raise ValueError
14
        t[ind] = 0
15
        m -= 1
16
        continue
17
    elif t[ind] <= m+1:
18
        t[ind] = 0
19
        m -= 2
20
        continue
21
    else:
22
        if ind == 2 or ind == 4 or m <= 3: # to speed up the forward step
23
            t[ind] -= m
24
            m -= 1
25
        else:
26
            t[ind] -= (m-1)
27
            m -= 1

Final code for 6 keys#

Since x^2+x=? may have 0 or 2 solutions, the backward step should keep a list of solutions. Thus batch_compress_backward is added to the Transform class.

1
from tqdm import tqdm
2
from dataclasses import dataclass
3
from functools import lru_cache
4
from concurrent.futures import ProcessPoolExecutor
5
import random, os
6
from uov import uov_1p_pkc as uov
7

8
F = GF(2**8, name='z', modulus=x^8 + x^4 + x^3 + x + 1)
9
z = F.gen()
10
R = PolynomialRing(F, 'xbar')
11
xbar = R.gen()
12

13
elms = [F.from_integer(i) for i in range(2**8)]
14

15
def gen(n, m):
16
    U = random_matrix(F, n, n)
17
    while U.det() == 0:
18
        U = random_matrix(F, n, n)
19
    pk = []
20
    for _ in range(m):
21
        M = random_matrix(F, n, n)
22
        M[:m, :m] = 0
23
        M = U.transpose() * M * U
24
        pk.append(M)
25
    return pk
26

27
def from_uov_vector(elms):
28
    m = uov.m
29
    bs = elms.to_bytes(m, 'big')
30
    vs = [F.from_integer(i) for i in bs]
31
    return vs
32

33
def from_uov(pk):
34
    v = uov.v
35
    m = uov.m
36
    Ms = [[[0] * n for _ in range(n)] for _ in range(m)]
37
    m1 = uov.unpack_mtri(pk, v)
38
    m2 = uov.unpack_mrect(pk[uov.p1_sz:], v, m)
39
    m3 = uov.unpack_mtri(pk[uov.p1_sz + uov.p2_sz:], m)
40

41
    for i in range(v):
42
        for j in range(i, v):
43
            vec = from_uov_vector(m1[i][j])
44
            for k in range(m):
45
                Ms[k][i][j] = vec[k]
46

47
    for i in range(v):
48
        for j in range(m):
49
            vec = from_uov_vector(m2[i][j])
50
            for k in range(m):
51
                Ms[k][i][j + v] = vec[k]
52

53
    for i in range(m):
54
        for j in range(i, m):
55
            vec = from_uov_vector(m3[i][j])
56
            for k in range(m):
57
                Ms[k][i + v][j + v] = vec[k]
58

59
    Ms = [matrix(F, M) for M in Ms]
60
    return Ms
61

62
def to_uov(sol):
63
    output = []
64
    for s in sol:
65
        output.append([x.to_integer() for x in s.list()])
66
    output = [bytearray(x) for x in output]
67
    return output
68

69
def F_sqrt(x):
70
    return x.square_root()
71

72
quad_sols = {}
73
for u in elms:
74
    for t in elms:
75
        quad_sols[(u, t)] = []
76
    for x in elms:
77
        t = u * x + x**2
78
        quad_sols[(u, t)].append(x)
79

80
def F_quad_solve(u, t):
81
    # solve x^2 + u * x + t = 0
82
    return quad_sols[(u, t)]
83

84
class EQ:
85
    def __init__(self, Ms: list, v):
86
        self.Ms = Ms
87
        self.v = v
88

89
    def evaluate(self, sigs):
90
        result = 0
91
        for M, sig in zip(self.Ms, sigs):
92
            result += sig * M * sig
93
        result += self.v
94
        return result
95

96
    def evaluate_M(self, sigs):
97
        result = 0
98
        for M, sig in zip(self.Ms, sigs):
99
            result += sig * M * sig
100
        return result
101

102
    def compress(self, Ts):
103
        new_Ms = []
104
        for M, T in zip(self.Ms, Ts):
105
            new_Ms.append(T.transpose() * M * T)
106
        return EQ(new_Ms, self.v)
107

108
    @property
109
    def ns(self):
110
        return [M.nrows() for M in self.Ms]
111

112
    def __mul__(self, other):
113
        other = F(other)
114
        new_Ms = [other * M for M in self.Ms]
115
        new_v = other * self.v
116
        return EQ(new_Ms, new_v)
117

118
    def __add__(self, other):
119
        assert isinstance(other, EQ)
120
        assert len(self.Ms) == len(other.Ms)
121
        new_Ms = [M + other.Ms[i] for i, M in enumerate(self.Ms)]
122
        new_v = self.v + other.v
123
        return EQ(new_Ms, new_v)
124

125
    def  __rmul__(self, other):
126
        other = F(other)
127
        return self.__mul__(other)
128

129
class ProblemData:
130
    def __init__(self, eqs: list[EQ]):
131
        self.eqs = eqs
132
        self.ns = eqs[0].ns
133
        self.k = len(self.ns)
134
        self.m = len(eqs)
135
        for i in range(1, self.m):
136
            assert eqs[i].ns == self.ns
137

138
    def check(self, sigs):
139
        vs = [eq.evaluate(sigs) for eq in self.eqs]
140
        return all(v == 0 for v in vs)
141

142
    def get_M(self, i, index):
143
        return self.eqs[i].Ms[index]
144

145
    def get_key(self, index):
146
        Ms = [self.get_M(i, index) for i in range(self.m)]
147
        return Ms
148

149
    def compress(self, Ts):
150
        new_eqs = [eq.compress(Ts) for eq in self.eqs]
151
        return ProblemData(new_eqs)
152

153
class Transform:
154
    def __init__(self):
155
        pass
156

157
    def forward(self, problem: ProblemData):
158
        raise NotImplementedError
159

160
    def compress(self, Ts):
161
        raise NotImplementedError
162

163
    def backward(self, x):
164
        raise NotImplementedError
165

166
    def compress_backward(self, x_compressed):
167
        raise NotImplementedError
168

169
    def batch_compress_backward(self, x_compressed_list):
170
        raise NotImplementedError
171

172
class Transform0(Transform):
173
    def __init__(self, index):
174
        self.index = index
175
        self.eq = None
176

177
    def forward(self, problem: ProblemData):
178
        Ms = problem.get_key(self.index)
179
        Ms = [M + M.transpose() for M in Ms]
180
        assert all(M[0].is_zero() for M in Ms)
181
        m = len(Ms)
182
        eqs = problem.eqs[:]
183
        for i0 in range(m):
184
            if Ms[i0][0, 0] != 0:
185
                break
186
        eqs[0], eqs[i0] = eqs[i0], eqs[0]
187
        coeffs = [eqs[i].Ms[self.index][0, 0] for i in range(problem.m)]
188
        eqs[0] = eqs[0] * coeffs[0].inverse()
189
        for i in range(1, problem.m):
190
            eqs[i] = eqs[i] + (-coeffs[i]) * eqs[0]
191
        assert eqs[0].Ms[self.index][0, 0] == 1
192
        for i in range(problem.m):
193
            eqs[i].Ms[self.index] = eqs[i].Ms[self.index][1:, 1:]
194
        self.eq = eqs.pop(0)
195
        return ProblemData(eqs)
196

197
    def compress(self, Ts):
198
        self.compressed_eq = self.eq.compress(Ts)
199
        new_Ts = Ts[:]
200
        new_Ts[self.index] = block_matrix(F, [[identity_matrix(1), 0], [0, Ts[self.index]]])
201
        return new_Ts
202

203
    def backward(self, x):
204
        u = self.eq.evaluate(x)
205
        x0 = F_sqrt(u)
206
        v0 = x[self.index]
207
        v1 = vector(F, [x0] + v0.list())
208
        x_new = x[:]
209
        x_new[self.index] = v1
210
        return x_new
211

212
    def compress_backward(self, x_compressed):
213
        u = self.compressed_eq.evaluate(x_compressed)
214
        x0 = F_sqrt(u)
215
        v0 = x_compressed[self.index]
216
        v1 = vector(F, [x0] + v0.list())
217
        x_new = x_compressed[:]
218
        x_new[self.index] = v1
219
        return x_new
220

221
    def batch_compress_backward(self, x_compressed_list):
222
        return [self.compress_backward(x) for x in x_compressed_list]
223

224

225
class Transform1(Transform):
226
    def __init__(self, index, T):
227
        self.index = index
228
        self.T = T
229

230
    def forward(self, problem):
231
        return ProblemData([self.apply_T(eq) for eq in problem.eqs])
232

233
    def apply_T(self, eq: EQ):
234
        Ms = eq.Ms[:]
235
        Ms[self.index] = self.T.transpose() * Ms[self.index] * self.T
236
        return EQ(Ms, eq.v)
237

238
    def compress(self, Ts):
239
        new_Ts = Ts[:]
240
        new_Ts[self.index] = self.T * Ts[self.index]
241
        return new_Ts
242

243
    def backward(self, x):
244
        x0 = x[self.index]
245
        x_new = x[:]
246
        x_new[self.index] = self.T * x0
247
        return x_new
248

249
    def compress_backward(self, x_compressed):
250
        return x_compressed
251

252
    def batch_compress_backward(self, x_compressed_list):
253
        return x_compressed_list[:]
254

255
class Transform2(Transform):
256
    def __init__(self, index):
257
        self.index = index
258
        self.eq = None
259
        self.pad = 0
260

261
    def forward(self, problem: ProblemData):
262
        Ms = problem.get_key(self.index)
263
        m = len(Ms)
264
        self.pad = Ms[0].nrows() - 1
265
        eqs = problem.eqs[:]
266
        for i0 in range(m):
267
            if Ms[i0][0, 0] != 0:
268
                break
269
        eqs[0], eqs[i0] = eqs[i0], eqs[0]
270
        coeffs = [eqs[i].Ms[self.index][0, 0] for i in range(problem.m)]
271
        eqs[0] = eqs[0] * coeffs[0].inverse()
272
        for i in range(1, problem.m):
273
            eqs[i] = eqs[i] + (-coeffs[i]) * eqs[0]
274
        assert eqs[0].Ms[self.index][0, 0] == 1
275
        new_eqs = []
276
        for i in range(problem.m):
277
            eq = eqs[i]
278
            Ms = eq.Ms[:]
279
            Ms.pop(self.index)
280
            new_eqs.append(EQ(Ms, eq.v))
281
        self.eq = new_eqs.pop(0)
282
        return ProblemData(new_eqs)
283

284
    def compress(self, Ts):
285
        self.compressed_eq = self.eq.compress(Ts)
286
        new_T = matrix(F, [[1]] + [[0] for _ in range(self.pad)])
287
        new_Ts = Ts[:self.index] + [new_T] + Ts[self.index:]
288
        return new_Ts
289

290
    def backward(self, x):
291
        u = self.eq.evaluate(x)
292
        x0 = F_sqrt(u)
293
        v1 = vector(F, [x0] + [0] * self.pad)
294
        x_new = x[:self.index] + [v1] + x[self.index:]
295
        return x_new
296

297
    def compress_backward(self, x_compressed):
298
        u = self.compressed_eq.evaluate(x_compressed)
299
        x0 = F_sqrt(u)
300
        v1 = vector(F, [x0])
301
        x_new = x_compressed[:self.index] + [v1] + x_compressed[self.index:]
302
        return x_new
303

304
    def batch_compress_backward(self, x_compressed_list):
305
        return [self.compress_backward(x) for x in x_compressed_list]
306

307
class Transform3(Transform):
308
    def __init__(self, index):
309
        self.index = index
310
        self.eq = None
311
        self.u = None
312

313
    def forward(self, problem: ProblemData):
314
        Ms = problem.get_key(self.index)
315
        m = len(Ms)
316
        for i0 in range(m):
317
            if Ms[i0][0, 0] != 0:
318
                break
319
        eqs = problem.eqs[:]
320
        eqs[0], eqs[i0] = eqs[i0], eqs[0]
321
        coeffs = [eqs[i].Ms[self.index][0, 0] for i in range(problem.m)]
322
        eqs[0] = eqs[0] * coeffs[0].inverse()
323
        for i in range(1, problem.m):
324
            eqs[i] = eqs[i] + (-coeffs[i]) * eqs[0]
325
        M_is = [eqs[i].Ms[self.index] for i in range(problem.m)]
326
        for i in range(1, m):
327
            assert M_is[i][0, 0] == 0
328
            assert (M_is[i] + M_is[i].transpose())[0].is_zero()
329
        self.u = (M_is[0]+M_is[0].transpose())[0][1:]
330
        for i in range(problem.m):
331
            eqs[i].Ms[self.index] = eqs[i].Ms[self.index][1:, 1:]
332
        self.eq = eqs.pop(0)
333
        return ProblemData(eqs)
334

335
    def compress(self, Ts):
336
        self.compressed_eq = self.eq.compress(Ts)
337
        self.compressed_u = self.u * Ts[self.index]
338
        new_Ts = Ts[:]
339
        new_Ts[self.index] = block_matrix(F, [[identity_matrix(1), 0], [0, Ts[self.index]]])
340
        return new_Ts
341

342
    def compress_backward_2(self, x_compressed):
343
        v0 = self.compressed_eq.evaluate(x_compressed)
344
        u0 = self.compressed_u * x_compressed[self.index]
345
        x0_sols = F_quad_solve(u0, v0)
346
        sols = []
347
        for x0 in x0_sols:
348
            v1 = vector(F, [x0] + x_compressed[self.index].list())
349
            x_new = x_compressed[:]
350
            x_new[self.index] = v1
351
            sols.append(x_new)
352
        return sols
353

354
    def batch_compress_backward(self, x_compressed_list):
355
        ret = []
356
        for x in x_compressed_list:
357
            ret.extend(self.compress_backward_2(x))
358
        return ret
359

360
def solve_T_even(Ms):
361
    Ms2 = [M + M.transpose() for M in Ms]
362
    n = Ms2[0].nrows()
363
    k = len(Ms2)
364
    if k >= n+1:
365
        raise ValueError
366
    for e in elms:
367
        M0 = Ms2[0] + e * Ms2[1]
368
        if M0.det() != 0:
369
            break
370
    M1 = Ms2[1] * M0.inverse()
371
    M2 = Ms2[2] * M0.inverse()
372
    M3 = Ms2[3] * M0.inverse()
373
    for e1 in elms:
374
        for e2 in elms:
375
            MM = M1 + e1 * M2 + e2 * M3
376
            r = MM.charpoly().roots()
377
            if len(r) == 0:
378
                continue
379
            r0 = r[0][0]
380
            b = (MM-r0*identity_matrix(F, n)).left_kernel().basis()[0]
381
            if b[0] == 0:
382
                continue
383
            b = b / b[0]
384
            T = identity_matrix(F, n)
385
            T[0] = b
386
            Ms3 = [T*M*T.transpose() for M in Ms]
387
            C = zero_matrix(F, k, n-1)
388
            for i in range(k):
389
                M_tmp = Ms3[i] + Ms3[i].transpose()
390
                C[i] = M_tmp[0][1:]
391
            if C.rank() >= n-1:
392
                continue
393
            bs = C.right_kernel().matrix()
394
            T2 = block_matrix(F, [[identity_matrix(1),0], [0,bs]])
395
            return (T2*T).transpose()
396

397
def solve_T_small(Ms):
398
    n = Ms[0].nrows()
399
    k = len(Ms)
400
    if k >= n+1:
401
        raise ValueError
402

403
    C = zero_matrix(F, k, n-1)
404
    for i in range(k):
405
        M_tmp = Ms[i] + Ms[i].transpose()
406
        C[i] = M_tmp[0][1:]
407
    bs = C.right_kernel().matrix()
408
    T2 = block_matrix(F, [[identity_matrix(1),0], [0,bs]])
409
    return T2.transpose()
410

411
def solve_T_odd(Ms):
412
    Ms2 = [M + M.transpose() for M in Ms]
413
    n = Ms2[0].nrows()
414
    k = len(Ms2)
415
    if k >= n+1:
416
        raise ValueError
417
    M0 = Ms2[0]
418
    M1 = Ms2[1]
419
    M2 = Ms2[2]
420
    for e1 in elms:
421
        for e2 in elms:
422
            MM = M0 + e1 * M1 + e2 * M2
423
            b = MM.left_kernel().basis()[0]
424
            if b[0] == 0:
425
                continue
426
            b = b / b[0]
427
            T = identity_matrix(F, n)
428
            T[0] = b
429
            Ms3 = [T*M*T.transpose() for M in Ms]
430
            C = zero_matrix(F, k, n-1)
431
            for i in range(k):
432
                M_tmp = Ms3[i] + Ms3[i].transpose()
433
                C[i] = M_tmp[0][1:]
434
            if C.rank() >= n-1:
435
                continue
436
            bs = C.right_kernel().matrix()
437
            T2 = block_matrix(F, [[identity_matrix(1),0], [0,bs]])
438
            return (T2*T).transpose()
439

440
def solve_T_even_semi_linear(Ms):
441
    Ms2 = [M + M.transpose() for M in Ms]
442
    n = Ms2[0].nrows()
443
    k = len(Ms2)
444
    if k >= n+1:
445
        raise ValueError
446
    for e0 in elms:
447
        M0 = Ms2[0] + e0 * Ms2[1]
448
        if M0.det() != 0:
449
            break
450
    M1 = Ms2[1] * M0.inverse()
451
    M2 = Ms2[2] * M0.inverse()
452
    M3 = Ms2[3] * M0.inverse()
453
    for e1 in elms:
454
        for e2 in elms:
455
            MM = M1 + e1 * M2 + e2 * M3
456
            r = MM.charpoly().roots()
457
            if len(r) == 0:
458
                continue
459
            r0 = r[0][0]
460
            U0 = r0*(Ms[0]+e0*Ms[1])+Ms[1]+e1 * Ms[2]+e2 * Ms[3]
461
            for b in (MM-r0*identity_matrix(F, n)).left_kernel().basis():
462
                if b[0] == 0:
463
                    continue
464
                if b*U0*b != 0:
465
                    continue
466
                b = b / b[0]
467
                T = identity_matrix(F, n)
468
                T[0] = b
469

470
                Ms3 = [T*M*T.transpose() for M in Ms]
471
                Ms3_copy = Ms3[:]
472
                for i0 in range(k):
473
                    if Ms3_copy[i0][0, 0] != 0:
474
                        break
475
                Ms3_copy[0], Ms3_copy[i0] = Ms3_copy[i0], Ms3_copy[0]
476
                Ms3_copy[0] = Ms3_copy[0] * Ms3_copy[0][0, 0].inverse()
477
                for i in range(1, k):
478
                    if Ms3_copy[i][0, 0] == 0:
479
                        continue
480
                    Ms3_copy[i] = Ms3_copy[i] - Ms3_copy[i][0, 0] * Ms3_copy[0]
481

482
                C = zero_matrix(F, k-1, n-1)
483
                for i in range(1, k):
484
                    M_tmp = Ms3_copy[i] + Ms3_copy[i].transpose()
485
                    C[i-1] = M_tmp[0][1:]
486
                if C.rank() >= n-1:
487
                    continue
488
                bs = C.right_kernel().matrix()
489
                T2 = block_matrix(F, [[identity_matrix(1),0], [0,bs]])
490
                return (T2*T).transpose()
491

492
def solve_T_odd_semi_linear(Ms):
493
    Ms2 = [M + M.transpose() for M in Ms]
494
    n = Ms2[0].nrows()
495
    k = len(Ms2)
496
    if k >= n:
497
        raise ValueError
498
    M0 = Ms2[0]
499
    M1 = Ms2[1]
500
    M2 = Ms2[2]
501
    for e1 in elms:
502
        for e2 in elms:
503
            MM = M0 + e1 * M1 + e2 * M2
504
            b = MM.left_kernel().basis()[0]
505
            if b[0] == 0:
506
                continue
507
            b = b / b[0]
508
            T = identity_matrix(F, n)
509
            T[0] = b
510

511
            U0 = Ms[0] + e1 * Ms[1] + e2 * Ms[2]
512
            if b*U0*b != 0:
513
                continue
514

515
            Ms3 = [T*M*T.transpose() for M in Ms]
516
            Ms3_copy = Ms3[:]
517
            for i0 in range(k):
518
                if Ms3_copy[i0][0, 0] != 0:
519
                    break
520
            Ms3_copy[0], Ms3_copy[i0] = Ms3_copy[i0], Ms3_copy[0]
521
            Ms3_copy[0] = Ms3_copy[0] * Ms3_copy[0][0, 0].inverse()
522
            for i in range(1, k):
523
                if Ms3_copy[i][0, 0] == 0:
524
                    continue
525
                Ms3_copy[i] = Ms3_copy[i] - Ms3_copy[i][0, 0] * Ms3_copy[0]
526

527
            C = zero_matrix(F, k-1, n-1)
528
            for i in range(1, k):
529
                M_tmp = Ms3_copy[i] + Ms3_copy[i].transpose()
530
                C[i-1] = M_tmp[0][1:]
531
            if C.rank() >= n-1:
532
                continue
533
            bs = C.right_kernel().matrix()
534
            T2 = block_matrix(F, [[identity_matrix(1),0], [0,bs]])
535
            return (T2*T).transpose()
536

537
def solve_T_small_semi_linear(Ms):
538
    Ms2 = [M + M.transpose() for M in Ms]
539
    n = Ms2[0].nrows()
540
    k = len(Ms2)
541

542
    Ms3_copy = Ms[:]
543
    for i0 in range(k):
544
        if Ms3_copy[i0][0, 0] != 0:
545
            break
546
    Ms3_copy[0], Ms3_copy[i0] = Ms3_copy[i0], Ms3_copy[0]
547
    Ms3_copy[0] = Ms3_copy[0] * Ms3_copy[0][0, 0].inverse()
548
    for i in range(1, k):
549
        if Ms3_copy[i][0, 0] == 0:
550
            continue
551
        Ms3_copy[i] = Ms3_copy[i] - Ms3_copy[i][0, 0] * Ms3_copy[0]
552

553
    C = zero_matrix(F, k-1, n-1)
554
    for i in range(1, k):
555
        M_tmp = Ms3_copy[i] + Ms3_copy[i].transpose()
556
        C[i-1] = M_tmp[0][1:]
557
    bs = C.right_kernel().matrix()
558
    T2 = block_matrix(F, [[identity_matrix(1),0], [0,bs]])
559
    return T2.transpose()
560

561
def make_one_linear(problem: ProblemData, index):
562
    ns = problem.ns
563
    Ms = problem.get_key(index)
564
    if problem.m <= 2:
565
        T = solve_T_small(Ms)
566
    elif ns[index] % 2 == 0:
567
        T = solve_T_even(Ms)
568
    else:
569
        T = solve_T_odd(Ms)
570
    if T is None:
571
        raise ValueError
572
    transform = Transform1(index, T)
573
    problem = transform.forward(problem)
574
    return problem, transform
575

576
def make_one_semi_linear(problem: ProblemData, index):
577
    ns = problem.ns
578
    Ms = problem.get_key(index)
579
    if problem.m <= 3:
580
        T = solve_T_small_semi_linear(Ms)
581
    elif ns[index] % 2 == 0:
582
        T = solve_T_even_semi_linear(Ms)
583
    else:
584
        T = solve_T_odd_semi_linear(Ms)
585
    if T is None:
586
        raise ValueError
587
    transform = Transform1(index, T)
588
    problem = transform.forward(problem)
589
    return problem, transform
590

591
def eliminate(problem: ProblemData, index):
592
    transform = Transform0(index)
593
    problem = transform.forward(problem)
594
    return problem, transform
595

596
def eliminate_semi_linear(problem: ProblemData, index):
597
    transform = Transform3(index)
598
    problem = transform.forward(problem)
599
    return problem, transform
600

601
def remove_one_key(problem: ProblemData, index):
602
    transform = Transform2(index)
603
    problem = transform.forward(problem)
604
    return problem, transform
605

606
n = 112
607
m = 40
608

609
names = ['Miku', 'Ichika', 'Minori', 'Kohane', 'Tsukasa', 'Kanade']
610
pks_uov = [uov.expand_pk(uov.shake256(name.encode(), 43576)) for name in names]
611
pks = [from_uov(pk) for pk in pks_uov]
612
message = b'SEKAI'
613
t = uov.shake256(message, 44)
614
target = [F.from_integer(i) for i in t]
615

616
eqs = [EQ(list(Ms), t) for Ms, t in zip(list(zip(*pks)), target)]
617

618
problem_orig = ProblemData(eqs[:])
619

620
# random shuffle
621
for i in range(256):
622
    ind = random.randrange(1, len(eqs))
623
    elm = random.choice(elms)
624
    eqs[ind] = eqs[ind] + elm * eqs[0]
625
    eqs[0], eqs[ind] = eqs[ind], eqs[0]
626
    if i > 128 and eqs[0].v != 0:
627
        break
628

629
last_eq = eqs.pop(0)
630
last_eq = last_eq * last_eq.v.inverse()
631
for i in range(len(eqs)):
632
    eqs[i] = eqs[i] + (-eqs[i].v) * last_eq
633

634
problem0 = ProblemData(eqs[:m])
635
rest_eqs = eqs[m:]
636

637
def walk(problem, index):
638
    print(problem.ns, problem.k, problem.m, index)
639
    n = problem.ns[index]
640
    m = problem.m
641
    if n == 0:
642
        raise ValueError
643
    if n < m:
644
        return [remove_one_key(problem, index)]
645
    elif n > m+1:
646
        if index != 2 and index != 4:
647
            prob1, transform1 = make_one_semi_linear(problem, index)
648
            prob2, transform2 = eliminate_semi_linear(prob1, index)
649
        else:
650
            prob1, transform1 = make_one_linear(problem, index)
651
            prob2, transform2 = eliminate(prob1, index)
652
        return [(prob1, transform1), (prob2, transform2)]
653
    elif n == m+1 or n == m:
654
        prob1, transform1 = make_one_linear(problem, index)
655
        prob2, transform2 = eliminate(prob1, index)
656
        prob3, transform3 = remove_one_key(prob2, index)
657
        return [(prob1, transform1), (prob2, transform2), (prob3, transform3)]
658

659
seq = [4,5,5,5,4,4,4,2,3,2,3,3,3,2,2,1,1,1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0]
660

661
path = []
662
problem = problem0
663

664
for ind in seq:
665
    path.extend(walk(problem, ind))
666
    problem = path[-1][0]
667

668
print(problem.ns, problem.k, problem.m, ind)
669
eq0 = problem.eqs[0]
670
assert eq0.v == 0
671
M = eq0.Ms[0]
672
M = M / M[0][0]
673
M0 = M[1:, 1:]
674
u0 = (M+M.transpose())[0][1:]
675

676
Ts = [identity_matrix(F, 5)]
677

678
for problem, transform in reversed(path):
679
    Ts = transform.compress(Ts)
680

681
problem0_compressed = problem0.compress(Ts)
682
rest_eqs_compressed = [eq.compress(Ts) for eq in rest_eqs]
683
last_eq_compressed = last_eq.compress(Ts)
684

685
def try_solve():
686
    x0 = random_vector(F, 4)
687
    val = (x0 * M0 * x0)
688
    u = u0 * x0
689
    xsols = F_quad_solve(u, val)
690
    if len(xsols) == 0:
691
        return 0, None
692

693
    sols = []
694
    for i in range(len(xsols)):
695
        x = vector(F, [xsols[i]] + x0.list())
696
        sols.append([x])
697

698
    for problem, transform in reversed(path):
699
        if len(sols) == 0:
700
            return 0, None
701
        # assert problem.check(sol), (problem.ns, problem.k, problem.m, len(sol), [len(x) for x in sol], type(transform))
702
        # sol = transform.backward(sol)
703
        sols = transform.batch_compress_backward(sols)
704
    # assert problem0.check(sol)
705
    # assert problem0_compressed.check(sol)
706
    for sol in sols:
707
        find = True
708
        for i, eq in enumerate(rest_eqs_compressed):
709
            if eq.evaluate(sol) != 0:
710
                find = False
711
        if not find:
712
            continue
713
        r = last_eq_compressed.evaluate_M(sol)
714
        if r == 0:
715
            print("Sad...")
716
            continue
717
        v = last_eq_compressed.v
718
        scale = F_sqrt(v/r)
719
        sol = [scale * x for x in sol]
720
        print("Found solution")
721
        return len(sols), sol
722
    return len(sols), None
723

724
pbar = tqdm()
725
def worker_init():
726
    set_random_seed(os.getpid()+random.randint(0, 2**24))
727

728
with ProcessPoolExecutor(max_workers=12, initializer=worker_init) as executor:
729
    tasks = []
730
    for i in range(512):
731
        tasks.append(executor.submit(try_solve))
732

733
    while True:
734
        task = tasks.pop(0)
735
        cnt, sol = task.result()
736
        if sol is not None:
737
            break
738
        pbar.update(cnt)
739
        tasks.append(executor.submit(try_solve))
740
    for task in tasks:
741
        task.cancel()
742
pbar.close()
743

744
real_sol = [T*x for T, x in zip(Ts, sol)]
745
assert problem0.check(real_sol)
746
assert all(eq.evaluate(real_sol) == 0 for eq in rest_eqs)
747
assert problem_orig.check(real_sol)
748

749
sol = to_uov(real_sol)
750
print("".join(s.hex() for s in sol))
751

752
from pwn import xor
753
t = xor(*[uov.pubmap(s, pk) for s, pk in zip(sol, pks_uov)])
754
print(t.hex())
755
print(uov.shake256(message, 44).hex())

Conclusion#

My rating for the (not so) recent hard crypto is unfairy-ring > zerodaycrypto > lance-hard > flcg. The solution is a combination of many tricks and optimizations, which makes it quite complicated and hard to really solve it during a two-day CTF. But the key idea is still close to Cheng’s paper and can’t be further improved to 5 keys.

Neobeo has another independent code implementation written in C++, which is more efficient than my sage implementation. It’s very cool because I have no idea how to do GF(256) and these matrix operations in C++.

Finally, let’s return fairy-ring and break it with no reused keys.

fairy-ring sol

H. Miura, Y. Hashimoto, T. Takagi, Extended algorithm for solving underdefined multivariate quadratic equations ↩
C.M. Cheng, Y. Hashimoto, H. Miura, T. Takagi, A polynomial-time algorithm for solving a class of underdetermined multivariate quadratic equations over fields of odd characteristics ↩ ↩² ↩³
H. Furue, S. Nakamura, T. Takagi, Improving Thomae-Wolf algorithm for solving underdetermined multivariate quadratic polynomial problem. ↩
Y. Hashimoto, An improvement of algorithms to solve under-defined systems of multivariate quadratic equations. ↩