Ricerca CTF 2023 作問者Writeup [Crypto編]

この記事は、2023年4月22日に弊社が開催したRicerca CTF 2023の公式Writeupです。今回はCryptoカテゴリの問題のうち、warmupを除く問題の解法を紹介します。

配布ファイルやスクリプトは以下のGitHubリポジトリを参照してください。

https://github.com/RICSecLab/ricerca-ctf-2023-public

  また、他のジャンルのWriteupは以下の投稿から読むことができます。

• Reversing編
• Pwn編
• Web編
• Misc, Forensics編

Rotated Secret Analysis

Writeup著者：ptr-yudai

問題概要

配布ファイルは以下のようなシンプルなスクリプトと、その実行結果になります。

import os
from Crypto.Util.number import bytes_to_long, getPrime, isPrime

flag = os.environ.get("FLAG", "fakeflag").encode()

while True:
  p = getPrime(1024)
  q = (p << 512 | p >> 512) & (2**1024 - 1) # bitwise rotation (cf. <https://en.wikipedia.org/wiki/Bitwise_operation#Rotate>)
  if isPrime(q): break

n = p * q
e = 0x10001
m = bytes_to_long(flag)

c = pow(m, e, n)

print(f'{n=}')
print(f'{e=}')
print(f'{c=}')

RSA暗号になってより、平文$m$に対する暗号$c = m^e \mod n$および公開鍵$e, n$が与えられます。ただし、素数$p, q$の作り方が特徴的で、$p$は$q$を512ビットローテートした値になっています。

解法

RSAの素数$p, q$に何かしらの関係があるとき、$n$を素因数分解できてしまうことが多くあります。今回の問題では、512ビットに収まる整数$a, b$を使って、

$p = 2^512 a + b$

$q = 2^512 b + a$

と表すことができます。したがって、

$n = (2^512 a + b) (2^512 b + a) = 2^1024 ab + 2^512 (a^2 + b^2) + ab$

となります。項ごとにオーダーが異なることが分かります。図で確認してみましょう。

$ab$の下位512ビットと上位512ビットは、それぞれ$n$の下位512ビットと上位512ビット（図の赤色の部分）から計算できます。$2^512(a^2+b^2)+ab$の繰り上がりを考慮すると、$n$の上位512ビット$ab$は完全ではありません。しかし、繰り上がりは発生するかしないかの2通りなので、完全な$ab$の値を求められます。

また、$ab$が分かると$a^2+b^2$の値も

$a^2 + b^2 = 2^{-512}(n - (2^1024 + 1) ab)$

という計算で取り出せるため、$a^2+b^2$と$ab$が求まります。さらに、

$(a+b)^2 = (a^2 + b^2) + 2(ab)$

より、$a+b$も求まり、ここから$a, b$を計算できます。

したがって、$n$から$a, b$の両方の値を取り出せるため、$p, q$が分かります。

ソルバ

from Crypto.Util.number import long_to_bytes
from math import isqrt

exec(open("../distfiles/output.txt").read())

P_BITS = 1024
A_BITS = P_BITS // 2

ab_upper = int(bin(n)[2:][:-A_BITS*3], 2)
ab_lower = int(bin(n)[2:][-A_BITS:], 2)

a, b = None, None
for i in [0, 1]:
  ab = ((ab_upper - i) << A_BITS) + ab_lower
  # a^2+b^2
  a2pb2 = (n - (ab << P_BITS) - ab) >> A_BITS
  # (a+b)^2
  apb2 = a2pb2 + 2 * ab
  if apb2 < 0 or isqrt(apb2) ** 2 != apb2: continue

  apb = isqrt(a2pb2 + 2 * ab)
  d = apb**2 - 4*ab
  if d < 0 or isqrt(d) ** 2 != d: continue

  a = (apb + isqrt(d)) // 2
  b = (apb - isqrt(d)) // 2

  assert a * b == ab
  break

assert a is not None and b is not None

p = (a << A_BITS) + b
q = (b << A_BITS) + a

assert p * q == n

print(long_to_bytes(pow(c, pow(e, -1, (p - 1) * (q - 1)), n)).decode())

 

RSALCG

