SekaiCTF 2024 Writeup - Sceleri's Blog

I played as Twilight Light and solved 3.5 crypto challenges. The crypto is very hard and I can’t fully understand why everything works. So this writeup will mainly focus on how I found these properties heuristically.

はやぶさ#

The challenge is a standard Falcon implementation with only degree 64(too small), which makes it possible to break the lattice using BKZ/flatter.

The verification process of Falcon involves taking a polynomial $f(x)$ and checking whether the coefficients of $f(x)$ and $f(x)h(x)-\mathrm{hashed}$ are all small. So we can build the lattice like this,

\begin{pmatrix} I_{64}&M\\ 0&pI_{64} \end{pmatrix}

where the $i$ -th row of $M$ if the coeficients of $x^ih(x)$ . After the reduction, we only need to find a cvp to [0]*64+hashed.

1
from falcon import falcon
2
from falcon.encoding import compress, decompress
3
from falcon.ntt import mul_zq
4
from flag import flag
5
import json
6
from pwn import *
7

8
def flatter(M):
9
    from subprocess import check_output
10
    from re import findall
11
    z = "[[" + "]\n[".join(" ".join(map(str, row)) for row in M) + "]]"
12
    ret = check_output(["flatter"], input=z.encode())
13
    return matrix(M.nrows(), M.ncols(), map(int, findall(b"-?\\d+", ret)))
14

15
def babai_cvp(B, t, perform_reduction=True):
16
    if perform_reduction:
17
        B = B.LLL(delta=0.75)
18

19
    G = B.gram_schmidt()[0]
20
    b = t
21
    for i in reversed(range(B.nrows())):
22
        c = ((b * G[i]) / (G[i] * G[i])).round()
23
        b -= c * B[i]
24
    return t - b
25

26
sk = falcon.SecretKey(64)
27
pk = falcon.PublicKey(sk)
28

29
io = remote("hayabusa.chals.sekai.team", 1337, ssl=True)
30

31
io.recvuntil(b"h = ")
32
h = json.loads(io.recvline().strip().decode())
33
pk.h = h
34
sk.h = h
35

36
q = 12 * 1024 + 1
37

38
def one_hot(i):
39
    ret = [0]*64
40
    ret[i] = 1
41
    return ret
42

43
B0 = matrix([mul_zq(one_hot(i), h) for i in range(64)])
44
B = block_matrix(ZZ, [[identity_matrix(64), -B0], [zero_matrix(64), identity_matrix(64)*q]])
45
B = flatter(B)
46

47
while 1:
48
    salt = b"\x61" * 40
49

50
    hashed = sk.hash_to_point(b"Can you break me", salt)
51

52
    v_h = vector([0]*64+hashed)
53

54
    re = v_h - babai_cvp(B, v_h, perform_reduction=False)
55
    print("cvp:", re)
56

57
    v0 = vector(GF(q), re[:64])
58

59
    print(list(re[:64]))
60

61
    print(mul_zq(list(re[:64]), h))
62
    print(hashed)
63

64
    fake_sig = b"\x36" + salt + compress(list(re[:64]), 122-41)
65

66
    if pk.verify(b"Can you break me", fake_sig):
67
        print("well done!!")
68
        break
69

70
io.sendline(fake_sig.hex())
71
io.interactive()

マスタースパーク#

This is a CSIDH challenge with non-standard parameters. It’ll first ask you to provide a prime number with some restrictions and then perform a key exchange with a point on the curve. However, one of the points is multiplied with a secret.

The fun fact is that isogeny will keep the group structure, so secret*P=Q still holds after the isogeny.

The rest part is straightforward, just get the modulus of secret for many small p and recover it using CRT. For CSIDH, the order of the curve is p+1=4*prod(prime_list), which is exactly the small primes we provided. Therefore, Pohlig–Hellman is enough to solve the discrte log. However, the sign of the discrete log may flip due to the x-axis can’t full determine a point on the curve, and the solution is just finding some big p and enumerate all possible flips.

1
from pwn import *
2
import os
3

