DownUnderCTF 2021 | yadlp

#joseph

#downunderctf2021

def G_add(A, B):
    x1, y1 = A
    x2, y2 = B
    return ((x1*x2 + D*y1*y2) % p, (x1*y2 + x2*y1 + 2*y1*y2) % p)

def G_mul(A, k):
    out = (1, 0)
    while k > 0:
        if k & 1:
            out = G_add(out, A)
        A = G_add(A, A)
        k >>= 1
    return out

def rand_element():
    while True:
        x = randint(1, p-1)
        d = x^2 * (D + 1) - D
        if (x & 1 == d & 1) and kronecker(d, p) == 1:
            y = (x + sqrt(Zmod(p)(d))) * inverse_mod(D, p) % p
            return (x, y)

D = 13337
p = 17568142778435152362975498611159042138909402642078949814477371651322179417849164549408357464774644525711780515232117470272550677945089719112177956836141583
assert p.nbits() >= 512
assert ((p-1)//2).is_prime() # safe prime

FLAG = open('flag.txt', 'rb').read().strip()
assert len(FLAG) % 8 == 0
M = [int.from_bytes(FLAG[i:i+8], 'big') for i in range(0, len(FLAG), 8)]

G = [rand_element() for _ in M]
c = (1, 0)
for m, gi in zip(M, G):
    c = G_add(c, G_mul(gi, m))

print(f'{D = }')
print(f'{p = }')
print(f'{G = }')
print(f'{c = }')

 x^2 + 2xy - Dy^2 \equiv 1 \mod pという曲線上のdlpができますか、という問題。

 c = m_0G_0 + m_1G_1 + \dots + m_kG_k m_iを求めるのが問題だけど、適当なbasepoint  Gを考えると c = m_0x_0G + m_1x_1G + \dots + m_kx_kGとなるのでdlpしてからLLL

あるいはたくさん G, G', G''', \dotsを作って行列を解いても良い

方針1

少し変形すると (x+y)^2 - (D+1)y^2 \equiv 1 \mod pという形になるのでさらにもう少しずらして

 x'^2 - D'y^2 \equiv 1 \mod p とするとこれはPell's Equation

Pell's Equationでは次の性質の良い式が成り立つ

 (a^2 - Db^2)(c^2 - Dd^2) = (ac + Dbd)^2 - D(ad + bc)^2

ここから加法公式が導かれる

わからん……

 \mathbb{F_{p^2}} \simeq \mathbb{F_p}\lbrack x \rbrack (x^2 - D)で、 (a, b) \to a+bxという関係が成り立つらしい

http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.66.8688&rep=rep1&type=pdf

方針2

雑にやると p+1が位数っぽいことがわかる(適当な点 Pをとって実験してみると (p+1)P = \mathcal{O}になるので)

 p + 1は小さい素因数をたくさん持つので、Pohlig-Hellman Attacksageのdiscrete logが対応していない群の離散対数問題を解くテクをやると解ける

def G_inv(A):
    return G_mul(A, p)


def bsgs(g, h, n):
    """return x s.t. g^x = h"""
    m = ceil(sqrt(n))
    table = {}
    for j in range(m):
        table[G_mul(g, j)] = j
    factor = G_mul(G_inv(g), m)
    gamma = h
    for i in range(m):
        if gamma in table:
            j = table[gamma]
            ret = i * m + j
            assert G_mul(g, ret) == h
            return ret
        gamma = G_add(gamma, factor)

O = (1, 0)

def calc_order(g):
    order = p + 1
    for pi, e in list(factor(order)):
        for i in range(1, e+1):
            if G_mul(g, order // pi) == O:
                order //= pi
    return order


def pohlig_hellman(g, h):
    a_list = []
    b_list = []
    order = calc_order(g)
    for pi, e in list(factor(order)):
        gi = G_mul(g, order // pi ** e)
        hi = G_mul(h, order // pi ** e)
        gamma = G_mul(gi, pi ** (e - 1))
        xk = 0
        for k in range(e):
            hk = G_mul(G_add(G_mul(G_inv(gi), xk), hi), (pi ** (e - 1 - k)))
            dk = bsgs(gamma, hk, pi)
            xk = xk + pi ** k * dk
        xi = xk
        a_list.append(xi)
        b_list.append(pi ^ e)
    return crt(a_list, b_list)


ds = [pohlig_hellman(c, gi) for gi in G]

# [q, m1, m2, ..., m6]
mat = matrix(ZZ, 7, 7)
mat[0, 0] = -(p+1)
for i in range(1, 7):
    mat[i, 0] = ds[i-1]
    mat[i, i] = 1
C = mat.LLL()
ms = mat.solve_left(C[0])
b"".join(long_to_bytes(m) for m in ms)