R3CTF 2026 Writeup - Sceleri's Blog

这周跟 Project SEKAI 一起打了 R3CTF 2026，做了三道 Crypto。

teRRibleRing#

Description#

Notice that there are 3 “R”s in the title :)

1
from sage.stats.distributions.discrete_gaussian_integer import DiscreteGaussianDistributionIntegerSampler
2
from Crypto.Util.number import *
3
from secrets import randbelow
4
from secret import flag
5

6
p = 0x8000000b
7
PR.<x> = PolynomialRing(Zmod(p))
8
f = x^512 + 355853415*x^511 + ...
9

10
def F0():
11
    D = DiscreteGaussianDistributionIntegerSampler(sigma=5.0)
12
    s = PR([D() for _ in range(512)])
13
    e = PR([D() for _ in range(512)])
14
    a = PR([randbelow(p) for _ in range(512)])
15
    b = (a*s + e) % f
16
    return [a.list(), b.list()]
17

18
def F1():
19
    a = PR([randbelow(p) for _ in range(512)])
20
    b = PR([randbelow(p) for _ in range(512)])
21
    return [a.list(), b.list()]
22

23
with open("samples.txt", "w") as s:
24
    s.write("output =" + str([F0() if i == "0" else F1() for i in bin(bytes_to_long(flag))[2:].zfill(len(flag)*8)]))

Solution#

这道题给了 344 个 degree 512 的 decisional RLWE 样本。其中 $p$ 大约 $2^{31}$ ，而噪声的 std 只有 5。显然跑 344 个 1024 维的 LLL 是不现实的，所以只能是 f 有后门。

假如 $f=g\cdot h$ ，那么我们可以将 $\mathbb{Z}_p[x]/(f)$ 投影到 $\mathbb{Z}_p[x]/(g)$ ，从而在更小的维度解决这个 RLWE。但是这个投影相当于计算 $a\mod g$ ，而 $g$ 的参数可能非常大，导致噪声也会被放大。因此这个方法只有在 $g$ 的系数足够小的情况下才有用。

尝试分解 $f$ 可以发现它有一堆小因子，但它们的系数都很大。所以我们先试着猜一些因子乘起来能得到一个系数较小的多项式。假如真的存在这么一个因子，我们可以推测出如下的性质：

这个因子的度数不能太大，不然 LLL 还是跑不完，这里我们假设度数不超过 128，对应的格维度是 256。
这个因子的系数不能太大。根据我们前一条假设度数不超过 128， $x^{512}\mod g$ 会差不多把 $g$ 四次方，所以就算是十几的系数都会爆炸，必须得是非常小非常稀疏的才有戏。

现在我们希望从 $f$ 的 75 个因子里面找出乘积很小的组合。显然枚举是不现实的，但是我们可以这么想：一个因子乘上一个多项式后系数很小的概率是很低的，所以我们可以先猜 $g$ 有哪些因子，再用 LLL 反过来找这个乘上什么多项式后系数很小。这个和 LFSR 找一个稀疏的递推关系是类似的。

1
def find_short_poly(f, max_deg=128):
2
    coeffs = f.coefficients(sparse=False)
3
    coeffs = [ZZ(c) for c in coeffs]
4
    M = []
5
    for i in range(max_deg - len(coeffs) + 1):
6
        row = [0] * i + coeffs + [0] * (max_deg - len(coeffs) - i)
7
        M.append(row)
8
    M = matrix(ZZ, M)
9
    M = block_matrix([[M], [identity_matrix(max_deg)*p]])
10
    M = M.LLL(algorithm='flatter')
11
    return M[0]

考虑到 $p$ 大约 $2^{31}$ 而系数只有个位数，我们可以认为对所有度数大于 12 的因子都不会出现伪阳。于是成功发现一个度数 16 的因子的倍数里面有 $x^{127} - x^{46} + x^{19} - x^8 + 1$ 。

接下来就是把 RLWE 投影过来跑 LLL 了。这里直接简单用 Kannan Embedding 跑 cvp 即可。

\begin{pmatrix} e & 0 & b \\ 0 & I & A\\ 0 & 0 & pI \end{pmatrix}

这里 $A$ 是 $a$ 的系数矩阵，然后 $e$ 取一个和噪声差不多大的数。但是 flatter 没法解出这个 SVP，所以我换到了 BLASter 这个库，支持更多的 BKZ，每个可以 6 秒跑完。最后花了 20 多分钟跑完 344 个 bit，然后重构出 flag 即可。

1
from Crypto.Util.number import *
2
from tqdm import tqdm
3

4
p = 0x8000000b
5
PR.<x> = PolynomialRing(Zmod(p))
6
f = x^512 + 355853415*x^511 + ...
7

8
R.<x> = PolynomialRing(ZZ)
9
r = x^127 - x^46 + x^19 - x^8 + 1
10

11
rp = r.change_ring(Zmod(p))
12

13
def lift_centered(poly):
14
    return R([c.lift_centered() for c in poly.list()])
15

16
def create_ideal_lattice(a, f):
17
    x = a.parent().gen()
18
    n = f.degree()
19
    M = []
20
    for i in range(n):
21
        s = (a * x**i) % f
22
        row = lift_centered(s).list()
23
        M.append(row)
24
    return matrix(ZZ, M)
25

26
def flatter(M):
27
    from subprocess import check_output
28
    from re import findall
29
    z = "[[" + "]\n[".join(" ".join(map(str, row)) for row in M) + "]]"