問題・Writeup著者：ptr-yudai

問題概要

問題名の通り、RSA暗号とLCGを組み合わせた暗号化を使っています。 LCGは疑似乱数を生成するスキームで、定数$a, b, n$とシード$s$に対して、乱数列$\{r_i\}$を

$$r_0 = s \ \ \ \ (s < n) \\ r_i \equiv ar_{i-1} + b \mod n \ \ \ \ (i \in \mathbb{N})$$

と定義します。

また、RSAでは素数$p,q$を使い$n = pq$とし、$e$を65537と定義しています。

RSALCGでは、平文$m$を$m \oplus (r_i^e \mod n)$により暗号化します。 すなわち、LCGが生成した乱数$r_i$をRSAで暗号化し、その結果と平文をXORすることで最終的な暗号文を生成します。

公開鍵は$n, e, a, b$で、秘密鍵は$s$です。 また、1回目に既知平文、2回目にフラグ、3回目に別の既知平文を暗号化した結果が与えられます。

解法

Franklin-Reiter Related Message Attack

今回は3回暗号化するため、$r_1, r_2, r_3$が使われます。

$r_1 \equiv as + b \mod n$

$r_2 \equiv a^2s + ab + b \mod n$

$r_3 \equiv a^3s + a^2b + ab + b \mod n$

また、1回目と3回目に暗号化される平文は分かっているので、XORを戻して

$c_1 \equiv r_1^e \equiv (as + b)^e \mod n$

$c_3 \equiv r_3^e \equiv (a^3s + a^2b + ab + b)^e \mod n$

は判明しています。

ここで、

$r_3 \equiv a^3s + a^2b + ab + b = a^2(as + b) + (ab + b) = a^2r_1 + (ab + b)$

です。 $A=a^2, B=ab+b, r_1=m$とおくと、

$r_3 \equiv Am + B  \mod n \ \ \ \ (B \neq 0)$

の形になります。$m, Am + B$の暗号文を持っているとき、Franklin-Reiter Related Message Attackが適用できます。

剰余多項式$f_1(x), f_3(x)$を

$f_1(x) \equiv c_1 - x^e \mod n$

$f_3(x) \equiv c_3 - (Ax+B)^e \mod n$

とおくと、$f_1(x), f_3(x)$はともに$m$を解に持つため、

$f_1(x) = (x - m)h_1(x)$

$f_3(x) = (x - m)h_3(x)$

と因数分解できるはずです。 したがって、$\gcd(f_1(x), f_3(x)) = x - m$より、平文$m$が求まります。

Half-GCD

$f_1(x), f_3(x)$のGCDを計算する際に問題が発生します。 今回$e=65537$ですので、$f_1(x), f_3(x)$はともに65537を次数として持ちます。このような非常に大きい次数の多項式では、ユークリッドの互除法を使っても（現実的に終わるものの）計算に時間がかかります。

次数が高い多項式のGCDを取るとき、より高速なアルゴリズムとしてHalf-GCDがあります。 今回のパターンでは、ユークリッドの互除法では数十分から数時間かかる計算が、Half-GCDでは数十秒から数分で終わります。

ソルバ

以下は、Franklin-Reiter Related Message AttackとHalf-GCDを組み合わせたソルバになります。

# Half-GCD code from
# <https://scrapbox.io/crypto-writeup-public/half-GCD>
def pdivmod(u, v):
    """
    polynomial version of divmod
    """
    q = u // v
    r = u - q*v
    return (q, r)

def hgcd(u, v, min_degree=10):
    """
    Calculate Half-GCD of (u, v)
    f and g are univariate polynomial
    <http://web.cs.iastate.edu/~cs577/handouts/polydivide.pdf>
    """
    x = u.parent().gen()

    if u.degree() < v.degree():
        u, v = v, u

    if 2*v.degree() < u.degree() or u.degree() < min_degree:
        q = u // v
        return matrix([[1, -q], [0, 1]])

    m = u.degree() // 2
    b0, c0 = pdivmod(u, x^m)
    b1, c1 = pdivmod(v, x^m)

    R = hgcd(b0, b1)
    DE = R * matrix([[u], [v]])
    d, e = DE[0,0], DE[1,0]
    q, f = pdivmod(d, e)

    g0 = e // x^(m//2)
    g1 = f // x^(m//2)

    S = hgcd(g0, g1)
    return S * matrix([[0, 1], [1, -q]]) * R

