Birds of Randomness — HackTheBox

Description

"You think you found your soulmate, but you lost the train ticket that your love gave you to meet again. The only thing that you remember is the coordinates of the departing station and your love disappearing into the steam of the train."

The server generates a random 256-bit "ticket number" as a seed for a Wichmann-Hill PRNG (three coupled LCGs with moduli 30269, 30307, 30323). It rotates the PRNG until all three components (x, y, z) are prime, computes scalar d = x·y·z, and returns the EC point d·G as the "departing station coordinates." It then does the same for the "destination" — which we must guess in at most 6 attempts.

Files provided: source.py, love_letter.txt (flavor text warning not to do stego on forologia.jpg).

Analysis

Elliptic Curve Parameters

E: y² = x³ + 2x + 3  over  F_p
p = 17101937747109687265202713197737423  (~113 bits)
G = (3543321030468950376213178213609418, 14807290861072031659976937040569354)
n = 17101937747109687496599931614463506  (curve order)

Critical Vulnerability: Smooth Curve Order

The curve order factors completely into small primes:

n = 2 × 3 × 7² × 1487 × 3761 × 176489 × 439693 × 3113111 × 43054831

The largest prime factor is only ~4.3×10⁷. This makes the ECDLP trivially solvable via the Pohlig-Hellman algorithm with Baby-step Giant-step in each subgroup. Total work is roughly √(largest_factor) ≈ 6500 EC operations — instant.

PRNG Structure (Wichmann-Hill 1982)

x = (171 * x) % 30269
y = (172 * y) % 30307
z = (170 * z) % 30323

The PRNG rotates until all three components are simultaneously prime, then outputs d = x·y·z as the EC scalar. Since each component is bounded by ~30000, the scalar d < 30268 × 30306 × 30322 ≈ 2.78×10¹³, much smaller than the curve order ~1.7×10³⁴.

Attack Plan

1. Pohlig-Hellman on the departing point P = d₁·G to recover d₁
2. Factor d₁ into a prime triple (x₁, y₁, z₁) with each ≤ 30322
3. Try all 6 permutations of the triple as PRNG state, continue PRNG to predict destination
4. 6 guesses = 3! permutations — guaranteed to find the correct one

Solution

Step 1: Pohlig-Hellman ECDLP

For each prime power pᵢᵉⁱ dividing n:

1. Project both G and P into the subgroup of order pᵢᵉⁱ by multiplying by n/pᵢᵉⁱ
2. Solve the DLP in the small subgroup via BSGS (for prime powers, use p-adic lifting)
3. Combine all residues via CRT to get d₁ mod n

Since d₁ < 2.78×10¹³ ≪ n, the result is the actual scalar (or its negation n - d₁).

Step 2: Factor d₁ into Prime Triple

d₁ = x₁ · y₁ · z₁ where all are primes in [2, 30322]. Iterate over primes p₁ ≤ 30268, check divisibility, then check if d₁/p₁ factors into two primes both ≤ 30322. This yields exactly one unordered triple with 6 ordered permutations.

Step 3: Predict Destination

For each of the 6 permutations (x₁, y₁, z₁):

1. Set PRNG state to (x₁, y₁, z₁)
2. Call goToNextStation() — rotate until all-prime, compute d₂·G
3. Submit the resulting point as a guess

One of the 6 will match the server's actual PRNG state.

#!/usr/bin/env python3
"""
Birds of Randomness - HackTheBox Crypto Challenge
Attack: Pohlig-Hellman ECDLP (smooth order) → recover PRNG state → predict next output.
"""

import socket
import re
import math
import time
from ecdsa import ellipticcurve as ecc
from Crypto.Util.number import isPrime

# Curve parameters (from source)
p_curve = 17101937747109687265202713197737423
a_coeff = 2
b_coeff = 3
ec_order = 17101937747109687496599931614463506
Gx = 3543321030468950376213178213609418
Gy = 14807290861072031659976937040569354

E = ecc.CurveFp(p_curve, a_coeff, b_coeff)
G = ecc.Point(E, Gx, Gy, ec_order)
INFINITY = ecc.INFINITY

# Order factorization: 2 * 3 * 7^2 * 1487 * 3761 * 176489 * 439693 * 3113111 * 43054831
ORDER_FACTORS = [
    (2, 1), (3, 1), (7, 2), (1487, 1), (3761, 1),
    (176489, 1), (439693, 1), (3113111, 1), (43054831, 1),
]

MAX_D = 30268 * 30306 * 30322  # ~2.78e13