4
context.log_level = 'debug'
5
used_primes = set()
6

7
def sample_p():
8
    global used_primes
9
    def sample_new_prime(r):
10
        while True:
11
            p = random_prime(r)
12
            if p == 2:
13
                continue
14
            if p not in used_primes:
15
                return p
16

17
    while True:
18
        plist = []
19
        for _ in range(9):
20
            p = sample_new_prime(512)
21
            plist.append(p)
22
        if len(set(plist)) != 9:
23
            continue
24
        p0 = sample_new_prime(2**16)
25
        if p0 in plist:
26
            continue
27
        plist.append(p0)
28
        plist.append(p0)
29
        q = 4*prod(plist) - 1
30
        if is_prime(q) and ((q + 1) // 4) % 2 == 1 and q.bit_length() <= 96:
31
            used_primes = used_primes.union(plist)
32
            print(q.bit_length())
33
            return q
34

35
io = process(['sage', 'challenge.sage'])
36

37
def get_montgomery(Fp2, G):
38
    A = (G[1]**2 - G[0]**3 - G[0]) / (G[0]**2)
39
    return EllipticCurve(Fp2, [0, A, 0, 1, 0])
40

41
modlist = []
42

43
def get_PQ(p):
44
    Fp = GF(p)
45
    Fp2.<j> =  GF(p ^ 2, modulus=x ^ 2 + 1)
46

47
    io.recvuntil(b'input your prime number or secret > ')
48
    io.sendline(str(p).encode())
49
    P0 = eval(io.recvline().strip())
50
    Q0 = eval(io.recvline().strip())
51
    E = get_montgomery(Fp2, P0)
52
    P = E(*P0)
53
    Q = E(*Q0)
54
    dlog = Q.log(P)
55
    order = P.order()//4
56
    print(order.bit_length())
57
    modlist.append((dlog, order))
58

59
for _ in range(6):
60
    get_PQ(sample_p())
61

62
rs, ps = zip(*modlist)
63
rs = list(rs)
64
ps = list(ps)
65

66
for u in range(2**len(rs)):
67
    new_rs = [(-1)**((u >> i) & 1) * r for i, r in enumerate(rs)]
68
    secret = crt(new_rs, ps)
69
    if secret.bit_length() <= 256:
70
        print(secret)
71
        io.sendline(str(secret).encode())
72

73
io.interactive()

Squares vs. Cubes#

The 3 hardest challenges are authored by Neobeo. All of them are intersesting and I really like(not solving) them. Here is Neobeo’s writeup.

Description#

1
from Crypto.Util.number import bytes_to_long, getPrime
2
from secrets import token_bytes, randbelow
3
from flag import FLAG
4

5
padded_flag = bytes_to_long(FLAG + token_bytes(128 - len(FLAG)))
6

7
p, q, r = getPrime(512), getPrime(512), getPrime(512)
8
N = e = p * q * r
9
phi = (p - 1) * (q - 1) * (r - 1)
10
d = pow(e, -1, phi)
11

12
# Genni likes squares and SBG likes cubes. Let's calculate their values
13
value_for_genni = p**2 + (q + r * padded_flag)**2
14
value_for_sbg   = p**3 + (q + r * padded_flag)**3
15

16
x0 = randbelow(N)
17
x1 = randbelow(N)
18

19
print(f'{N = }')
20
print(f'{x0 = }')
21
print(f'{x1 = }')
22
print('\nDo you prefer squares or cubes? Choose wisely!')
23

24
# Generate a random k and send v := (x_i + k^e), for Oblivious Transfer
25
# This will allow you to calculate either Genni's or SBG's value
26
# I have no way of knowing who you support. Your secret is completely safe!
27
v = int(input('Send me v: '))
28

29
m0 = (pow(v - x0, d, N) + value_for_genni) % N
30
m1 = (pow(v - x1, d, N) + value_for_sbg) % N
31
print(f'{m0 = }')
32
print(f'{m1 = }')

First attempt#

It is an OT challenge with two complicated values. So the first question is how to get the flag if I have these two values.

My method is to calculate $(p^2+C^2)^3-(p^3+C^3)^2=p(...)$ and gcd with $N$ to get $p$ . But here’s the problem: I can’t eliminate $v-x0$ and $v-x1$ when trying to get this expression. The only thing I can control is relations like (v-x0)=-(v-x1) and get the sum of two equations.

But it suggests that $p$ maybe the first target of the challenge, so instead of $\mod{N}$ , we may consider it as $\mod{p}$ .

This is a small trick to make the analysis easier. Now these values can be considered univarate with $C=q+r\times\mathrm{padded\_flag}$ .

Step 1#

Then I find this blog(I was finding scripts for multivarate coppersmith at that time and somehow he wrote a challenge about OT) and the unintended solution seems useful for this too.

If we write this out, it will look like this,

(M_1-p^3-C^3)^N\equiv x0-x1\pmod{N}

However, $e$ is change to $N$ , thus Half GCD won’t work.

Luckily, we also have another equation $p^2+C^2\equiv M_0\pmod{N}$ , which can be used to do something like pow(M_1-p^3-C^3, N, p^2+C^2-M_0). Remember that our target is to get p, so we are actually running the pow under mod $p$ , and it looks like

\begin{align*} C^2-M_0&\equiv 0\pmod{p}\\ (M_1-C^3)&\equiv x0-x1\pmod{p} \end{align*}

So $C$ is solvable and gcd(N, C^2-M_0) gives us $p$ .

After extracting $p$ , we can rerun the power but this time under mod $N$ with the true $M_1-P^3$ value. This time we get $C\pmod{N}$ , considering that $C\approx2^{1536}$ , we can throw mod $N$ and get $C$ .

So the first step ends with know the value of $p$ and $C$ . Now we can’t get more from the OT and we have to solve $q,r$ using $C, N$ .

Step 2#

Then I stucked for 3 hours and didn’t solve it before the end of contest, that’s why I claim solving 0.5 challenges. I was hoping the coppersmith would work and asked Neobeo whether the length of flag matters. He says 7 is enough, but he uses a different coppersmith. So the rest of the challenge please check his writeup.

Here’s the code of the first step,

1
from pwn import *
2
from Crypto.Util.number import bytes_to_long, getPrime
3
import sys
4
sys.setrecursionlimit(10**6)
5

6
def flatter(M):
7
    from subprocess import check_output
8
    from re import findall
9
    z = "[[" + "]\n[".join(" ".join(map(str, row)) for row in M) + "]]"
10
    ret = check_output(["flatter"], input=z.encode())
11
    return matrix(M.nrows(), M.ncols(), map(int, findall(b"-?\\d+", ret)))
12

13
def solve(N, m0, m1, x0, x1):
14
    R.<c> = PolynomialRing(Zmod(N))
15
    poly = c**2 - m0
16
    r1 = m1 - c**3
17

18
    def poly_pow(r, n, mod):
19
        if n == 0:
20
            return R(1)
21
        if n == 1:
22
            return r
23
        if n % 2 == 0:
24
            u = poly_pow(r, n // 2, mod)
25
            u = u * u
26
            _, u = u.quo_rem(mod)
27
            return u
28
        else:
29
            t = poly_pow(r, (n - 1)//2, mod)
30
            t = (t**2) * r
31
            _, t = t.quo_rem(mod)
32
            return t
33

34
    r2 = poly_pow(r1, N, poly)-(x0-x1)
35
    e, f = r2.coefficients()
36
    y = -e/f
37
    p = ZZ(gcd(N, poly(c=y)))
38

39
    assert N % p == 0 and p != 1
40
    qr = N // p
41

42
    mp0 = m0 - p**2
43
    mp1 = m1 - p**3
44

45
    poly2 = c**2 - mp0
46
    r1 = mp1 - c**3
47

48
    r2 = poly_pow(r1, N, poly2)-(x0-x1)
49
    e, f = r2.coefficients()
50
    y = ZZ(-e/f)
51
    assert poly2(c=y) == 0
52
    return qr, y
53

54
def sample():
55
    io = process(["python3", "chall.py"])
56
    # io = remote("squares-vs-cubes.chals.sekai.team", 1337, ssl=True)
57
    N = int(io.recvline().strip().split(b" = ")[1])
58
    x0 = int(io.recvline().strip().split(b" = ")[1])
59
    x1 = int(io.recvline().strip().split(b" = ")[1])
60
    io.recvuntil(b"Send me v: ")
61
    io.sendline(str(x0).encode())
62
    m0 = int(io.recvline().strip().split(b" = ")[1])
63
    m1 = int(io.recvline().strip().split(b" = ")[1])
64
    io.close()
65
    return solve(N, m0, m1, x0, x1)

zerodaycrypto#

It will be my favourite crypto of this year. ~~Though I would say it is a reverse engineering after releasing the hint.~~

1
assert __import__('re').fullmatch(rb'SEKAI{\w{32}}', flag:=input().encode()) and [pow(int.from_bytes(flag[6:-1], 'big') + i, -1, 2**255-19) >> 170 for i in range(1+3+3+7)] == [29431621415867921698671444, 12257315102018176664717361, 6905311467813097279935853, 13222913226749606936127836, 25445478808277291772285314, 9467767007308649046406595, 33796240042739223741879633, 520979008712937962579001, 31472015353079479796110447, 38623718328689304853037278, 17149222936996610613276307, 21988007084256753502814588, 11696872772917077079195865, 6767350497471479755850094]

Some Basic Analysis#

The challenge gives us high bits of 14 distinct $(x+i)^{-1}\pmod{p}$ and asks to recover the $x$ . The equations are as follows,

(x+i)a_i-1\equiv0\pmod{p}

We know the high bits of $a_i$ , and it will later be replaced from $a_i$ to $N_i+a_i$ for simplicity. The problem looks like a multivarate coppersmith, but it has fewer known bits than standard multivarate coppersmith required. Therefore, we need to find better low degree polynomials.

Recover the Lattice#

Now let’s reverse engineering the lattice. The hint decribes a linear space called SBG, and one of the important property of SBG is that it is sqaure-free for all variables. It also tells that the intended solution involves a 232-dimensional lattice reduction.

Surprisingly, the paper also claims that $SBG_{10}$ has exact 232 rank, which is the same as lattice. So we can assume that the row of lattice represents the basis decomposition of the polynomials and the rest steps are the same as the coppersmith method.

The next question is what are the polynomials of the lattice.

Recall that the key of the coppersmith is using product of polynomials to construct a multiple of $p^k$ . Following SBG’s idea, it should be sqaure-free, thus here it will be something like,

((x+i)a_i-1)((x+j)a_j-1)\equiv0\pmod{p^2}

and it shouldn’t contain $x$ , so we need to eliminate $x$ using $2d-2$ variables because SBG says that $k\geq 2d-2$ . For this case, we need $4$ variables to get a degree $3$ symmetric polynomial that is a multiple of $p^2$ .

I wrote a function F to interpolate the desired beetle, and it turns out it’s a sum of $(4,3)$ beetle and a $(4,2)$ beetle. Dig deeper we can find it is always a sum or subtraction of a $(2d-2,d)$ beetle and a $(2d-2,d-1)$ beetle.

The hint also says that SBG is closed under translation, so it will be fine to replace $a_i$ to $N_i+a_i$ .

Now we’ve fully understand the detail of the lattice.

Solve Multivarate Polynomial in Real Field#

The next step is solving it in real field.

In general, we will get some polynomials and use Groebner basis to solve it. But Groebner basis is extremely slow, so we need some tricks to solve it.

When I check the lattice, I find that except the last row, the others are correct in real field. It is quite weird and When I asked Neobeo, he says that it is even true for $1500\times1500$ lattice. Whatever, just right_kernel gives me the $a_i$ and that’s solved.

1
from re import fullmatch, findall
2
import os
3
import itertools
4

5
def flatter(M):
6
    from subprocess import check_output
7
    z = "[[" + "]\n[".join(" ".join(map(str, row)) for row in M) + "]]"
8
    ret = check_output(["flatter"], input=z.encode())
9
    return matrix(M.nrows(), M.ncols(), map(int, findall(b"-?\\d+", ret)))
10

11
flag = b"SEKAI{" + b"a" * 32 + b"}"
12
assert fullmatch(rb'SEKAI{\w{32}}', flag)
13

14
# print([pow(int.from_bytes(flag[6:-1], 'big') + i, -1, 2**255-19) >> 170 for i in range(1+3+3+7)])
15
# [29431621415867921698671444, 12257315102018176664717361, 6905311467813097279935853, 13222913226749606936127836, 25445478808277291772285314, 9467767007308649046406595, 33796240042739223741879633, 520979008712937962579001, 31472015353079479796110447, 38623718328689304853037278, 17149222936996610613276307, 21988007084256753502814588, 11696872772917077079195865, 6767350497471479755850094]
16

17
q = 2**255-19
18
oc = int.from_bytes(os.urandom(32), 'big')
19

20
# oracles = [pow(oc + i, -1, q) for i in range(20)]
21

22
def build_sbg(k, d, i_list, x_list):
23
    assert k>= 2*d-2
24
    assert len(i_list) == k
25

26
    def build_row(xi, i):
27
        return [xi*(pow(i, u)) for u in range(d)] + [pow(i, u) for u in range(k-d)]
28

29
    # return matrix([build_row(x, i) for x, i in zip(x_list, i_list)])
30
    return matrix([build_row(x_list[i], i) for i in i_list])
31

32
li = [13,9,3,4,5,6,7,8,0,17]
33
# assert (build_sbg(10, 6, li, oracles).det()-build_sbg(10, 5, li, oracles).det())%(q**5) == 0
34

35
def F(x,y,z,w):
36
    li = [x,y,z,w]
37

38
    plist = []
39
    real_plist = []
40
    for a, b in [(0,1), (0,2), (0,3), (1,2), (1,3), (2,3)]:
41
        u1, u2 = set([0,1,2,3]).difference([a,b])
42
        poly1 = ((x**2*alist[a]*alist[b]+x*(li[a]*alist[a]+li[b]*alist[b]-2)))
43
        plist.append(poly1*alist[u1])
44
        plist.append(poly1*alist[u2])
45
        poly2 = ((x+li[a])*alist[a]-1)*((x+li[b])*alist[b]-1)
46
        real_plist.append(poly2*alist[u1])
47
        real_plist.append(poly2*alist[u2])
48

49
    ms = sum(p0*randint(10, 2**32) for p0 in plist).monomials()
50
    def to_coef(poly):
51
        coefs = poly.coefficients()
52
        monos = poly.monomials()
53
        t = [0] * len(ms)
54
        for i, m in enumerate(ms):
55
            if m in monos:
56
                t[i] = coefs[monos.index(m)]
57
        return t
58

59
    B = matrix([to_coef(poly) for poly in plist])
60
    v = B.left_kernel().basis()[0]
61
    poly = 0
62
    for i in range(len(v)):
63
        poly += v[i]*real_plist[i]
64
    return poly
65

66
R.<x,a0,a1,a2,a3,a4,a5,a6,a7,a8,a9,i,j,k,l> = PolynomialRing(ZZ)
67
alist = [a0,a1,a2,a3,a4,a5,a6,a7,a8,a9]
68
# highbits = [n>>170 for n in oracles]
69
highbits = [29431621415867921698671444, 12257315102018176664717361, 6905311467813097279935853, 13222913226749606936127836, 25445478808277291772285314, 9467767007308649046406595, 33796240042739223741879633, 520979008712937962579001, 31472015353079479796110447, 38623718328689304853037278, 17149222936996610613276307, 21988007084256753502814588, 11696872772917077079195865, 6767350497471479755850094]
70
true_vlist = [a+n*(2**170) for a,n in zip(alist, highbits)]
71
# redacted_vlist = {f"a{i}": n & ((2**170)-1) for i, n in enumerate(oracles[:10])}
72

73
def create_basis(n,d):
74
    if d==0:
75
        return [(set(), R(1))]
76
    if d==1:
77
        return [(set([i]), a) for i, a in enumerate(alist)]
78
    S = list(range(0, n-d+2))
79
    T = list(range(n-d+2, n))
80
    # subset of S size d
81
    basis = []
82
    for subset in itertools.combinations(S, d):
83
        A = build_sbg(2*d-2, d, list(subset)+T, alist)
84
        poly = A.det()
85
        basis.append((set(subset), poly))
86
    return basis
87

88
all_basis = []
89
for d in range(7):
90
    bs = create_basis(10, d)
91
    all_basis.extend(bs)
92
monomials = [c[0] for c in all_basis]
93
# basis_true_value = vector(ZZ, [(c[1])(**redacted_vlist) for c in all_basis])
94
vec_basis = vector(R, [c[1] for c in all_basis])
95

96
assert len(all_basis) == 232
97

98
def create_eqs(d):
99
    assert d >= 2
100
    n = 10
101
    S = list(range(0, n-d+2))
102
    T = list(range(n-d+2, n))
103
    eqs = []
104
    for subset in itertools.combinations(S, d):
105
        poly = build_sbg(2*d-2, d, list(subset)+T, true_vlist).det() - (-1)**d * build_sbg(2*d-2, d-1, list(subset)+T, true_vlist).det()
106
        # assert poly(**redacted_vlist)%(q**(d-1)) == 0
107
        eqs.append(poly*(q**(6-d)))
108
    return eqs
109

110
all_eqs = []
111
for d in range(2, 7):
112
    all_eqs += create_eqs(d)
113

114
def decompose_sbg(poly):
115
    if poly == 0:
116
        return [0] * len(monomials)
117

118
    def find_ids():
119
        deg = poly.degree()
120
        for mono in poly.monomials():
121
            if mono.degree() == deg:
122
                vs = mono.variables()
123
                ids = [alist.index(v) for v in vs]
124
                assert len(ids) == deg
125
                if set(ids) in monomials:
126
                    return set(ids)
127

128
    ids = find_ids()
129
    assert ids is not None, poly.monomials()
130
    for s, p in all_basis:
131
        if s == ids:
132
            di, rem = poly.quo_rem(p)
133
            decomposed = decompose_sbg(rem)
134
            assert di.monomials() == [1]
135
            decomposed[monomials.index(ids)] += di
136
            return decomposed
137

138
decomposed_eqs = [vector(ZZ, decompose_sbg(eq)) for eq in all_eqs]
139
M = matrix(decomposed_eqs)
140
print(len(decomposed_eqs))
141

142
# assert all(eq*basis_true_value%(q**5) == 0 for eq in decomposed_eqs)
143

144
scale_d = {
145
    0: 1,
146
    1: 2**170,
147
    2: 2**341,
148
    3: 2**512,
149
    4: 2**688,
150
    5: 2**870,
151
    6: 2**1054
152
}
153
# print([b.bit_length() for b in basis_true_value])
154
scale = [scale_d[len(m)] for m in monomials]
155
assert len(scale) == 232
156

157
M = block_matrix([[M], [(q**5)*identity_matrix(232)[:11]]])
158

159
print(M.nrows(), M.ncols())
160

161
for i, s in enumerate(scale):
162
    M.rescale_col(i, s)
163

164
M = flatter(M)
165
M = M.change_ring(QQ)
166
for i, s in enumerate(scale):
167
    M.rescale_col(i, 1/s)
168
M = M.change_ring(ZZ)
169

170
M = M[:230]
171
# assert M * basis_true_value == vector(ZZ, [0]*230)
172

173
bs = M.right_kernel().basis()
174

175
print(bs[0][:10])
176
print(bs[1][:10])
177

178
a00 = bs[0][1]
179
x00 = pow(a00 + highbits[0]*(2**170), -1, q)
180

181
print(x00)
182
print([pow(x00 + i, -1, 2**255-19) >> 170 for i in range(1+3+3+7)])