def pgcd(u, v):
    """
    fast implementation of polynomial GCD
    using hgcd
    """
    if u.degree() < v.degree():
        u, v = v, u

    if v == 0:
        return u

    if u % v == 0:
        return v

    if u.degree() < 10:
        while v != 0:
            u, v = v, u % v
        return u

    R = hgcd(u, v)
    B = R * matrix([[u], [v]])
    b0, b1 = B[0,0], B[1,0]
    r = b0 % b1
    if r == 0:
        return b1

    return pgcd(b1, r)

with open("../distfiles/output.txt") as f:
    exec(f.readline())
    exec(f.readline())
    exec(f.readline())
    C1 = bytes.fromhex(f.readline())
    C2 = bytes.fromhex(f.readline())
    C3 = bytes.fromhex(f.readline())

# s1 = a*s + b mod n
# s2 = a^2*s + a*b + b mod n
# s3 = a^3*s + a^2*b + a*b + b mod n

# f1 = s1^e - c1 mod n
# f2 = s3^e - c2 mod n

M1 = b"The quick brown fox jumps over the lazy dog"
M3 = b"<https://translate.google.com/?sl=it&tl=en&text=ricerca>"
c1 = int.from_bytes(C1, 'big') ^^ int.from_bytes(M1, 'big')
c3 = int.from_bytes(C3, 'big') ^^ int.from_bytes(M3, 'big')

print("[+] Creating polynomial...")
F = Zmod(n)
PR.<x> = PolynomialRing(F)
e = 65537

f = (a*x + b)^e - c1
g = (a^3*x + a^2*b + a*b + b)^e - c3

print("[+] Calculating half-GCC...")
h = pgcd(f, g)
s0 = F(-h.monic()[0])

s1 = (a * s0 + b) % n
s2 = (a * s1 + b) % n
key = pow(s2, e, n)
m = int(key) ^^ int.from_bytes(C2, 'big')

print(int.to_bytes(m, 128, 'big').lstrip(b'\x00'))

 

dice-vs-kymn

問題著者：ptr-yudai, keymoon / Writeup著者：ptr-yudai

問題概要

この問題は、keymoonの振ったサイコロの目を77回連続で当てることでフラグを得られます。

サイコロの目は次の手順で生成されます。

1. 64ビット程度の乱数$x$を生成する。
2. ユーザーが整数$a (\neq 0), b (\neq 0)$を設定する。
3. $x$が楕円曲線$y^2 = x^3 + ax + b \mod p$上の点として有効になるように、256ビット程度のランダムな素数$p$を生成する。
4. $G=(x, y)$（yはlift_xにより決定されるどちらかの座標）とする。
5. 乱数$k$で$Q=kG$を計算し、$G_y, Q_y$を開示する。
6. $k + 1 \mod p$がkeymoon側のサイコロの目
7. ユーザーが自由にサイコロの目を設定する。

keymoonのサイコロの目と自分のサイコロの目の和が7になれば次のラウンドに進めます。つまり、毎ラウンドで$k$の値を当てる問題になります。もちろん、この問題は楕円曲線上の離散対数問題である上、$p$の値がわからず楕円曲線が構築できないため解けません。

解法

等分多項式の計算

$Q=kG$の$k \mod 6$の値がわかれば、サイコロのどの目が出るかがわかります。 したがって、ベースポイント$G$が$7G = G$と巡回するように楕円曲線のパラメータ$a, b$を設定する必要があります。今回の場合$6G=\mathcal{O}$とも表せます。

では、$nP = \mathcal{O}$を満たすためにはどうすれば良いでしょうか。 このような点は$n$-ねじれ点（$n$-torsion point）と呼ばれ、楕円曲線の加法公式から$n$-ねじれ点の$x$座標を根に持つ多項式が計算できます。また、この多項式を等分多項式と呼びます。