30
    ret = check_output(["sage", "src/app.py", "-b42", "-t1", "-P2"], input=z.encode(), cwd="/home/sceleri/BLASter")
31
    return matrix(M.nrows(), M.ncols(), map(int, findall(b"-?\\d+", ret)))
32

33
a, b, s, e = F0()
34

35
def check(data):
36
    a, b = data
37
    a = PR(a)
38
    b = PR(b)
39
    ar = a % rp
40
    br = b % rp
41
    coef_variance = [0] * 127
42

43
    for i in range(512):
44
        c = x ** i % r
45
        for j in range(127):
46
            coef_variance[j] += c[j] ** 2
47

48
    coef_variance = [round(sqrt(c * 65536)) for c in coef_variance]
49

50
    M = create_ideal_lattice(ar, rp)
51

52
    target = matrix(ZZ, [0]*127+lift_centered(br).list())
53

54
    n = M.nrows()
55
    M = block_matrix(ZZ, [[identity_matrix(ZZ, n), M], [zero_matrix(ZZ, n, n), p*identity_matrix(ZZ, n)]])
56

57
    M = block_matrix(ZZ, [[M, 0], [target, 512]])
58

59
    M = flatter(M)
60

61
    v = M[0]
62
    cnt = 0
63
    for c in v:
64
        if abs(c) < 1024:
65
            cnt += 1
66
    if cnt > 200:
67
        return 0
68
    else:
69
        return 1
70

71
output = []
72
exec(open("samples.txt").read())
73

74
c = []
75

76
def pack_byte(ulist):
77
    ret = 0
78
    for i in range(8):
79
        ret |= (ulist[i] << (7-i))
80
    return ret
81

82
def decode_flag(c):
83
    out = b""
84
    for i in range(0, len(c), 8):
85
        if i + 8 > len(c):
86
            break
87
        out += bytes([pack_byte(c[i:i+8])])
88
    return out
89

90
from tqdm.contrib.concurrent import process_map
91

