SECCON 14 quals: yukari infinity Writeup

内容の間違いがないように確認しながら作っていますが、専門外のためもし間違い等ありましたらコメント等で指摘していただけますと幸いです。

yukari infinity

authored by chocorusk

everlasting fuwafuwa

#!/usr/bin/env python3

from Crypto.PublicKey import RSA
from Crypto.Util.number import getPrime, isPrime

with open("flag.txt", "r") as f:
    FLAG = f.read()

for _ in range(128):
    p = getPrime(1024)
    print("p =", p)

    q = int(input("q: "))
    assert p != q
    assert q.bit_length() == 1024
    assert isPrime(q)

    n = p * q
    e = getPrime(64)
    d = pow(e, -1, (p - 1) * (q - 1))

    try:
        cipher = RSA.construct((n, e, d))
    except:
        print("error!")
        continue
    print("key setup successful")
    exit()

print(FLAG)

以前投稿したyukariの続編です。
https://qiita.com/potyamaaaa/items/406d5344e2549a517d3c
仕組みとして変わっているのは、

assert q.bit_length() == 1024

と素数$q$が1024ビットのみを受け付けるところです。
さて、このような厳しい条件下で128回も突破できるのでしょうか？

解法

本問題の着眼点

問題は次のように動作しています。

1. 1024bitの素数$p$を出力
2. 入力として$p \neq q$かつq.bit_length() == 1024を満たす$q$を与える
3. $n = pq, e = \text{getPrime(64)}, d \equiv e^{-1}\pmod{(p-1)(q-1)}$に従ってRSAパラメータ生成
4. RSA.construct()に$n,e,d$を渡し、例外が出れば成功

yukariのように、$q = 2kp + 1 \pmod p$としてしまうと、1024ビットの制限に引っかかってしまい、大変厳しいです。
ということで、再度pycryptodomeの実装を見てみると、RSA.constructが$n, e, d$から$p,q$の素因数分解を内部で試している処理が気になります。
実装（抜粋）は以下の通りで、この処理はRabin (1979) による非決定的な因数分解手法に基づいています。
論文: "Digitalized signatures and public-key functions as intractable as factorization"
コード: https://github.com/Legrandin/pycryptodome/blob/master/lib/Crypto/PublicKey/RSA.py#L539

        else:
            # Compute factors p and q from the private exponent d.
            # We assume that n has no more than two factors.
            # See 8.2.2(i) in Handbook of Applied Cryptography.
            ktot = d * e - 1
            # The quantity d*e-1 is a multiple of phi(n), even,
            # and can be represented as t*2^s.
            t = ktot
            while t % 2 == 0:
                t //= 2
            # Cycle through all multiplicative inverses in Zn.
            # The algorithm is non-deterministic, but there is a 50% chance
            # any candidate a leads to successful factoring.
            # See "Digitalized Signatures and Public Key Functions as Intractable
            # as Factorization", M. Rabin, 1979
            spotted = False
            a = Integer(2)
            while not spotted and a < 100:
                k = Integer(t)
                # Cycle through all values a^{t*2^i}=a^k
                while k < ktot:
                    cand = pow(a, k, n)
                    # Check if a^k is a non-trivial root of unity (mod n)
                    if cand != 1 and cand != (n - 1) and pow(cand, 2, n) == 1:
                        # We have found a number such that (cand-1)(cand+1)=0 (mod n).
                        # Either of the terms divides n.
                        p = Integer(n).gcd(cand + 1)
                        spotted = True
                        break
                    k *= 2
                # This value was not any good... let's try another!
                a += 2
            if not spotted:
                raise ValueError("Unable to compute factors p and q from exponent d.")
            # Found !
            assert ((n % p) == 0)
            q = n // p

この部分ではまず、
$$
k_{tot} = de-1
$$
を計算します。さて、$k_{tot}$は$\phi(n)$の倍数かつ偶数なので、
$$
k_{tot} = t \cdot 2^s
$$
と表せます(ただし、$t$は奇数です)。