class TrainRoute:
    def __init__(self, seed_tuple):
        self.seed = seed_tuple

    def rotate(self):
        x, y, z = self.seed
        x = (171 * x) % 30269
        y = (172 * y) % 30307
        z = (170 * z) % 30323
        self.seed = x, y, z

    def goToNextStation(self):
        while True:
            self.rotate()
            x, y, z = self.seed
            if isPrime(x) and isPrime(y) and isPrime(z):
                d = x * y * z
                new_point = d * G
                return int(new_point.x()), int(new_point.y())


def point_key(P):
    if P == INFINITY:
        return "INF"
    return (int(P.x()), int(P.y()))


def bsgs_subgroup(base, target, n):
    """Baby-step Giant-step: find d in [0, n-1] such that target = d * base."""
    if n <= 1:
        return 0
    m = math.isqrt(n) + 1
    baby = {}
    jP = INFINITY
    for j in range(m):
        baby[point_key(jP)] = j
        jP = jP + base
    mP = m * base
    neg_mP = (
        ecc.Point(E, int(mP.x()), (-int(mP.y())) % p_curve, ec_order)
        if mP != INFINITY else INFINITY
    )
    gamma = target
    for i in range(m + 1):
        k = point_key(gamma)
        if k in baby:
            return (i * m + baby[k]) % n
        gamma = gamma + neg_mP
    return None


def pohlig_hellman(base, target, order, factors):
    """Solve target = d * base using Pohlig-Hellman + CRT."""
    residues, moduli = [], []
    for pi, ei in factors:
        pe = pi**ei
        cofactor = order // pe
        base_sub = cofactor * base
        target_sub = cofactor * target
        if ei == 1:
            di = bsgs_subgroup(base_sub, target_sub, pi)
        else:
            di = 0
            gamma = target_sub
            base_p = (order // pi) * base
            for k in range(ei):
                exp = pe // (pi ** (k + 1))
                if exp == 0:
                    exp = 1
                h = exp * gamma
                dk = bsgs_subgroup(base_p, h, pi)
                di += dk * (pi**k)
                gamma = gamma + ((-dk * (pi**k)) * base_sub)
        residues.append(di)
        moduli.append(pe)

    # CRT
    M = 1
    for m in moduli:
        M *= m
    result = 0
    for ri, mi in zip(residues, moduli):
        Mi = M // mi
        yi = pow(Mi, -1, mi)
        result += ri * Mi * yi
    return result % M


def factor_d_into_primes(d):
    """Find all (x, y, z) with d = x*y*z, all prime, in PRNG ranges."""
    results = []
    primes = [n for n in range(2, 30323) if isPrime(n)]
    primes_set = set(primes)
    for px in primes:
        if px > 30268:
            break
        if d % px != 0:
            continue
        remainder = d // px
        for py in primes:
            if py > 30306:
                break
            if remainder % py != 0:
                continue
            pz = remainder // py
            if 2 <= pz <= 30322 and pz in primes_set:
                results.append((px, py, pz))
    return results


def recv_until(sock, marker):
    data = b""
    while marker not in data:
        chunk = sock.recv(4096)
        if not chunk:
            break
        data += chunk
    return data.decode()


def solve():
    host = "154.57.164.66"
    port = 30535

    sock = socket.socket(socket.AF_INET, socket.SOCK_STREAM)
    sock.settimeout(30)
    sock.connect((host, port))

    data = recv_until(sock, b"x coordinate:")
    coord_match = re.search(r"\((\d+),\s*(\d+)\)", data)
    Px, Py = int(coord_match.group(1)), int(coord_match.group(2))
    P_point = ecc.Point(E, Px, Py, ec_order)

    # Step 1: Pohlig-Hellman ECDLP
    d1 = pohlig_hellman(G, P_point, ec_order, ORDER_FACTORS)
    if d1 > MAX_D:
        d1 = ec_order - d1  # try negation

    # Step 2: Factor into prime triple
    factorizations = factor_d_into_primes(d1)

    # Step 3: Predict destination for each permutation
    candidates = []
    seen = set()
    for x1, y1, z1 in factorizations:
        tr = TrainRoute((x1, y1, z1))
        dest = tr.goToNextStation()
        if dest not in seen:
            seen.add(dest)
            candidates.append(dest)

    # Step 4: Submit guesses (up to 6)
    for dx, dy in candidates[:6]:
        sock.sendall(f"{dx}\n".encode())
        resp = recv_until(sock, b"y coordinate:")
        sock.sendall(f"{dy}\n".encode())
        resp = recv_until(sock, b"x coordinate:") if True else sock.recv(4096).decode()
        if "HTB{" in resp:
            flag = re.search(r"HTB\{[^}]+\}", resp)
            print(f"FLAG: {flag.group()}")
            break

    sock.close()


if __name__ == "__main__":
    solve()