SECCON CTF 2023 Quals writeup - チョコラスクのブログ

SECCON CTF 2023 QualsにチームBunkyoWesternsで参加しました。

国内2位で国内決勝に進出することができました。

今年はCryptoがとても難しかったですが、10solvesの問題を含む4問を解くことができて良かったです。

[crypto] plai_n_rsa
[crypto] RSA 4.0
[crypto] CIGISICGICGICG
[crypto] mystic_harmony

[crypto] plai_n_rsa

次のソースコードとその実行結果が与えられます。

import os

from Crypto.Util.number import bytes_to_long, getPrime

flag = os.getenvb(b"FLAG", b"SECCON{THIS_IS_FAKE}")
assert flag.startswith(b"SECCON{")
m = bytes_to_long(flag)
e = 0x10001
p = getPrime(1024)
q = getPrime(1024)
n = p * q
e = 65537
phi = (p-1)*(q-1)
d = pow(e, -1, phi)
hint = p+q
c = pow(m,e,n)

print(f"e={e}")
print(f"d={d}")
print(f"hint={hint}")
print(f"c={c}")

RSA暗号で、 $n$ が与えられない代わりに、秘密鍵 $d$ および $hint=p+q$ の値が与えられます。

$d=e^{-1}\bmod{\phi(n)}$ なので、

$\displaystyle{ ed=1+k\phi(n) }$

を満たす正整数 $k$ があります。ここで $d\lt \phi(n)$ なので、 $k\lt e$ が成り立ちます。

よって、 $k$ を全探索することができ、 $\phi(n)$ の候補が得られます。さらに $\phi(n)=(p-1)(q-1)=n-hint+1$ より、 $n$ の候補が得られます。

各 $n$ の候補について、 $m=c^{d}\bmod{n}$ により復号を試し、結果が SECCON から始まるものがフラグです。

from Crypto.Util.number import *
from output import *

for k in range(1, e):
    if (e*d-1)%k!=0:
        continue
    phi = (e*d-1)//k
    n = phi+hint-1
    m = pow(c, d, n)
    flag = long_to_bytes(m)
    if b'SECCON' in flag:
        print(flag)
# b'SECCON{thank_you_for_finding_my_n!!!_GOOD_LUCK_IN_SECCON_CTF}'

[crypto] RSA 4.0

次のソースコードとその実行結果が与えられます。

import os

from Crypto.Util.number import bytes_to_long, getStrongPrime

m = bytes_to_long(os.getenvb(b"FLAG", b"FAKEFLAG{THIS_IS_FAKE}"))
e = 0x10001
p = getStrongPrime(1024, e=e)
q = getStrongPrime(1024, e=e)
n = p * q
assert m < n
Q = QuaternionAlgebra(Zmod(n), -1, -1)
i, j, k = Q.gens()
enc = (
    1 * m
    + (3 * m + 1 * p + 337 * q) * i
    + (3 * m + 13 * p + 37 * q) * j
    + (7 * m + 133 * p + 7 * q) * k
) ** e
print(f"{n = }")
print(f"{e = }")
print(f"{enc = }")

関係式 $i^{2}=-1$ , $j^{2}=-1$ , $k=ij=-ji$ で定まる $\mathbb{Z}/n\mathbb{Z}$ 上の四元数環 $Q$ における、RSA暗号の問題です。ここで、平文 $a+bi+cj+dk$ は、フラグ $m$ と $p$ , $q$ の1次式になっています。

$Q$ は、

$\displaystyle{ a+bi+cj+dk \mapsto \begin{pmatrix} a+bi & c+di \\ -c+di & a-bi \end{pmatrix} = M }$

により行列環に埋め込むことができます (Wikipedia)。

この行列の対角化を考えます。

$\displaystyle{ M = P \begin{pmatrix} \alpha & 0 \\ 0 & \beta \end{pmatrix} P^{-1} }$

ここで、 $\alpha, \beta$ は $M$ の固有多項式

$\displaystyle{ \det\begin{pmatrix} \lambda-(a+bi) & -(c+di) \\ -(-c+di) & \lambda-(a-bi) \end{pmatrix} = \lambda^2-2a\lambda+(a^2+b^2+c^2+d^2) }$

の根であり、