次に、候補各$a(a=2,4,\dots,98)$ごとに
$$
cand_i \equiv a^{t\cdot2^i}
$$
を順に計算します。もし、
$$
cand_i^2 \equiv 1\pmod n, cand_i \not\equiv \pm 1\pmod n
$$
を満たす非自明な平方根を得るなら、
$$
(cand_i - 1)(cand_i+1) \equiv 0 \pmod n
$$
より、$\text{gcd}(cand_i \pm 1, n)$によって非自明な因子が得られるため、因数分解が成功します。

従って、この実装の試行範囲に対して因数分解を回避するには、各$a$と各$i$について
$$
cand_i^2 \equiv 1\pmod n, cand_i \equiv \pm 1\pmod n
$$
と非自明な平方根が現れないようにできればよいです。

ここで、$cand_i \equiv \pm 1\pmod n$が起きるためには、CRTより
$$
a^{t\cdot2^i} \equiv \pm1 \pmod p, a^{t\cdot2^i} \equiv \pm1\pmod q
$$
を同時に満たす必要があります。さらに、それらの符号はどちらも同符号である必要があります。

例えば、$p=3, q=5$という簡単なパターンを考えてみましょう。
同符号の場合は

\begin{align}
cand \equiv 1 \pmod 3 , cand \equiv 1 \pmod 5\\
cand \equiv 2  \pmod 3 , cand \equiv 4 \pmod 5
\end{align}

のとき

\begin{align}
cand \equiv 1 \pmod {15}\\
cand \equiv 14 \equiv -1 \pmod {15}
\end{align}

から、

\begin{align}
cand^2 &= 1^2  \equiv 1 \pmod {15}\\
cand^2 &= 14^2  \equiv 1 \pmod {15}
\end{align}

と条件を満たしています。
一方で、異符号の場合は

\begin{align}
cand \equiv 1 \pmod 3 , cand \equiv 4 \pmod 5\\
cand \equiv 2  \pmod 3 , cand \equiv 1 \pmod 5
\end{align}

のとき

\begin{align}
cand \equiv 4 \pmod {15}\\
cand \equiv 11  \pmod {15}
\end{align}

から、

\begin{align}
cand^2 &= 4^2  \equiv 1 \pmod {15}\\
cand^2 &= 11^2  \equiv 1 \pmod {15}
\end{align}

となり、条件を満たせていません。
つまり、非自明な平方根を出さないためには、同符号である必要があります。

さて、このアルゴリズムでは、$k_{tot} = t \cdot 2^s$のうち、$t$を求めています。つまり、2べきの部分を辿っていることがわかります。
ということは、
$$
V = v_2(p-1)
$$
を計算し、この$V$の値によって上記の条件を満たすような$q$を探すのが今回の目的です。

$v_2(n)$は2-adic valuationで、$n$を割り切る2のべき乗の指数を表します

V=1のとき

$V=1$、つまり
$$
V = v_2(p-1) = 1
$$
の時を考えましょう。
この時、2-Sylow 部分群は{$-1, 1$}のみとなります。以上から、$V=1$のとき、2べきは1になるか-1になるかの二値で記述でき、0/1的な情報に落とすことができます。
特に奇数$t$に対しては、
$$
\left ( \frac{a^t}{p}\right) = \left (\frac{a}{p}\right)
$$
が成立するため、法$p$における符号は$a$の平方剰余、すなわちLegendre記号で表せます。同様に法$q$でも$\left(\frac{a}{q}\right)$で表せるので、試される各$a$について
$$
\left (\frac{a}{p}\right) = \left (\frac{a}{q}\right)
$$
となるように設計し、最後にCRTにより一意に合成することによってRSA.constructが失敗する$q$を構成することができます。

V >= 2のとき