SageMathの場合、楕円曲線Eに対してE.division_polynomial(n)として等分多項式$f_{n}(x) \in \mathbb{F}_{p}[x]$を求められますが、今回は$a, b$が未定で$p$が未知なので、このメソッドは使えません。

代わりに、$a,b$を変数とした整数等分多項式$f_{6}(a,b) \in \mathbb{Z}[a,b]$を計算します。$\mathbb{Z}$上で等分多項式の根が求められれば、その値は$\mathbb{F}_p$上でも等分多項式の根になっているはずです。

等分多項式を$\mathbb{Z}$上で実装してみましょう。

from functools import cache

@cache
def division_polynomial(a, b, x, m):
    f = lambda m: division_polynomial(a, b, x, m)
    if m == 0:
        return 0
    elif m == 1 or m == 2:
        return 1
    elif m == 3:
        return 3*x^4 + 6*a*x^2 + 12*b*x - a^2
    elif m == 4:
        return 2*(x^6 + 5*a*x^4 + 20*b*x^3 - 5*a^2*x^2 - 4*a*b*x - 8*b^2 - a^3)
    elif m % 2 == 0:
        n = m // 2
        return (f(n+2)*f(n-1)^2 - f(n-2)*f(n+1)^2) * f(n)
    else:
        n = (m-1) // 2
        F = 4*(x^3 + a*x + b)
        if n % 2 == 1:
            return f(n+2)*f(n)^3 - F^2*f(n-1)*f(n+1)^3
        else:
            return F^2*f(n+2)*f(n)^3 - f(n-1)*f(n+1)^3

PR.<a, b, x> = ZZ[]
f6 = division_polynomial(a, b, x, 6)
print(f6)

実行すると、次のような多項式になります。

 与えられた$x$を代入して生成された多項式の根$a, b \in \mathbb{Z}$を求めるのが問題です。 上記の多項式をfactor関数で因数分解すると、次のようになります。

 もっとも簡単な式は

$h_x(a, b) = 3x^4 + 6ax^2 - a^2 + 12bx$

ですが、これは3等分多項式なので使えません。したがって、

$h_x(a, b) = x^{12} + 22ax^{10} - 165a^2x^8 + 220bx^9 - 92a^3x^6 - 528abx^7 - 185a^4x^4 + 264a^2bx^5 - 1776b^2x^6 - 90a^5x^2 - 80a^3bx^3 - 960ab^2x^4 - 3a^6 - 132a^4bx - 624a^2b^2x^2 - 320b^3x^3 - 96a^3b^2 - 896ab^3x - 512b^4$

を採用します。 $h_x(a,b)$が恒等的に0になるような$a, b$を求める必要がありますが、$a, b$は$x$に依存するためこのような$a, b$を求めるのは困難です。

見通しをよくするために、まずは3等分、4等分について考えてみます。

$f_3(x,a,b) = 3x^4 + 6ax^2 - a^2 + 12bx$ $h_4(x,a,b) = 2x^6 + 10ax^4 - 10a^2x^2 + 40bx^3 - 2a^3 - 8abx - 16b^2$

いま$a, b$は$x$の関数ですが、楕円曲線の式より$a,b$の次数をそれぞれ2次,3次に限定してみましょう。つまり、$a(x)=a'x^2, b(x)=b'x^3$と置きます。これにより等分多項式を$x$の累乗と$x$を含まない式に因数分解できるので、$a,b$に関する$x$と独立した問題として扱えます。

$f_3(a,b)=x^4(-a^2 + 6a + 12b + 3)$

$f_4(a,b)=x^6(-2a^3 - 10a^2 - 8ab - 16b^2 + 10a + 40b + 2)$

これらの式はいずれも整数解を持ちます。そこで、先程採用した式$h_x$でも同様に$a,b$を設定してみましょう。このとき解くべき式は次のような形をしています。

$h_x(a,b)=x^{12}(-3a^6 - 90a^5 - 132a^4b +\cdots)$

8次式は解の公式を持ちません。係数のオーダーが小さいことから解も小さいことに期待して、解を探索してみましょう。

PR.<a, b, x> = ZZ[]
f3 = division_polynomial(a, b, x, 3)
f6 = division_polynomial(a, b, x, 6)
ff = (f6 / f3)(x=1)