$\displaystyle{ P=\begin{pmatrix} r & s \\ t & u \end{pmatrix}=\begin{pmatrix} c+di & c+di \\ \alpha-(a+bi) & \beta-(a+bi) \end{pmatrix} }$

です。

このとき、

$\displaystyle{ \begin{aligned} M^e &= P \begin{pmatrix} \alpha^e & 0 \\ 0 & \beta^e \end{pmatrix} P^{-1} \\\\ &= \frac{1}{ru-st} \begin{pmatrix} ru\alpha^e-st\beta^e & -rs(\alpha^e-\beta^e) \\ tu(\alpha^e-\beta^e) & ru\beta^e-st\alpha^e \end{pmatrix} \end{aligned} }$

です。

よって、 $(a+bi+cj+dk)^{e}=a'+b'i+c'j+d'k$ とおくと

$\displaystyle{ (c'+d'i)tu=-(-c'+d'i)rs }$

です。 $rs=(c+di)^{2}$ , $tu=c^{2}+d^{2}$ を用いて整理すると

$\displaystyle{ cd'=c'd }$

がわかります。対称性より、 $bc'=b'c$ もわかります。

$b=3m+p+337q$ などを用いると、

$\displaystyle{ \begin{aligned} c_{11}m+c_{12}p+c_{13}q &= 0 \pmod{n} \\ c_{21}m+c_{22}p+c_{23}q &= 0 \pmod{n} \end{aligned} }$

のような2本の関係式が得られます。 $m$ を消去すると、

$\displaystyle{ c_{1}p+c_{2}q=0 \pmod{n} }$

のような式が得られます。 $c_{2}$ は $p$ の倍数なので、

$\displaystyle{ p=\gcd(c_{2}, n) }$

により $p$ が求まり、 $n$ を素因数分解できます。

$\alpha, \beta$ は $\mathbb{Z}/n\mathbb{Z}$ 係数の2次方程式の解なので、 $\mathbb{F}_{p^{2}}\times \mathbb{F}_{q^{2}}$ の元とみなせます。さらに、 $\alpha, \beta\not\in \mathbb{F}_{p}$ の場合、 $\alpha, \beta$ は $\mathbb{F}_{p}$ 上共役です。よって、四元数環の元の位数は $(p^{2}-1)(q^{2}-1)$ の約数なので、 $enc$ を $e^{-1}\bmod (p^{2}-1)(q^{2}-1)$ 乗することで復号できます。

n = 2487202132 ...
e = 65537
a = 3859728242 ...
b = 6689400988 ...
c = 1425535224 ...
d = 1693964190 ...

from Crypto.Util.number import *

c11, c12, c13 = 3*c-3*b, c-13*b, 337*c-37*b
c11, c12, c13 = c11%n, c12%n, c13%n
c21, c22, c23 = 3*d-7*b, d-133*b, 337*d-7*b
c21, c22, c23 = c21%n, c22%n, c23%n
c1, c2 = c21*c12-c11*c22, c21*c13-c11*c23
c1, c2 = c1%n, c2%n
p = GCD(c2, n)
q = n//p

Q = QuaternionAlgebra(Zmod(n), -1, -1)
i, j, k = Q.gens()
enc = a+b*i+c*j+d*k
d = inverse(e, (p**2-1)*(q**2-1))
m = enc**d
print(long_to_bytes(int(m[0])))
# b'SECCON{pr0m153_m3!d0_n07_3ncryp7_p_0r_q_3v3n_w17h_r54_4.0}'

実験が下手でなかなか規則性に気づけなかったのですが、行列表示して対角化を考えればとりあえず成分の関係式が2本出てくるはずだと思い、頑張って計算してみることにしました。めちゃめちゃ単純な関係式が導かれてびっくりしましたが、何とか解けてよかったです。

[crypto] CIGISICGICGICG

次のソースコードとその実行結果が与えられます。

import os
from functools import reduce
from secrets import randbelow


flag = os.getenvb(b"FLAG", b"FAKEFLAG{THIS_IS_FAKE}")
p1 = 21267647932558653966460912964485513283
a1 = 6701852062049119913950006634400761786
b1 = 19775891958934432784881327048059215186
p2 = 21267647932558653966460912964485513289
a2 = 10720524649888207044145162345477779939
b2 = 19322437691046737175347391539401674191
p3 = 21267647932558653966460912964485513327
a3 = 8837701396379888544152794707609074012
b3 = 10502852884703606118029748810384117800