一方で、$V\geq 2$のとき平方剰余の比較のみでは不十分です。
実際、
$$
V = v_2(p-1) \geq 2
$$
のとき、$(\mathbb{Z}/p\mathbb{Z})^*$は位数$p-1$の巡回群であるから、位数$2^2, 2^3, \cdots, 2^V$の元が存在し得ます。さて、$V=1$のときに使用したLegendre記号はaが平方か否かの2択しか区別できないため、高次の2べき情報まで揃えることはできません。
例えば、$p=17$の時を考えてみましょう。
$4, -1$はともに平方剰余ですが、
$$
ord_{17}(4) = 4, ord_{17}(-1) = 2
$$
は異なります。そして、
$$
v_2(ord_{17}(4) = 2), v_2(ord_{17}(-1)) = 1
$$
となります。実際、SageMathによって次のように確認できます。

print(kronecker(4, 17))
print(kronecker(-1, 17))
print("====")
print(mod(4,17).multiplicative_order())
print(mod(-1,17).multiplicative_order())

1
1
====
4
2

よって$V\geq2$では、Legendre記号の一致だけでは$v_2(ord_p(a))$がどの値かを揃えることができません。
以上より、$V \geq 2$の場合には、平方剰余性だけでなく、2べき剰余性まで一致させる必要があります。

この高次の2べき剰余性を扱うためには、$2^V$乗根の単位を含む円分体
$$
K=\mathbb{Q}\left(\zeta_{2^{V}}\right)
$$
を基底として用います。円分体では、素数 $r\nmid 2^V$ の分解は $r\bmod 2^V$で決まり、特に
$$
r\equiv 1\pmod{2^V}
$$
なら $r$ は $K$ で完全分解します。
よって
$$
p\equiv 1\pmod{2^V}
$$
なら、$p$ の上にノルム $p$ をもつ素イデアル $\mathfrak p\mid(p)$ が存在します。
以降は、必要に応じて$\mathfrak p$ が主イデアルになるような$p$を選び、
$$
\mathfrak p=(\pi_p)
$$
と書ける場合を考えます。主イデアルのノルムの定義より
$$
\left|N_{K/\mathbb{Q}}(\pi_p)\right|=N_{K/\mathbb{Q}}((\pi_p)) =N_{K/\mathbb{Q}}(\mathfrak p) = p
$$
が成立します。
次に、ある合同条件の下でノルムが素数となる元を得るため、$\mathcal O_K$ の元 $t\in\mathcal O_K$ を用いて
$$
\pi=\pi_p+M t
$$
と動かし、
$$
q=\left|N_{K/\mathbb{Q}}(\pi)\right|
$$
が素数となる$\pi$を探索します。
こうして、$\pi$を合同条件（$\bmod M$）で縛ったままノルムが素数になるような元を拾うことができます。

そして$v_2(ord_q(a))$ を $v_2(ord_p(a))$ に一致させるためには、円分体 $K$ が $2^V$ 乗根の単位を含むことを利用して Kummer 拡大

$$
L_a = K\left(\sqrt[2^V]{a}\right)
$$

を考えると楽になります。
このとき、素イデアルにおける Frobenius は高次剰余記号を与えるものとして理解でき、$a$ がどの段階まで $2^j$ 乗剰余になるかが、拡大 $L_a/K$ における分解挙動として符号化されます。
したがって、集合$A=${$2,4,\dots,98$}に対して合成体
$$
L=\prod_{a\in A} L_a
$$
を考え、$p$ と同じ Frobenius 条件を満たすように $q$ を選べれば、各 $a\in A$ について
$$
v_2(ord_q(a))=v_2(ord_p(a))
$$
を同時にみたすように設計できるはずです。

これらの考えをもとに実装したものが以下になります。
ただし、実装の関係で上記の説明とずれている個所があります。

• $V \geq 6$のとき$V=5$として計算
• もし見つからなったら総当たり
  また、実装力がなく成功率は低いです。
  作問者Writeupが上がったら参考にしながら別解という形で書き加えます。