for aa in range(0x1000):
    for bb in range(0x1000):
        if ff(a=aa, b=bb) == 0:
            print(f"Found: (a, b) = ({aa}, {bb})")
        elif ff(a=aa, b=-bb) == 0:
            print(f"Found: (a, b) = ({aa}, {-bb})")
        elif ff(a=-aa, b=bb) == 0:
            print(f"Found: (a, b) = ({-aa}, {bb})")
        elif ff(a=-aa, b=-bb) == 0:
            print(f"Found: (a, b) = ({-aa}, {-bb})")

すぐに解が見つかります。

$ sage test.sage
Found: (a, b) = (-3, 2)
Found: (a, b) = (-15, -22)

確認してみましょう。

print(f3(a=-3*x^2, b=2*x^3))
print(f6(a=-3*x^2, b=2*x^3))
print(f3(a=-15*x^2, b=-22*x^3))
print(f6(a=-15*x^2, b=-22*x^3))

探索で見つかった(-3, 2)は3等分多項式の解にもなってしまっているので使えません。(-15, -22)は、6等分多項式の解ですが3等分多項式の解にはなっていないため、理想的です。

0
0
-576*x^4
0

したがって、$a=-15x^2, b=-22x^3$と置くと、与えられた$x$が$p$に関わらず6等分点を生成することが分かりました。

2G,3G,4G,5Gの特性

等分多項式を満たす$a,b$を与えれば、$Q=kG$の$G$と$Q_y$をサーバーから貰えることになります。

$k\equiv 1 \mod 6$のとき、$Q=G$より、$Q_y=G_y$を確認することで判断できます。 $k\equiv 0 \mod 6$のとき、$Q=\mathcal{O}$より、$Q_y=1$を確認することで判断できます。（無限遠点$Q=(0:1:0)$に対してSageMathはQ[1]と記述すると$Z$に関係なく$Y$の値を返す。）

では、$k \mod 6$が2,3,4,5のときはどのように判定すれば良いでしょうか。

6-ねじれ点に対して、今回の楕円曲線上で$y$座標がどのような特性を持つか調べるため、シンボリックに$2G$から$5G$を計算してみましょう。 剰余$p$がわかっていない状態で楕円曲線の演算をしなくてはならないので、ここでは有理環$\mathbb{Q}$上で点$G(x,y)$がどこに移るかを調べます。

PR.<x, y> = PolynomialRing(QQ)
a = -15*x^2
b = -22*x^3

def reduce_zz(f):
    poly = 0
    for coeff, biv in list(f):
        while biv % y^2 == 0:
            biv /= y^2
            biv = biv.numerator()
            biv *= x^3 + a*x + b
        poly += coeff * biv
    return poly

def reduce_qq(f):
    return reduce_zz(f.numerator()) / reduce_zz(f.denominator())

x1, y1 = x, y

x2 = ((3*x1^2 + a) / (2*y1))^2 - 2*x1
y2 = (3*x1^2 + a) / (2*y1) * (x1 - x2) - y1
x2, y2 = reduce_qq(x2), reduce_qq(y2)
print("[ 2G ]")
print("x:", x2)
print("y:", y2)

for k in range(3, 6):
    x1, y1 = x2, y2
    x2 = ((y1 - y) / (x1 - x))^2 - x1 - x
    y2 = ((y1 - y) / (x1 - x)) * (x1 - x2) - y1
    x2, y2 = reduce_qq(x2), reduce_qq(y2)
    print(f"[ {k}G ]")
    print("x:", x2)
    print("y:", y2)

実行結果は次のようになります。

[ 2G ]
x: -3*x
y: (-3456*x^3)/(-288*y)
[ 3G ]
x: -2*x
y: 0
[ 4G ]
x: -3*x
y: 1/3*y
[ 5G ]
x: x
y: -y

つまり、$\mathbb{F}_p$の$p$にかかわらず、次のように点の座標が遷移することがわかります。

$$G = (x, y)\\ 2G = (-3x, 12x^3/y)\\ 3G = (-2x, 0)\\ 4G = (-3x, y/3)\\ 5G = (x, -y)$$