92
c = []
93
for i in tqdm(range(len(output) // 8)):
94
    u = process_map(check, [output[i*8+j] for j in range(8)])
95
    c.extend(u)
96
    print(c)
97
    print(decode_flag(c))

rECp1cG#

Description#

Quiet steps, old notes:
https://eprint.iacr.org/2007/099.pdf

1
from hashlib import sha256
2
from random import SystemRandom
3
from secrets import token_hex
4
import signal
5

6
from secret import flag
7

8
p_bits = 1024
9
k = 21
10
d_bits = 451
11
Delta = 1 << d_bits
12
solve_timeout = 888
13

14
rand = SystemRandom()
15

16
def is_prime(x):
17
    if x < 2:
18
        return False
19
    small = (2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43,
20
             47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97)
21
    for q in small:
22
        if x == q:
23
            return True
24
        if x % q == 0:
25
            return False
26

27
    odd = x - 1
28
    power = 0
29
    while odd % 2 == 0:
30
        odd //= 2
31
        power += 1
32

33
    bases = (2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43,
34
             47, 53, 59, 61)
35
    for base in bases:
36
        if base % x == 0:
37
            continue
38
        y = pow(base, odd, x)
39
        if y == 1 or y == x - 1:
40
            continue
41
        for _ in range(power - 1):
42
            y = pow(y, 2, x)
43
            if y == x - 1:
44
                break
45
        else:
46
            return False
47
    return True
48

49
def make_prime(bits):
50
    while True:
51
        x = rand.getrandbits(bits)
52
        x |= 1 << (bits - 1)
53
        x |= 3
54
        if is_prime(x):
55
            return x
56

57
def add(P, Q, a, p):
58
    if P is None:
59
        return Q
60
    if Q is None:
61
        return P
62

63
    x1, y1 = P
64
    x2, y2 = Q
65
    if x1 == x2 and (y1 + y2) % p == 0:
66
        return None
67

68
    if P == Q:
69
        slope = (3 * x1 * x1 + a) * pow(2 * y1, -1, p) % p
70
    else:
71
        slope = (y2 - y1) * pow(x2 - x1, -1, p) % p
72

73
    x3 = (slope * slope - x1 - x2) % p
74
    y3 = (slope * (x1 - x3) - y1) % p
75
    return x3, y3
76

77
def random_point(a, b, p):
78
    while True:
79
        x = rand.randrange(1, p)
80
        y2 = (x * x * x + a * x + b) % p
81
        if y2 != 0 and pow(y2, (p - 1) // 2, p) == 1:
82
            y = pow(y2, (p + 1) // 4, p)
83
            if rand.getrandbits(1):
84
                y = (-y) % p
85
            return x, y
86

87
def main():
88
    key_tag = token_hex(16)
89
    p = make_prime(p_bits)
90

91
    while True:
92
        a = rand.randrange(1, p)
93
        b = rand.randrange(1, p)
94
        if (4 * a * a * a + 27 * b * b) % p != 0:
95
            break
96

97
    while True:
98
        g = random_point(a, b, p)
99
        u = random_point(a, b, p)
100
        points = [u]
101
        for _ in range(k - 1):
102
            u = add(u, g, a, p)
103
            if u is None:
104
                break
105
            points.append(u)
106
        if len(points) == k and len({P[0] for P in points}) == k:
107
            break
108

109
    states = []
110
    for x, _ in points:
111
        while True:
112
            shown = x - rand.randrange(-Delta, Delta + 1)
113
            if 0 <= shown < p:
114
                states.append(shown)
115
                break
116

117
    size = (p.bit_length() + 7) // 8
118
    nums = (a, b, g[0], g[1], points[0][0], points[0][1])
119
    material = b"".join(int(x).to_bytes(size, "big") for x in nums)
120
    key = sha256(material + b"|" + key_tag.encode()).digest()
121

122
    pad = b""
123
    ctr = 0
124
    while len(pad) < len(flag):
125
        pad += sha256(key + ctr.to_bytes(4, "big")).digest()
126
        ctr += 1
127
    ct = bytes(x ^ y for x, y in zip(flag, pad)).hex()
128

129
    print(f"p = {p}")
130
    print(f"a = {a}")
131
    print(f"b = {b}")
132
    print(f"Delta = {Delta}")
133
    print(f"G = {g}")
134
    print(f"states = {states}")
135
    print()
136
    print("# submit P0.x", flush=True)
137

138
    signal.alarm(solve_timeout)
139
    try:
140
        answer = input().strip()
141
    except EOFError:
142
        return
143

144
    if len(answer) > 4096:
145
        print("# wrong", flush=True)
146
        return
147

148
    try:
149
        recovered_x = int(answer, 0)
150
    except ValueError:
151
        print("# wrong", flush=True)
152
        return
153

154
    if recovered_x != points[0][0]:
155
        print("# wrong", flush=True)
156
        return
157

158
    print("# ok", flush=True)
159
    print(f"key_tag = {key_tag!r}")
160
    print(f"ct = {ct!r}", flush=True)
161

162
if __name__ == "__main__":
163
    main()

Solution#

一眼看上去就是论文题，但是如果顺着题目里给的论文找就被阴了。这里我们需要找的是只给 $x_i$ MSB 的 Coppersmith 构造，然后应该搜索的是 ECHNP。在 Google Scholar 看 cite 可以找到 2026 年的最新工作¹。然后就简单了，ocr 成 latex 交给 AI 去抄那一堆多项式构造即可，而且代入测试数据检查 $\mkern{-20mu}\mod p^2$ 是否为 0 就可以验证正确性。

实际跑的时候发现对于 451 bits 的 delta，LLL reduce 后的多项式大概会落在 $-3p^2\sim 3p^2$ 的范围内，但是只要不是 0 就没法用 grobner 解。这个 Coppersmith 有 11 个变量，所以直接暴力猜测是不行的（除非你愿意跑 1000 次 408 维的 LLL），所以需要一点小优化就是 coppersmith base polys 里面有很多都包含了 $x$ ，所以正确做法是在设置 bounds 的时候把 $x$ 的 bound 除以 8，这样就只要猜 8 次即可。

下面的代码是清理过的，同时把 GPT 写的 grobner 换成了 msolve。²

1
from sage.all import *
2
from hashlib import sha256
3
from itertools import combinations
4
import argparse
5
import ast
6
import re
7
import shutil
8
import subprocess
9
from pwn import context, remote, process
10

11
n = 10
12
d = 2
13
t = 1
14
mid = n
15

16
def ec_add(P, Q, a, p):
17
    if P is None:
18
        return Q
19
    if Q is None:
20
        return P
21
    x1, y1 = P
22
    x2, y2 = Q
23
    if x1 == x2 and (y1 + y2) % p == 0:
24
        return None
25
    if P == Q:
26
        lam = (3 * x1 * x1 + a) * inverse_mod(2 * y1, p) % p
27
    else:
28
        lam = (y2 - y1) * inverse_mod(x2 - x1, p) % p
29
    x3 = (lam * lam - x1 - x2) % p
30
    y3 = (lam * (x1 - x3) - y1) % p
31
    return (ZZ(x3), ZZ(y3))
32

33
def ec_neg(P, p):
34
    if P is None:
35
        return None
36
    return (P[0], (-P[1]) % p)
37

38
def ec_mul(k, P, a, p):
39
    R = None
40
    Q = P
41
    k = ZZ(k)
42
    while k:
43
        if k & 1:
44
            R = ec_add(R, Q, a, p)
45
        Q = ec_add(Q, Q, a, p)
46
        k >>= 1
47
    return R
48

49
def as_zz_data(data):
50
    out = dict(data)
51
    for k in ("p", "a", "b", "Delta"):
52
        out[k] = ZZ(out[k])
53
    out["G"] = (ZZ(out["G"][0]), ZZ(out["G"][1]))
54
    out["states"] = [ZZ(x) for x in out["states"]]
55
    if "points" in out:
56
        out["points"] = [(ZZ(P[0]), ZZ(P[1])) for P in out["points"]]
57
    if "errors" in out:
58
        out["errors"] = [ZZ(e) for e in out["errors"]]
59
    return out
60

61
def signed_mod(x, m):
62
    x = ZZ(x) % m
63
    if x > m // 2:
64
        x -= m
65
    return ZZ(x)
66

67
class ECXHNPSolver:
68
    def __init__(self, data, check=False, verbose=True):
69
        self.data = as_zz_data(data)
70
        self.p = self.data["p"]
71
        self.a = self.data["a"]
72
        self.b = self.data["b"]
73
        self.Delta = self.data["Delta"]
74
        self.G = self.data["G"]
75
        self.states = self.data["states"]
76
        self.check = False
77
        self.verbose = verbose
78
        if len(self.states) < 2*n + 1:
79
            raise ValueError("need 21 states for n=10")
80

81
        names_z = ["x"] + ["z%d" % i for i in range(1, n + 1)]
82
        names_y = ["x"] + ["y%d" % i for i in range(1, n + 1)]
83
        self.Rz = PolynomialRing(ZZ, names_z, order="degrevlex")
84
        self.Ry = PolynomialRing(ZZ, names_y, order="degrevlex")
85
        self.RQy = PolynomialRing(QQ, names_y, order="degrevlex")
86
        self.xz = self.Rz.gen(0)
87
        self.zs = [None] + [self.Rz.gen(i) for i in range(1, n + 1)]
88
        self.xy = self.Ry.gen(0)
89
        self.ys = [None] + [self.Ry.gen(i) for i in range(1, n + 1)]
90

91
        self.Q = [None] + [ec_mul(i, self.G, self.a, self.p) for i in range(1, n + 1)]
92
        self.h0 = self.states[mid]
93
        self.A = [None]
94
        self.B = [None]
95
        self.C = [None]
96
        self.D = [None]
97
        self.H = [None]
98
        self.F = [None]
99
        self._L_cache = {}
100
        self._L0_triples = {}
101
        self._tilde_cache = {}
102
        self._over_cache = {}
103
        self._build_base_polys()
104

105
    def log(self, s):
106
        if self.verbose:
107
            print(s, flush=True)
108

109
    def _build_base_polys(self):
110
        x = self.xz
111
        for i in range(1, n + 1):
112
            xq, yq = self.Q[i]
113
            tilde_h = self.states[mid + i] + self.states[mid - i]
114
            Ai = self.h0 - xq
115
            Bi = -2 * (self.a + 3 * xq**2)
116
            Ci = Ai * Bi - 4 * yq**2
117
            Di = tilde_h - 2 * xq
118
            self.A.append(ZZ(Ai))
119
            self.B.append(ZZ(Bi))
120
            self.C.append(ZZ(Ci))
121
            self.D.append(ZZ(Di))
122
            self.H.append((x + Ai)**2 * self.zs[i] + Bi*x + Ci)
123
        for i in range(1, n + 1):
124
            self.F.append(self.to_y(self.H[i]))
125

126
    def _check_base_polys(self):
127
        for i in range(1, n + 1):
128
            self.assert_z(self.H[i], self.p, "H_%d" % i)
129
            self.assert_y(self.F[i], self.p, "F_%d" % i)
130
        self.log("[check] H_i/F_i roots mod p: ok")
131

132
    def assert_z(self, poly, mod, label):
133
        val = ZZ(poly(*self.zroot))
134
        if val % mod != 0:
135
            raise AssertionError("%s(root) mod %s != 0" % (label, mod))
136

137
    def assert_y(self, poly, mod, label):
138
        val = ZZ(poly(*self.yroot))
139
        if val % mod != 0:
140
            raise AssertionError("%s(root) mod %s != 0" % (label, mod))
141

142
    def to_y(self, poly):
143
        args = [self.xy] + [self.ys[i] + self.D[i] for i in range(1, n + 1)]
144
        return self.Ry(poly(*args))
145

146
    def reduce_poly(self, poly, mod, ring):
147
        terms = {}
148
        for exp, coeff in poly.dict().items():
149
            c = ZZ(coeff) % mod
150
            if c:
151
                terms[exp] = c
152
        return ring(terms)
153

154
    def H_vector(self, js):
155
        js = tuple(js)
156
        vec = []
157
        x = self.xz
158
        for m, jm in enumerate(js):
159
            term = self.zs[jm]
160
            for r, jr in enumerate(js):
161
                if r != m:
162
                    term *= self.H[jr]
163
            vec.append(term)
164
        for m, jm in enumerate(js):
165
            term = x * self.zs[jm]
166
            for r, jr in enumerate(js):
167
                if r != m:
168
                    term *= self.H[jr]
169
            sub = self.Rz(0)
170
            for u, ju in enumerate(js):
171
                if u == m:
172
                    continue
173
                part = self.Rz(self.B[ju])
174
                for r, jr in enumerate(js):
175
                    if r != u:
176
                        part *= self.H[jr]
177
                sub += part
178
            vec.append(term - sub)
179
        return vec
180

181
    def M_matrix(self, js, mod):
182
        x = self.xz
183
        rows = []
184
        for shift in (0, 1):
185
            for m, jm in enumerate(js):
186
                poly = self.Rz(x**shift)
187
                for r, jr in enumerate(js):
188
                    if r != m:
189
                        poly *= (x + self.A[jr])**2
190
                row = [ZZ(poly.monomial_coefficient(x**k)) % mod for k in range(2*len(js))]
191
                rows.append(row)
192
        return Matrix(Zmod(mod), rows)
193

194
    def L(self, js, m):
195
        js = tuple(sorted(js))
196
        key = (js, m)
197
        if key in self._L_cache:
198
            return self._L_cache[key]
199
        l = len(js)
200
        mod = self.p ** (l - 1)
201
        W = self.M_matrix(js, mod).inverse()
202
        hv = self.H_vector(js)
203
        poly = self.Rz(0)
204
        for c, h in zip(list(W[m]), hv):
205
            poly += ZZ(c) * h
206
        poly = self.reduce_poly(poly, mod, self.Rz)
207
        self._L_cache[key] = poly
208
        if self.check:
209
            self.assert_z(poly, mod, "L_%d%s" % (m, js))
210
        return poly
211

212
    def coeff_xz(self, poly, xpow, zset):
213
        exp = [0] * (n + 1)
214
        exp[0] = xpow
215
        for j in zset:
216
            exp[j] = 1
217
        return ZZ(poly.monomial_coefficient(self.Rz.monomial(*exp)))
218

219
    def overline_L0(self, js):
220
        js = tuple(sorted(js))
221
        if js in self._over_cache:
222
            return self._over_cache[js]
223
        poly = self.L(js, 0)
224
        for pair in combinations(js, d):
225
            c = self.coeff_xz(poly, 2*d - 1, pair)
226
            if c:
227
                poly -= c * self.L(pair, 2*d - 1)
228
        poly = self.reduce_poly(poly, self.p, self.Rz)
229
        self._over_cache[js] = poly
230
        if self.check:
231
            self.assert_z(poly, self.p, "overline_L0%s" % (js,))
232
        return poly
233

234
    def ordered_triples(self):
235
        triples = list(combinations(range(1, n + 1), d + 1))
236
        triples.sort(key=lambda c: tuple(1 if i in c else 0 for i in range(n, 0, -1)))
237

238
        pairs = list(combinations(range(1, n + 1), d))
239
        rows = []
240
        for tri in triples:
241
            L0 = self.L(tri, 0)
242
            self._L0_triples[tri] = L0
243
            rows.append([self.coeff_xz(L0, 2*d - 1, pair) % self.p for pair in pairs])
244

245
        M = Matrix(GF(self.p), rows).transpose()
246
        pivots = list(M.pivots())
247
        if len(pivots) < binomial(n, d):
248
            raise ValueError("could not find enough independent triples")
249
        pivots = pivots[:binomial(n, d)]
250
        pivot_set = set(pivots)
251
        ordered = [triples[i] for i in pivots] + [
252
            triples[i] for i in range(len(triples)) if i not in pivot_set
253
        ]
254
        self.triples = ordered
255
        self.pairs = pairs
256

257
        q = self.p**d
258
        Arows = []
259
        for tri in self.triples[:binomial(n, d)]:
260
            L0 = self._L0_triples[tri]
261
            Arows.append([self.coeff_xz(L0, 2*d - 1, pair) % q for pair in pairs])
262
        self.elim_A = Matrix(Zmod(q), Arows)
263
        _ = self.elim_A.inverse()
264
        self.log("[construct] selected %d invertible L0 triples for Construction II" % len(Arows))
265

266
    def tilde_L(self, i0, v):
267
        key = (i0, v)
268
        if key in self._tilde_cache:
269
            return self._tilde_cache[key]
270
        tri = self.triples[v]
271
        q = self.p**d
272
        target = vector(Zmod(q), [
273
            self.coeff_xz(self.L(tri, i0), 2*d - 1, pair) % q
274
            for pair in self.pairs
275
        ])
276
        e = self.elim_A.transpose().solve_right(target)
277
        poly = self.L(tri, i0)
278
        for u in range(binomial(n, d)):
279
            eu = ZZ(e[u])
280
            if eu:
281
                poly -= eu * self.L(self.triples[u], 0)
282
        poly = self.reduce_poly(poly, q, self.Rz)
283
        self._tilde_cache[key] = poly
284
        if self.check:
285
            self.assert_z(poly, q, "tilde_L_%d_%s" % (i0, tri))
286
        return poly
287

288
    def construct_lattice_polys(self):
289
        self.ordered_triples()
290
        x = self.xy
291
        polys = []
292
        labels = []
293
        p = self.p
294

295
        for i0 in range(0, 2*d - 1):
296
            polys.append(p**d * x**i0)
297
            labels.append("I(%d)" % i0)
298

299
        for j in range(1, n + 1):
300
            for i0 in range(0, 2):
301
                polys.append(p**d * x**i0 * self.ys[j])
302
                labels.append("IIa(%d,%d)" % (i0, j))
303
            polys.append(p**(d - 1) * self.F[j])
304
            labels.append("IIb(2,%d)" % j)
305

306
        for pair in combinations(range(1, n + 1), 2):
307
            for i0 in range(0, 2*d - 1):
308
                polys.append(p**(d - 2 + 1) * self.to_y(self.L(pair, i0)))
309
                labels.append("IIIa(%d,%s)" % (i0, pair))
310

311
        cutoff = binomial(n, d)
312
        for v, tri in enumerate(self.triples):
313
            if v < cutoff:
314
                polys.append(p * self.to_y(self.overline_L0(tri)))
315
                labels.append("IVa(0,%s)" % (tri,))
316
            else:
317
                polys.append(self.to_y(self.tilde_L(0, v)))
318
                labels.append("IVb(0,%s)" % (tri,))
319
        for i0 in range(1, t + 1):
320
            for v, tri in enumerate(self.triples):
321
                polys.append(self.to_y(self.tilde_L(i0, v)))
322
                labels.append("IVb(%d,%s)" % (i0, tri))
323

324
        expected_dim = (t + 1) * binomial(n, d + 1) + (2*d - 1) * sum(
325
            binomial(n, l) for l in range(d + 1)
326
        )
327
        if len(polys) != expected_dim:
328
            raise AssertionError("dimension mismatch: got %d, expected %d" % (len(polys), expected_dim))
329

330
        if self.check:
331
            for label, poly in zip(labels, polys):
332
                self.assert_y(poly, p**d, label)
333
            self.log("[check] all %d lattice polynomials vanish mod p^2: ok" % len(polys))
334
        self.polys = polys
335
        return polys
336

337
    def scale_poly_terms(self, poly, bounds):
338
        terms = {}
339
        for exp, coeff in poly.dict().items():
340
            c = ZZ(coeff)
341
            for ei, bi in zip(exp, bounds):
342
                if ei:
343
                    c *= bi**ei
344
            if c:
345
                terms[exp] = c
346
        return terms
347

348
    def build_lattice(self, bounds):
349
        polys = self.construct_lattice_polys()
350
        term_rows = [self.scale_poly_terms(f, bounds) for f in polys]
351
        monoms = sorted(set().union(*[set(r.keys()) for r in term_rows]))
352
        idx = {m: i for i, m in enumerate(monoms)}
353
        M = Matrix(ZZ, len(term_rows), len(monoms))
354
        for r, row in enumerate(term_rows):
355
            for exp, coeff in row.items():
356
                M[r, idx[exp]] = coeff
357
        self.monoms = monoms
358
        self.bounds = bounds
359
        self.log("[lattice] matrix %d x %d" % (M.nrows(), M.ncols()))
360
        return M
361

362
    def vector_to_poly(self, row):
363
        terms = {}
364
        for coeff, exp in zip(row, self.monoms):
365
            coeff = ZZ(coeff)
366
            if not coeff:
367
                continue
368
            scale = ZZ(1)
369
            for ei, bi in zip(exp, self.bounds):
370
                if ei:
371
                    scale *= bi**ei
372
            if coeff % scale != 0:
373
                return None
374
            terms[exp] = coeff // scale
375
        return self.Ry(terms)
376

377
    def reduced_polys(self, bounds, max_rows):
378
        M = self.build_lattice(bounds)
379
        self.log("[LLL] reducing...")
380
        B = self.reduce_matrix(M)
381
        self.log("[LLL] done")
382
        if max_rows is None or max_rows <= 0:
383
            row_limit = B.nrows()
384
        else:
385
            row_limit = min(max_rows, B.nrows()-1)
386
        out = []
387
        decoded = 0
388
        zero_over_zz = 0
389
        for i in range(row_limit):
390
            f = self.vector_to_poly(B.row(i))
391
            if f is None or f == 0:
392
                continue
393
            decoded += 1
394
            if self.check:
395
                val = ZZ(f(*self.yroot))
396
                if val == 0:
397
                    zero_over_zz += 1
398
                    out.append(f)
399
            else:
400
                out.append(f)
401
        if self.check:
402
            self.log(
403
                "[LLL] ZZ root check: %d/%d decoded rows vanish (checked %d/%d rows)"
404
                % (zero_over_zz, decoded, row_limit, B.nrows())
405
            )
406
        self.log("[LLL] collected %d candidate integer polynomials" % len(out))
407
        return out
408

409
    def reduce_matrix(self, M):
410
        rows = []
411
        for i in range(M.nrows()):
412
            rows.append("[" + " ".join(str(ZZ(x)) for x in M.row(i)) + "]")
413
        inp = "[" + "\n".join(rows) + "\n]\n"
414
        proc = subprocess.run(
415
            ["flatter"],
416
            input=inp,
417
            text=True,
418
            capture_output=True,
419
            check=True,
420
            timeout=840,
421
        )
422
        parsed = []
423
        for line in proc.stdout.splitlines():
424
            line = line.strip()
425
            if not line or line == "]":
426
                continue
427
            line = line.lstrip("[").rstrip("]")
428
            if line:
429
                parsed.append([ZZ(x) for x in line.split()])
430
        if len(parsed) == M.nrows():
431
            self.log("[LLL] used flatter")
432
            return Matrix(ZZ, parsed)
433

434
    def solve_root(self, bounds, max_rows):
435
        short_polys = self.reduced_polys(bounds, max_rows=max_rows)
436
        if len(short_polys) < n + 1:
437
            raise ValueError("not enough short polynomials")
438

439
        ideal = self.RQy.ideal(short_polys)
440
        self.log("[root] computing Groebner basis...")
441
        import signal
442
        signal.alarm(120)
443
        root = ideal.variety(algorithm="msolve", proof=False)
444
        print(root)
445
        if not root:
446
            raise ValueError("variety returned no roots")
447
        return ZZ(root[0][self.RQy("x")])
448

449
    def recover_p0_x(self, e0):
450
        x_mid = (self.h0 + e0) % self.p
451
        rhs = (x_mid**3 + self.a*x_mid + self.b) % self.p
452
        if kronecker(rhs, self.p) != 1:
453
            raise ValueError("recovered middle x is not on curve")
454
        y_mid = ZZ(power_mod(rhs, (self.p + 1) // 4, self.p))
455
        candidates = []
456
        tenG = ec_mul(mid, self.G, self.a, self.p)
457
        for y in (y_mid, (-y_mid) % self.p):
458
            Pmid = (x_mid, y)
459
            P0 = ec_add(Pmid, ec_neg(tenG, self.p), self.a, self.p)
460
            candidates.append(P0[0])
461
        near = [x for x in candidates if abs(ZZ(x) - self.states[0]) <= self.Delta]
462
        if near:
463
            return ZZ(near[0])
464
        return ZZ(candidates[0])
465

466
def parse_remote_data(blob):
467
    text = blob.decode(errors="replace")
468

469
    def grab_int(name):
470
        m = re.search(r"^%s\s*=\s*([0-9]+)\s*$" % re.escape(name), text, re.M)
471
        if not m:
472
            raise ValueError("could not parse %s" % name)
473
        return ZZ(m.group(1))
474

475
    def grab_literal(name):
476
        m = re.search(r"^%s\s*=\s*(.+)\s*$" % re.escape(name), text, re.M)
477
        if not m:
478
            raise ValueError("could not parse %s" % name)
479
        return ast.literal_eval(m.group(1))
480

481
    return {
482
        "p": grab_int("p"),
483
        "a": grab_int("a"),
484
        "b": grab_int("b"),
485
        "Delta": grab_int("Delta"),
486
        "G": grab_literal("G"),
487
        "states": grab_literal("states"),
488
    }
489

490
def solve_data(data, max_rows):
491
    solver = ECXHNPSolver(data, check=False)
492
    bounds = [solver.Delta//8] + [solver.Delta] * n
493
    e0 = solver.solve_root(bounds, max_rows=max_rows)
494
    p0x = solver.recover_p0_x(e0)
495
    print("P0.x =", p0x)
496
    return p0x
497

498
def main():
499
    context.log_level = "debug"
500
    io = process(["sage", "challenge.py"])
501
    # io = remote("challenge.ctf2026.r3kapig.com", 30388)
502
    blob = io.recvuntil(b"# submit P0.x")
503
    data = parse_remote_data(blob)
504
    print("[remote] parsed p_bits=%d, states=%d" % (ZZ(data["p"]).bit_length(), len(data["states"])))
505
    p0x = solve_data(data, max_rows=120)
506
    io.sendline(str(p0x).encode())
507
    io.recvuntil(b"key_tag = ")
508
    key_tag = ast.literal_eval(io.recvline().strip().decode())
509
    io.recvuntil(b"ct = ")
510
    ct = ast.literal_eval(io.recvline().strip().decode())
511
    ct = bytes.fromhex(ct)
512

513
    a = data["a"]
514
    b = data["b"]
515
    g = data["G"]
516
    p = data["p"]
517
    EC = EllipticCurve(GF(p), [a, b])
518
    p0y = EC.lift_x(p0x).y()
519

520
    def try_decrypt(a, b, g, p0x, p0y):
521
        nonlocal ct, key_tag
522
        size = (p.bit_length() + 7) // 8
523
        nums = (a, b, g[0], g[1], p0x, p0y)
524
        # print(nums)
525
        material = b"".join(int(x).to_bytes(size, "big") for x in nums)
526
        key = sha256(material + b"|" + key_tag.encode()).digest()
527
        pad = b""
528
        ctr = 0
529
        while len(pad) < len(ct):
530
            pad += sha256(key + ctr.to_bytes(4, "big")).digest()
531
            ctr += 1
532
        print(bytes(x ^^ y for x, y in zip(ct, pad)))
533

534
    try_decrypt(a, b, g, p0x, p0y)
535
    try_decrypt(a, b, g, p0x, -p0y % p)
536

537
if __name__ == "__main__":
538
    main()

Encrypted Activation#

Description#

The remote server uses a fixed key. The evaluation keys are provided in the attachments.

1
#!/usr/bin/env python3
2
from __future__ import annotations
3
import base64, json, os, random, sys
4
import fhe_core as fhe
5
from secret import FLAG
6
import signal
7

8
S = 4
9
N_DIGITS = 5
10
LUT_SIZE = S ** N_DIGITS
11
SETUP_DIR = os.path.join(os.path.dirname(__file__), "setup")
12
ROUNDS = 16
13

14
def _timeout(_signum, _frame):
15
    print("timeout")
16
    sys.stdout.flush()
17
    raise SystemExit(0)
18

19
def extract_radix(v: int, s: int, n: int):
20
    d = []
21
    for _ in range(n):
22
        d.append(v % s)
23
        v //= s
24
    return d
25

26
def combine_radix(d, s: int):
27
    v = 0
28
    for x in reversed(d):
29
        v = v * s + x
30
    return v
31

32
def write_setup():
33
    os.makedirs(SETUP_DIR, exist_ok=True)
34
    print("[*] generating keys...")
35
    sk, bsk, ksk = fhe.keygen()
36

37
    client_blob = fhe.serialize_client_key(sk)
38
    with open(os.path.join(SETUP_DIR, "client.bin"), "wb") as f:
39
        f.write(client_blob)
40
    print(f"[*] client.bin: {len(client_blob)/1e6:.1f} MB")
41

42
    bsk_blob = fhe.serialize_bsk(bsk)
43
    with open(os.path.join(SETUP_DIR, "bsk.bin"), "wb") as f:
44
        f.write(bsk_blob)
45
    print(f"[*] bsk.bin: {len(bsk_blob)/1e6:.1f} MB")
46

47
    ksk_blob = fhe.serialize_ksk(ksk)
48
    with open(os.path.join(SETUP_DIR, "ksk.bin"), "wb") as f:
49
        f.write(ksk_blob)
50
    print(f"[*] ksk.bin: {len(ksk_blob)/1e6:.1f} MB")
51
    print("[*] setup complete ->", SETUP_DIR)
52
    return
53

54
def load_setup():
55
    client_path = os.path.join(SETUP_DIR, "client.bin")
56
    if not (os.path.exists(client_path)):
57
        write_setup()
58
    with open(client_path, "rb") as f:
59
        client_blob = f.read()
60
    sk = fhe.parse_client_key(client_blob)
61
    return sk
62

63
def main() -> int:
64
    signal.signal(signal.SIGALRM, _timeout)
65
    signal.alarm(120)
66
    sk = load_setup()
67
    lut = [int(token) for token in open("lut", "r").read().split()]
68
    if len(lut) != LUT_SIZE:
69
        raise ValueError(f"expected {LUT_SIZE} LUT entries, got {len(lut)}")
70
    print("=== Encrypted Activation-Layer Inference ===")
71
    print('''Each ciphertext carries one base-%d activation symbol (2 bits).
72
Evaluate the published %d-bit activation table on the encrypted
73
input using the provided keys, and return the encrypted output
74
symbols.''' % (S, (S.bit_length()-1) * N_DIGITS))
75
    sys.stdout.flush()
76

77
    for idx in range(ROUNDS):
78
        x = fhe.rng.randrange(LUT_SIZE)
79
        in_digits = extract_radix(x, S, N_DIGITS)
80
        seeded_inputs = [fhe.encrypt_ciphertext(sk, d) for d in in_digits]
81
        payload = {
82
            "round": idx + 1,
83
            "ciphertext": [(int.from_bytes(seed, 'big'), ct.b) for seed, ct in seeded_inputs],
84
        }
85
        print(json.dumps(payload))
86
        sys.stdout.flush()
87

88
        line = sys.stdin.readline()
89
        if not line:
90
            return 1
91
        line = line.strip()
92
        parts = [t for t in line.replace(",", " ").split() if t]
93
        if len(parts) != N_DIGITS:
94
            print("wrong"); sys.stdout.flush(); return 0
95

96
        try:
97
            out_cts = [fhe.parse_lwe_ciphertext(base64.b64decode(t, validate=True)) for t in parts]
98
        except Exception:
99
            print("wrong"); sys.stdout.flush(); return 0
100

101
        bad = False
102
        for c in out_cts:
103
            if len(c.a) != fhe.N:
104
                bad = True; break
105
            zero_cnt = sum(1 for v in c.a if v == 0)
106
            if zero_cnt > 16:
107
                bad = True; break
108
            if len(set(c.a)) == 1:
109
                bad = True; break
110
        if bad:
111
            print("wrong"); sys.stdout.flush(); return 0
112

113
        out_digits = [fhe.decrypt_ciphertext(sk, c) for c in out_cts]
114
        if any(d >= S for d in out_digits):
115
            print("wrong"); sys.stdout.flush(); return 0
116

117
        y = combine_radix(out_digits, S)
118
        if y != lut[x]:
119
            print("wrong"); sys.stdout.flush(); return 0
120

121
    print(FLAG)
122
    sys.stdout.flush()
123
    return 0
124

125
if __name__ == "__main__":
126
    raise SystemExit(main())

Solution#

WARNING
这道题 95% 是 ai 做的，而且我不懂 TFHE，所以以下内容可能是错的，就当看个乐子。

TFHE 但是只给了加密解密，让 ai 先分析了 bsk 和 ksk，结论是和 TFHE 论文里的不同，任何非线性的计算只能使用 programmable bootstrapping 来实现。TFHE 原论文只给了 $2^{32}$ 的实现，但这道题是 $2^{64}$ ，所以还得自己写。于是我就先让 ai 实现了一个 programmable bootstrapping 的实现，然后让 ai 优化了一下速度。比如让 ai 把 TFHE 原论文用的 nayuki 和 fftw 都试着接进来，但是 64 bits 的用 double fft 精度问题太大了。总之 fft 也没带来多少提升，最终大约 250 ms 一次 bootstrap。

然后是处理这个 LUT 的问题。每次 bootstrap 只能使用 mod 5 的 LUT 变换。从信息论的角度来说，这么干得 1000 次以上才能完成 1024 的 LUT，但题目设置了平均每个 7.5 秒的时限，所以加上多线程也就最多几百次的 bootstrap。

然后 ai 就开始发力了，说 TFHE 是把信息放在高位的（这个看解密就可以看出来），但是 bootstrap 的算法其实对这个 mod 5 没有感知，所以可以变成 mod 1024 的 ct 之后再 LUT（我内心：啊？），但很可惜的是测试炸了，原因是噪声太大了和 negacyclic 的问题。

但是换成 mod 32 之后成功率就还可以。于是就把 1024 的 LUT 分成 32 个 mod 32 的 LUT，最后再用 mask 合并成 5 个 mod 5 的输出。这个 mask 的实现也很神秘，这个 fhe 没有乘法，所以实际是把 mask 和 LUT 得到的 5 个 2 bit 数分别加了起来，相当于一个 8 slot 的 torus 再进行 bootstrap LUT。

最后大概 400 多次 bootstrap 就可以完成，但是每次的成功率大概是 70%，然后找了一台 128 核的服务器碰运气就过了。

Xu, J., Sarkar, S., Wang, H., & Hu, L. (2026). New Results on Elliptic Curve Hidden Number Problem for ECDH Key Exchange: J. Xu et al. Journal of Cryptology, 39(2), 20. ↩
但还是用了 500 行，感觉除了没有用 sage 的 EC 之外都没法变短。如果没有 ai 在 ctf 的时候写 500 行是绝对的折磨，只能说感谢 ai。 ↩

teRRibleRing#

Description#

Solution#

rECp1cG#

Description#

Solution#

Encrypted Activation#

Description#

Solution#

Footnotes#