from __future__ import annotations

import random
import re
import sys

from dataclasses import dataclass
from sage.all import *
from pwn import *


ODD_PRIMES = (3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47)
P_ODD = 1
for _pr in ODD_PRIMES:
    P_ODD *= _pr


EVEN_BASES = tuple(range(2, 100, 2))

U_STAGE0 = (2,)
U_STAGE1 = (2, 6, 10, 14, 22, 26, 34, 38, 46, 58, 62, 74, 82, 86, 94)
U_STAGE2 = tuple(sorted(set(U_STAGE1 + (4, 8, 12, 16, 18, 20, 24, 28, 30, 32, 36, 40, 42, 44, 48, 50, 54, 56, 60, 64, 66, 68, 70, 72, 76, 78, 80, 84, 88, 90, 92, 96, 98))))

def v2(x: int) -> int:
    c = 0
    while x % 2 == 0:
        x //= 2
        c += 1
    return c


def u_val_mod_prime(p: int, a: int) -> int:
    V = v2(p - 1)
    P = Integer(p)
    exp = (P - 1) >> V
    x = Integer(a % p).powermod(exp, P)
    if x == 1:
        return 0
    u = 0
    while x != 1:
        x = (x * x) % P
        u += 1
        if u > V:
            raise RuntimeError("unexpected: exceeded V in 2-Sylow walk")
    return u


def u_vec(p: int) -> tuple[int, ...]:
    return tuple(u_val_mod_prime(p, a) for a in EVEN_BASES)


def _u_vec_target_map(p: int) -> dict[int, int]:
    uv = u_vec(p)
    return {a: uv[i] for i, a in enumerate(EVEN_BASES)}


def _u_matches_stage(q: int, V: int, target: dict[int, int], bases: tuple[int, ...]) -> bool:
    for a in bases:
        if u_val_mod_prime(q, a) != target[a]:
            return False
    return True


_SMALL_PRIMES = tuple(int(p) for p in primes_first_n(200))


def _passes_small_prime_sieve(n: int) -> bool:
    for p in _SMALL_PRIMES:
        if n % p == 0:
            return n == p
    return True