サーバーからは$y$座標が与えられるので、ひとまず0かどうかを調べれば$3G$は判定できます。

2G,4G,5Gの判別

ここまでの調査から、$2G$の$y$座標は$12x^3/y$、$4G$の$y$座標は$y/3$、$5G$の$y$座標は$-y$になることがわかりました。これは$\mathbb{Q}[x,y]$上での計算ですが、$\mathbb{F}_p[x,y]$にそのまま移しても除算を逆元の積を見れば同一視できます。

与えられた$Q_y$が$2G$のものか$4G$のものかを判定できれば、どれにも当てはまらないものが$5G$なので、まずは$2G$と$4G$の判定を考えます。

いま持っている情報は、

$x_1 = G_x, y_1 = G_y$と、$y_2=(2G)_y$、$y_4=(4G)_y$あるいは$y_5=(5G)_y$です。

まず、楕円曲線の公式から

$y_1^2 \equiv x_1^3 + ax_1 + b \equiv x_1^3 - 15x_1^3 - 22x_1^3 \equiv -36x_1^3 \mod p$

です。したがって、ある整数$k_1$について

$y_1^2 + 36x_1^3 = k_1p$

と表せます。...① 次に、サイコロの目が2の場合、すなわち$y_2$を持っている場合を考えます。このとき、

$y_2 \equiv \cfrac{12x_1^3}{y_1} \mod p$

より、ある整数$k_2$について

$y_2y_1 - 12x_1^3 = k_2p$

と表せます。...②

同様に、$y_4$を持っている場合、整数$k_4$について

$3y_4 - y_1 = k_3p$

が成り立ちます。...③

①②③より、どちらのものか不明な$y$座標$y_k$について

$s = \gcd(y_1^2 + 36x_1^3, y_ky_1 - 12x_1^3)$

$t = \gcd(y_1^2 + 36x_1^3, 3y_k - y_1)$

を計算すれば、$y_k=y_2$なら$p|s$、$y_k=y_4$なら$p|t$が成り立ちます。

かならずしも$s, t$が現実的な時間で素因数分解できるとは限りませんが、$p$を素因数に含む場合これらの値は十分に大きくなるはずです。 したがって、$s, t$を計算して十分に大きい値が出れば、それぞれの場合で$2G, 4G$であったことが分かり、どちらも小さい値になった場合は$5G$であったことが分かります。

ソルバ

ここまでの議論を整理すると、次のようにほぼ100%でサイコロの出目を判定できます。

from ptrlib import *
import os
from tqdm import tqdm

HOST = os.getenv("HOST", "localhost")
PORT = int(os.getenv("PORT", 5963))

# sock = Process(["sage", "../distfiles/server.sage"])
sock = Socket(HOST, PORT)

for round in tqdm(range(77)):
    Gx = int(sock.recvlineafter("x: "))
    a = -15*Gx**2
    b = -22*Gx**3
    sock.sendlineafter("a: ", a)
    sock.sendlineafter("b: ", b)

    commitment = sock.recvregex("\((\d+), (\d+)\)")
    Gy, Qy = int(commitment[0]), int(commitment[1])
    print(f'{round=} {Gy=}, {Qy=}')
    k1p = Gy**2 - Gx**3 - a*Gx - b

    if Qy == 0: # 3*G
        kymn_dice = 4

    elif Qy == 1: # 6*G
        kymn_dice = 1

    elif Qy == Gy: # 1*G
        kymn_dice = 2

    else:
        k2p = Qy * (-288*Gy) + 3456*Gx**3
        p = gcd(k1p, k2p)
        if p > 1<<200: # 2*G
            kymn_dice = 3

        else:
            k2p = Qy * 3 - Gy
            p = gcd(k1p, k2p)
            if p > 1<<200: # 4*G
                kymn_dice = 5

            else:
                kymn_dice = 6

    sock.sendlineafter("Your dice: ", 7 - kymn_dice)

    l = sock.recvline()
    if b'You lucky bastard' not in l:
        logger.error("Bad luck!")
        sock.close()
        break

print(sock.recvuntil(b'}'))