def prod(x: list[int]) -> int:
    return reduce(lambda a, b: a * b, x, 1)


def xor(x: bytes, y: bytes) -> bytes:
    return bytes([xi ^ yi for xi, yi in zip(x, y)])


class ICG:
    def __init__(self, p: int, a: int, b: int) -> None:
        self.p = p
        self.a = a
        self.b = b
        self.x = randbelow(self.p)

    def _next(self) -> int:
        if self.x == 0:
            self.x = self.b
            return self.x
        else:
            self.x = (self.a * pow(self.x, -1, self.p) + self.b) % self.p
            return self.x


class CIG:
    L = 256

    def __init__(self, icgs: list[ICG]) -> None:
        self.icgs = icgs
        self.T = prod([icg.p for icg in self.icgs])
        self.Ts = [self.T // icg.p for icg in self.icgs]

    def _next(self) -> int:
        ret = 0
        for icg, t in zip(self.icgs, self.Ts):
            ret += icg._next() * t
            ret %= self.T
        return ret % 2**self.L

    def randbytes(self, n: int) -> bytes:
        ret = b""
        block_size = self.L // 8
        while n > 0:
            ret += self._next().to_bytes(block_size, "big")[: min(n, block_size)]
            n -= block_size
        return ret


if __name__ == "__main__":
    random = CIG([ICG(p1, a1, b1), ICG(p2, a2, b2), ICG(p3, a3, b3)])
    enc_flag = xor(flag, random.randbytes(len(flag)))
    leaked = random.randbytes(300)
    print(f"{enc_flag = }")
    print(f"{leaked = }")

128bit程度の素数 $p_i$ と整数 $a_i$ , $b_i$ ( $i=1,2,3$ ) を用いて、次の漸化式で乱数生成を行っています。

$\displaystyle{ x_{j+1,i}=\frac{a_i}{x_{j,i}}+b_i \pmod{p_i} }$

$j$ 個目の出力は

$\displaystyle{ x_{j,1}p_2p_3+x_{j,2}p_3p_1+x_{j,3}p_1p_2\bmod{p_1p_2p_3} }$

の下256bitとなっています。300bytes (9個) の出力からseedを復元できれば、フラグが求まります。

中国剰余定理により、 $n=p_1p_2p_3$ として $\bmod{n}$ で考えることができます。

乱数生成式は1次分数変換なので、seedを $r$ とすると、 $i$ 個目の出力は $\frac{a_ir+b_i}{c_ir+d_i}\bmod{n}$ の下256bit、と表せます。よって、 $r_i$ ( $i=1,2,\ldots,9$ ) を出力、 $x_i$ を128bit程度の未知数として、

$\displaystyle{ \frac{a_ir+b_i}{c_ir+d_i}=2^{256}x_i+r_i \pmod{n} }$

と書けます。変形すると、

$\displaystyle{ r=\frac{b_i'+d_i'x_i}{a_i'+c_i'x_i} \pmod{n} }$

のようになります。

ここで $x_i$ に関する式と $x_j$ に関する式から $r$ を消去すると、

$\displaystyle{ a_{ij}+b_{ij}x_i+c_{ij}x_j=x_ix_j \pmod{n} }$

のようになります。 $x_ix_j$ は256bit程度で $n$ より十分小さいので、次のような行列

$\displaystyle{ \begin{pmatrix} a_{12} & a_{13} & a_{14} & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0\\ b_{12} & 0 & 0 & a_{23} & a_{24} & 0 & 0 & 1 & 0 & 0 & 0\\ 0 & b_{13} & 0 & b_{23} & 0 & a_{34} & 0 & 0 & 1 & 0 & 0\\ 0 & 0 & b_{14} & 0 & b_{24} & b_{34} & 0 & 0 & 0 & 1 & 0\\ n & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\ 0 & n & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\ 0 & 0 & n & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & n & 0 & 0 & 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & n & 0 & 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & 0 & n & 0 & 0 & 0 & 0 & 0\\ c_{12} & c_{13} & c_{14} & c_{23} & c_{24} & c_{34} & 0 & 0 & 0 & 0 & \infty \end{pmatrix} }$

に対しLLLアルゴリズムを適用すると、最後の行に次のようなベクトル

$\displaystyle{ (x_1x_2, x_1x_3, x_1x_4, x_2x_3, x_2x_4, x_3x_4, x_1, x_2, x_3, x_4, \infty) }$

が現れることが期待できます (実際は $1\leq i\lt j \leq 9$ すべてを使って行列を作ります)。実際は $x_i$ の列に $2^{128}$ 倍程度の重みを付ける (大きさを $x_ix_j$ の列と合わせる) ことで、欲しいベクトルが得られます。

from Crypto.Util.number import *
from output import *

p1 = 21267647932558653966460912964485513283
a1 = 6701852062049119913950006634400761786
b1 = 19775891958934432784881327048059215186
p2 = 21267647932558653966460912964485513289
a2 = 10720524649888207044145162345477779939
b2 = 19322437691046737175347391539401674191
p3 = 21267647932558653966460912964485513327
a3 = 8837701396379888544152794707609074012
b3 = 10502852884703606118029748810384117800
a = [a1,a2,a3]
b = [b1,b2,b3]
p = [p1,p2,p3]
n = p1*p2*p3

r = []
for i in range(9):
    l = leaked[i*32:(i+1)*32]
    l = bytes_to_long(l)
    r.append(l)

mats = []
for i in range(3):
    mat = []
    m = Matrix.identity(GF(p[i]), 2)
    for j in range(9):
        mat.append(m)
        m *= Matrix(GF(p[i]), [[b[i],a[i]],[1,0]])
    mats.append(mat)

mat = []
for i in range(9):
    mati = [[0,0],[0,0]]
    for i1 in range(2):
        for i2 in range(2):
            if i1==0:
                c = CRT_list([int(mats[j][i][i1,i2])*(n//p[j])%p[j] for j in range(3)], p)
            else:
                c = CRT_list([int(mats[j][i][i1,i2]) for j in range(3)], p)
            mati[i1][i2] = c
    mat.append(mati)

for i in range(9):
    mat[i][0][0] = (mat[i][0][0]-r[i]*mat[i][1][0])%n
    mat[i][0][1] = (mat[i][0][1]-r[i]*mat[i][1][1])%n

M = [[0]*46 for _ in range(46)]
k = 0
for i in range(9):
    ai,bi,ci,di = mat[i][0][0],mat[i][0][1],mat[i][1][0],mat[i][1][1]
    for j in range(i):
        aj,bj,cj,dj = mat[j][0][0],mat[j][0][1],mat[j][1][0],mat[j][1][1]
        v1 = (ai*bj-aj*bi)%n
        v2 = (2**256)*(di*aj-ci*bj)%n
        v3 = (2**256)*(cj*bi-ai*dj)%n
        v4 = (2**512)*(ci*dj-cj*di)%n
        v1 = v1*inverse(v4,n)%n
        v2 = v2*inverse(v4,n)%n
        v3 = v3*inverse(v4,n)%n
        M[45][k] = -v1
        M[i][k] = -v2
        M[j][k] = -v3
        M[9+k][k] = n
        k += 1
weight = 2**100
for i in range(9):
    M[i][36+i] = weight
M[45][45] = 2**384

res = Matrix(M).LLL()

xs = [x//weight for x in res[-1][36:45]]
a0,b0,c0,d0 = mat[0][0][0],mat[0][0][1],mat[0][1][0],mat[0][1][1]
r = (-b0+(2**256)*xs[0]*d0)*inverse(int(a0-(2**256)*xs[0]*c0), n)%n
rs = [r%pi for pi in p]

def xor(a, b):
    return bytes([x^^y for x,y in zip(a, b)])

flag = b''
for i in range(2,-1,-1):
    for j in range(3):
        rs[j] = inverse(int(rs[j]-b[j]),p[j])*a[j]%p[j]
    ret = 0
    for j in range(3):
        ret += rs[j]*(n//p[j])
    ret %= n
    ret %= 2**256
    m = enc_flag[32*i:min(len(enc_flag),32*(i+1))]
    flag = xor(m,long_to_bytes(ret)) + flag
print(flag)
# b'SECCON{ICG1c6iC6icgic6icgcIgIcg1C6ic6ICGICG1cGicG1C61CG1cG1c61cgIcg}'

[crypto] mystic_harmony

次のソースコードが与えられます。

import random
import Crypto.Cipher.AES as AES
from Crypto.Util.number import long_to_bytes
import hashlib
from flag import FLAG

R.<x, y> = PolynomialRing(GF(2))
size = 2^8
K.<alpha> = GF(size, modulus=x^8+x^4+x^3+x^2+1)

def make_human_world(human_world_size):
    H = 0
    for i in range(human_world_size):
        for j in range(human_world_size):
            H += (x^i) * (y^j) * (alpha^(random.randint(0,size-2)))
    return H

def make_spirit_world(H, spirit_world_size_param):
    Gx = prod(x-alpha^i for i in range(1,spirit_world_size_param+1))
    Gy = prod(y-alpha^i for i in range(1,spirit_world_size_param+1))
    return H % (Gx + Gy)

def make_disharmony(C, count):
    x_set = set()

    D = 0
    for i in range(count):
        r = random.randint(0,size-2)
        p = random.choice(list(C.dict().keys()))
        while p[0] in x_set:
            p = random.choice(list(C.dict().keys()))
        x_set.add(p[0])
        D += (x^p[0]) * (y^p[1]) * alpha^(r)
    return D

def make_key(D):
    key_seed = b""
    for pos, value in sorted(list(D.dict().items())):
        x = pos[0]
        y = pos[1]
        power = discrete_log(value, alpha, size-1)
        key_seed += long_to_bytes(x) + long_to_bytes(y) + long_to_bytes(power)
    m = hashlib.sha256()
    m.update(key_seed)
    return m.digest()

def get_polynomial_dict(C):
    res = C.dict()
    for key in res:
        res[key] = discrete_log(res[key], alpha, size-1)
    return res

human_world_size = 64
spirit_world_size_param = 32
disharmony_count = 16

H = make_human_world(human_world_size)
S = make_spirit_world(H, spirit_world_size_param)
World = H+S
D = make_disharmony(World, disharmony_count)
C = H+S+D

key = make_key(D)
cipher = AES.new( key, AES.MODE_ECB )

print("# Witch making the map! please wait.", flush=True)
witch_map = []
for i in range(spirit_world_size_param):
    row = []
    C_y = C(y=alpha^(i+1))
    for j in range(spirit_world_size_param):
        temp = C_y(x=alpha^(j+1))
        if temp == 0:
            row.append(None)
        else:
            row.append(discrete_log(temp, alpha, size-1))
    witch_map.append(row)
print("witch_map=", witch_map)
print("treasure_box=", cipher.encrypt(FLAG))

$S=H\bmod{(\prod_{i=1}^{32}(x-\alpha^{i})+\prod_{i=1}^{32}(y-\alpha^{i}))}$ より $S(\alpha^{i}, \alpha^{j})=H(\alpha^{i}, \alpha^{j})$ ( $1\leq i, j\leq 32$ ) であることに注意すると、結局次のような問題になります。

$K=\mathbb{F}_{2^{8}}=\mathbb{F}_{2}(\alpha)$ 係数の2変数多項式 $D(x,y)$ がある
$D(x,y)$ は $x$ について96次未満、 $y$ について64次未満であり、0でない項は16個である
$D(\alpha^{i}, \alpha^{j})$ ( $1\leq i, j\leq 32$ ) が与えられる
$D(x,y)$ を復元できればフラグが得られる

$G(x)=\prod_{i=1}^{32}(x-\alpha^{i})$ とします。各 $j$ に対し、 $D(x,\alpha^{j})\bmod{(x-\alpha^{i})}$ ( $i=1,\ldots,32$ ) が与えられているので、中国剰余定理により、 $D(x,\alpha^{j})\bmod{G(x)}$ が求まります。つまり、未知の多項式 $q_j(x)$ と、既知の32次未満の多項式 $r_j(x)$ を用いて

$\displaystyle{ D(x,\alpha^{j})=q_j(x)G(x)+r_j(x) }$

と表せます。

$D(x,\alpha^{j})$ たちは、次数が $j$ によらず共通な16個の項を除き係数が0となっているので、 $D(x,\alpha^{j})=\sum_{i}D_{ij}x^{i}$ と表したとき、行列 $(D_{ij})$ のrankは16以下となります。ここで、 $(D_{ij})$ のrankはちょうど16となる場合が多いと推測できます。

$(D_{ij})$ のrankがちょうど16と仮定し、係数が0でない項の次数を $i_1, \ldots, i_{16}$ とおくと、 $i$ が $i_1, \ldots, i_{16}$ のいずれかに一致する場合、

$\displaystyle{ \sum_{j=1}^{32}c_{j}D(x,\alpha^{j})=x^{i} }$

を満たす $c_j\in K$ が存在します。ここで $\bmod{G(x)}$ を考えると、

$\displaystyle{ \sum_{j=1}^{32}c_{j}r_{j}(x)=(x^{i}\bmod{G(x)}) }$

を満たす $c_j\in K$ が存在することが言えます。

$(D_{ij})$ のrankは16以下なので、17個の $j$ を取ってくると、それら $D(x,\alpha^{j})$ は1次従属となります。つまり、

$\displaystyle{ c_1D(x,\alpha^{j_1})+\cdots+c_{17}D(x,\alpha^{j_{17}})=0 }$

を満たす、すべてが0ではない $c_i\in K$ が存在します。 $\bmod{G(x)}$ を考えると、

$\displaystyle{ c_1r_{j_1}(x)+\cdots+c_{17}r_{j_{17}}(x)=0 }$

を満たす、すべてが0ではない $c_i\in K$ が存在することになります。よって、 $r_j(x)=\sum_{i}r_{ij}x^{i}$ と表したとき、行列 $(r_{ij})$ のrankも16以下となります。

$i$ が $i_1, \ldots, i_{16}$ のいずれでもない場合、 $i'=i, i_1,\ldots,i_{16}$ に対する $(x^{i'}\bmod{G(x)})$ たちは1次独立である場合が多いと推測できます。これを仮定すると、 $(r_{ij})$ のrankはちょうど16で、

$\displaystyle{ \sum_{j=1}^{32}c_jr_j(x)=(x^i\bmod{G(x)}) }$

を満たす $c_j\in K$ は存在しないと言えます。

以上より、 $i=0,1,\ldots,95$ のそれぞれについて

$\displaystyle{ \sum_{j=1}^{32}c_jr_j(x)=(x^i\bmod{G(x)}) }$

を満たす $c_j\in K$ が存在するかを、連立1次方程式を解いてチェックすることで、 $D(x,y)$ の0でない項の $x$ の次数を求めることができそうです。

同様に $D(x,y)$ の0でない項の $y$ の次数も求めることができます。これで、 $D(x,y)$ の未知の係数は $16\times 16$ 個となり、 $32\times 32$ 個の $(x,y)$ での値がわかっているので、連立1次方程式を解いて係数を求めることができます。数回サーバーに接続することで解けるケースを引くことができました。

from Crypto.Util.number import *
import Crypto.Cipher.AES as AES
import hashlib

human_world_size = 64
spirit_world_size_param = 32
disharmony_count = 16

witch_map= [[79, 234, 66, 99, 151, 3, 58, 164, 80, 54, 3, 252, 174, 48, 173, 241, 24, 217, 80, 46, 135, 23, 136, 118, 62, 151, 13, 13, 201, 95, 120, 212], ... ]
treasure_box= b'\x99{d\xfe\x06\xc9\x07\x7f\xca\xc8\x01\xff\x857l\x10\x170\x99\xeb\xc0\xc9\xb4Q\x8a0\xe18\xd3F\x95\x99\xba\x99Y\x81R\xc5\\\xfa\x0c\xb5\x10CI\x18\xe0\x96\xd9k\x1d\x884{\xd9\x01(\xac\xc0)\xd2\x1a81'

R.<x> = PolynomialRing(GF(2))
size = 2^8
K.<alpha> = GF(size, modulus=x^8+x^4+x^3+x^2+1)

def inv(a, mod):
    b, x, y=mod, 1, 0
    while b:
        t=a//b
        a-=t*b
        x-=t*y
        a, b=b, a
        x, y=y, x
    return x%mod

def crt(r1, q1, r2, q2):
    g = inv(q2, q1)
    k = g*q2%q1
    g *= K(k)^(-1)
    return (g*(r1-r2)%q1)*q2+r2

G = prod(x-alpha^i for i in range(1,spirit_world_size_param+1))

# xの次数を求める
rs = []
for j in range(spirit_world_size_param):
    r = 0
    q = 1
    for i in range(spirit_world_size_param):
        v = 0
        if witch_map[j][i] is not None:
            v = alpha^witch_map[j][i]
        r = crt(v, x-alpha^(i+1), r, q)
        q *= (x-alpha^(i+1))
    r = r.list()
    r += [0]*(spirit_world_size_param-len(r))
    rs.append(r)

rst = Matrix(K, rs).transpose()
indx = []
for i in range(spirit_world_size_param+human_world_size):
    v = (x^i%G).list()
    v += [0]*(spirit_world_size_param-len(v))
    try:
        w = rst.solve_right(vector(K, v))
    except:
        continue
    indx.append(i)

# yの次数を求める
rs = []
for j in range(spirit_world_size_param):
    r = 0
    q = 1
    for i in range(spirit_world_size_param):
        v = 0
        if witch_map[i][j] is not None:
            v = alpha^witch_map[i][j]
        r = crt(v, x-alpha^(i+1), r, q)
        q *= (x-alpha^(i+1))
    r = r.list()
    r += [0]*(spirit_world_size_param-len(r))
    rs.append(r)

rst = Matrix(K, rs).transpose()
indy = []
for i in range(spirit_world_size_param+human_world_size):
    v = (x^i%G).list()
    v += [0]*(spirit_world_size_param-len(v))
    try:
        w = rst.solve_right(vector(K, v))
    except:
        continue
    indy.append(i)

# 係数を求める
M = [[K(0)]*(disharmony_count*disharmony_count) for _ in range(spirit_world_size_param*spirit_world_size_param)]
for i in range(spirit_world_size_param):
    for j in range(spirit_world_size_param):
        for k in range(disharmony_count):
            for l in range(disharmony_count):
                M[i*spirit_world_size_param+j][k*disharmony_count+l] = alpha^((i+1)*indx[k]+(j+1)*indy[l])

v = [K(0)]*(spirit_world_size_param*spirit_world_size_param)
for i in range(spirit_world_size_param):
    for j in range(spirit_world_size_param):
        if witch_map[j][i] is not None:
            v[i*spirit_world_size_param+j] = alpha^witch_map[j][i]

c = Matrix(K, M).solve_right(vector(K, v))

R.<x, y> = PolynomialRing(GF(2))
D = 0
for i in range(disharmony_count):
    for j in range(disharmony_count):
        D += c[i*disharmony_count+j]*x^indx[i]*y^indy[j]

def make_key(D):
    key_seed = b""
    for pos, value in sorted(list(D.dict().items())):
        x = pos[0]
        y = pos[1]
        power = discrete_log(value, alpha, size-1)
        key_seed += long_to_bytes(x) + long_to_bytes(y) + long_to_bytes(power)
    m = hashlib.sha256()
    m.update(key_seed)
    return m.digest()

key = make_key(D)
cipher = AES.new( key, AES.MODE_ECB )
print(cipher.decrypt(treasure_box))
# b'SECCON{I_te4ch_y0u_secret_spell...---number_XIV---Temperance!!!}'

最初、0でない係数が少ない (係数のなすベクトルが小さい) ということでLLLを試していましたが、そもそも結果が0と1からなるベクトルにならず全然ダメでした。確かに最小のベクトルを求めてくれる保証は無さそうですが、あまりよくわかっていません。

$H$ や $S$ が出てくる意味がわからず、ヒントなのかなと思って $\bmod{G_x}$ について考えてみたところ、一気に考察が進んでとても面白かったです。

なお、 $H$ や $S$ の意図がよくわからなかったとツイートしたところ、Reed-Solomon 符号というものをもとにしているのでは、と教えていただきました。

The challenge is similar to Reed-Solomon code or error correcting code. H is the message, and S is added so that (H + S)(a^i, a^j) = 0. Then, any error E -> (H+S+E)(a^i,a^j) = E(…) -> can find E, so you can correct error in message!
— @grhkm@infosec.exchange (@grhkm21) September 17, 2023

確かに Wikipedia を見ると、 $H$ がメッセージ、 $H+S$ が符号、 $D$ がエラーに見えて、1変数多項式の場合の誤り訂正方法も書かれているので、真似すると解けたのかもしれないです。