def legendre_symbol(a: int, p: int) -> int:
    r = pow(a % p, (p - 1) // 2, p)
    return -1 if r == p - 1 else 1


def crt_pair(a1: int, m1: int, a2: int, m2: int) -> tuple[int, int]:
    g = int(gcd(m1, m2))
    if g != 1:
        raise ValueError("moduli must be coprime")
    t = ((a2 - a1) * pow(m1, -1, m2)) % m2
    x = a1 + m1 * t
    return x % (m1 * m2), m1 * m2


def crt_many(congruences: list[tuple[int, int]]) -> tuple[int, int]:
    x, m = congruences[0]
    x %= m
    for a, mod in congruences[1:]:
        x, m = crt_pair(x, m, a % mod, mod)
    return x, m


def find_q_v1(p: int, bits: int) -> int:
    if v2(p - 1) != 1:
        raise ValueError("p must satisfy v2(p-1)=1")

    lp2 = legendre_symbol(2, p)
    q_mod8 = 7 if lp2 == 1 else 3

    congruences: list[tuple[int, int]] = [(q_mod8, 8)]
    for l in ODD_PRIMES:
        target = legendre_symbol(l, p)
        want_q_over_l = target if (l % 4) == 1 else -target
        chosen = None
        for r in range(1, l):
            if legendre_symbol(r, l) == want_q_over_l:
                chosen = r
                break
        if chosen is None:
            raise RuntimeError("no residue found (unexpected)")
        congruences.append((chosen, l))

    x, mod = crt_many(congruences)
    lo = 1 << (bits - 1)
    hi = 1 << bits

    k = max(1, (lo - x + mod - 1) // mod)
    q = x + k * mod
    while q < hi:
        if q != p and q.bit_length() == bits and _passes_small_prime_sieve(q) and Integer(q).is_prime(proof=False):
            return q
        q += mod
    raise RuntimeError("failed to find 1024-bit q for V=1 (increase search range)")


@dataclass(frozen=True)
class CycloParams:
    coeff_bits: int
    tries: int
    modulus_int: int


def _default_cyclo_params(bits: int, deg: int, V: int) -> CycloParams:
    base = 1 << (V + 4)
    target_per_deg = max(1, bits // max(1, deg))

    if V == 5:
        prod = 1
        for pr in ODD_PRIMES:
            if pr == 47:
                break
            prod *= pr
        modulus_int = base * prod
        coeff_bits = 2
        tries = 250000
        return CycloParams(coeff_bits=coeff_bits, tries=tries, modulus_int=modulus_int)

    prod = 1
    for pr in ODD_PRIMES:
        if (base * prod * pr).bit_length() - 1 <= target_per_deg:
            prod *= pr
        else:
            break
    modulus_int = base * prod

    log2_mod = modulus_int.bit_length() - 1
    coeff_bits = max(2, target_per_deg - log2_mod + 3)

    tries = 60000 if deg <= 8 else 120000
    if deg >= 16:
        tries = 200000

    return CycloParams(coeff_bits=coeff_bits, tries=tries, modulus_int=modulus_int)


def _find_principal_prime_generator(K: NumberField, p: int, m: int) -> NumberFieldElement:
    z = K.gen()
    exp = (p - 1) // m
    rng = random.Random(p ^ (m << 32))

    for _ in range(512):
        a = rng.randrange(2, p - 1)
        r = pow(a, exp, p)
        if r == 1:
            continue
        if pow(r, m >> 1, p) == 1:
            continue
        Pideal = K.ideal(p, z - int(r))
        try:
            if not bool(Pideal.is_principal()):
                continue
            g = Pideal.gens_reduced()[0]
            if int(abs(ZZ(g.norm()))) == p:
                return g
        except Exception:
            continue

    fac = K.ideal(p).factor()
    for (Pideal, _) in fac:
        try:
            if bool(Pideal.is_principal()):
                g = Pideal.gens_reduced()[0]
                if int(abs(ZZ(g.norm()))) == p:
                    return g
        except Exception:
            continue
    raise RuntimeError("no principal prime ideal above p found")


def find_q_cyclotomic(p: int, bits: int) -> int:
    V_target = v2(p - 1)
    if V_target <= 1:
        raise ValueError("use find_q_v1 for V<=1")
    elif V_target <= 5:
        V_field = V_target
    else:
        V_field = 5

    m = 1 << V_field
    K = CyclotomicField(m)
    deg = K.degree()
    params = _default_cyclo_params(bits, deg, V_field)

    target_map = _u_vec_target_map(p)
    pi = _find_principal_prime_generator(K, p, m)

    z = K.gen()
    z_pows = [z**i for i in range(deg)]

    M = K.ideal(params.modulus_int)
    m_elem = M.gens_reduced()[0]

    rng = random.Random(p ^ (V_target << 16) ^ (V_field << 8))
    if V_field == 5:
        coeff_bit_candidates = (1, 2, 3)
    else:
        coeff_bit_candidates = tuple(max(2, params.coeff_bits + bump) for bump in range(-2, 3))

    log.info("[+] Computing q")
    for coeff_bits in coeff_bit_candidates:
        for _ in range(params.tries // 5):
            bound = 1 << (coeff_bits - 1)
            coeffs = [rng.getrandbits(coeff_bits) - bound for __ in range(deg)]
            t = sum(ZZ(coeffs[i]) * z_pows[i] for i in range(deg))

            cand = pi + m_elem * t
            q = int(abs(ZZ(cand.norm())))
            if q == p or q.bit_length() != bits:
                continue
            if not _passes_small_prime_sieve(q):
                continue
            if not Integer(q).is_prime(proof=False):
                continue
            if v2(q - 1) != V_target:
                continue
            if not _u_matches_stage(q, V_target, target_map, U_STAGE1):
                continue
            if not _u_matches_stage(q, V_target, target_map, U_STAGE2):
                continue
            if not _u_matches_stage(q, V_target, target_map, EVEN_BASES):
                continue
            return q

    raise RuntimeError("failed to find q (increase tries/adjust parameters)")


def find_q(p: int, bits: int = 1024) -> int:
    V = v2(p - 1)
    if V == 1:
        return find_q_v1(p, bits)
    if V >= 2:
        try:
            return find_q_cyclotomic(p, bits)
        except Exception:
            pass

    target_map = _u_vec_target_map(p)
    log.info("[+] Cant find q.")


    for i in range(200000):
        if i%10000 == 0:
            log.info(f"[*] i={i}")
        q = int(random_prime((1 << bits) - 1, lbound=1 << (bits - 1), proof=False))
        if q == p:
            continue
        if v2(q - 1) != V:
            continue
        if not _u_matches_stage(q, V, target_map, U_STAGE1):
            continue
        if not _u_matches_stage(q, V, target_map, EVEN_BASES):
            continue
        if u_vec(q) == tuple(target_map[a] for a in EVEN_BASES):
            return q
    raise RuntimeError("fallback random search failed")


P_RE = re.compile(r"^p\s*=\s*(\d+)\s*$")
HOST = "yukari-infinity.seccon.games"
PORT = 13910
io = remote(HOST, PORT)

i = 0
while True:
    try:
        line_b = io.recvline(timeout=context.timeout)
    except EOFError:
        break

    if not line_b:
        break

    try:
        sys.stdout.write(line_b.decode(errors="replace"))
    except Exception:
        sys.stdout.write(repr(line_b) + "\n")
    sys.stdout.flush()

    line = line_b.decode(errors="ignore").strip()
    m = P_RE.match(line)
    if not m:
        continue

    p = int(m.group(1))

    V = v2(p - 1)
    i = i+1
    log.info(f"[+] i={i}, V = {V}")

    q = find_q(p, bits=1024)

    if debug:
        log.info(f"[debug] computed q (V={V}) bits={q.bit_length()}")

    io.sendline(str(q).encode())

    if debug:
        log.info("[debug] sent q")

io.close()

[+] i=127, V = 1
[*] [debug] computed q (V=1) bits=1024
[*] [debug] sent q
q: error!
p = 91702515163758684036270912242390551837403332245304299652195656715046525761074972506862237728255372505453122187957954875517334311434135781983002797238929269726884347888816549245697078247719948694983312698528458364309768843113922545952418767822505952107428428729507737284761559914152260636827581934112495637739
[+] i=128, V = 1
[*] [debug] computed q (V=1) bits=1024
[*] [debug] sent q
q: error!
SECCON{67fd17bfdd32106ea5fff886bb26135dd697631f0aac89582a55b1330733fd6e}

Flag:SECCON{67fd17bfdd32106ea5fff886bb26135dd697631f0aac89582a55b1330733fd6e}

参考文献

• https://www.kurims.kyoto-u.ac.jp/~kenkyubu/bessatsu/open/B25/pdf/B25_008.pdf
• https://repository.kulib.kyoto-u.ac.jp/server/api/core/bitstreams/bae2694f-705c-447f-a387-6b81f24b57e5/content
• https://en.wikipedia.org/wiki/Chebotarev_density_theorem
• https://math.stackexchange.com/questions/103677/split-in-cyclotomic-field
• https://en.wikipedia.org/wiki/Field_norm
• https://pc1.math.gakushuin.ac.jp/~shin/html-files/Alg2/2022/13kummer.pdf
• 整数論２ 代数的整数論の基礎 (https://www.nippyo.co.jp/shop/book/6344.html)