pbctf 2021

October 10, 2021

The writeups are in progress, but here’s solution code!

GoodHash

The idea is to linearize field multiplication. Given two polynomials f and g, notice that their product expands into

f g = f_127 x_127 g + ... + f_0 x_0 g.

If we consider x_i g as a binary vector, it becomes clear that f g is a linear combination of 128 vectors, i.e. a matrix multiplication.

from Crypto.Cipher import AES
from Crypto.Cipher import _mode_gcm
from sage.crypto.util import ascii_to_bin, bin_to_ascii

K.<x> = GF(2^128, modulus=x^128 + x^7 + x^2 + x + 1)

def block2poly(block):
  assert len(block) == 16
  return K(list(str(ascii_to_bin(block.decode('raw_unicode_escape')))))

def poly2block(poly):
  bits = ZZ(poly.integer_representation()).bits()
  bits.extend([0]*(128 - len(bits)))
  return bin_to_ascii(bits).encode('raw_unicode_escape')

subkey = AES.new(b'goodhashGOODHASH', AES.MODE_ECB).encrypt(bytes(16))
H = block2poly(subkey)

def ghash(blocks, debug=False):
  if debug:
    return poly2block(sum([block2poly(block) * H^(len(blocks) - i) for i, block in enumerate(blocks)]))
  else:
    ghash = _mode_gcm._GHASH(subkey, _mode_gcm._ghash_portable)
    ghash.update(b''.join(blocks))
    return ghash.digest()

token = b'{"token": "db0042de22975348b132516275ab489a", "admin": false}'
token += bytes(3) + bytes(8) + int(61*8).to_bytes(8, 'big')
blocks = [token[i:i+16] for i in range(0, len(token), 16)]
polys = list(map(block2poly, blocks))

diffs = [0 for _ in polys]
diffs[4] = polys[4] + block2poly(bytes(8) + int(63*8).to_bytes(8, 'big'))
diffs[3] = polys[3] + block2poly(b'","admin":true}\x00')

monomials = [x^i for i in range(128) if i % 8 > 3 or i == 19]
M = column_matrix([vector(m*H^3) for m in monomials] + [vector(m*H^4) for m in monomials])
print(f'rank(M) = {M.rank()}')
v = vector(diffs[3]*H^2 + diffs[4]*H)
u = M.solve_right(v)

diffs[2] = sum(map(prod, zip(u[:len(u)//2], monomials)))
diffs[1] = sum(map(prod, zip(u[len(u)//2:], monomials)))

collision = list(map(poly2block, map(sum, zip(polys, diffs))))
print(b''.join(collision[:-1]).rstrip(b'\x00'))
assert ghash(blocks) == ghash(collision)

Seed Me

Notice that the n-th output of an LCG remains an affine transformation of the original seed, just with different coefficients:

a(a( ... (a(ax + c) + c) ... ) + c) + c ==
  a^n x + a^(n-1) c + a^(n-2) c + ... + ac + c

So our lucky seeds are really just the result of a matrix multiplication, which must be bounded in each entry.

[ : ]    [ 1  0 0 ... 0] [ x  ]   [ 0   ]   [ : ]
[ : ]    [ a1 q 0 ... 0] [ k1 ]   [ c1  ]   [ : ]
[ l ] <= [ a2 0 q ... 0] [ k2 ] + [ c2  ] < [ u ]
[ : ]    [ :  : :  :  :] [ :  ]   [ :   ]   [ : ]
[ : ]    [ an 0 0 ... q] [ kn ]   [ c15 ]   [ : ]

I’ve phrased this to translate nicely into rkm’s inequality solver, but the idea is to solve CVP where the target is the vector bounds’ component-wise average. This in turn reduces to performing lattice reduction on the matrix.

Plugging directly into the solver fails, but it’s almost correct! The algorithm finds a vector that succeeds for all but the first and last seeds, which fall slightly (like 0.05%) out of bounds. My ad hoc solution was to force LLL to prioritize individual rows by adjusting the weights. The solver already this based on the distance between each inequality’s bounds — in this case, the weights are equal. I manually modified them until the solver gave a valid solution.

if i == 0:
  ineq_weight = 4
elif i == 15:
  ineq_weight = 3
else:
  ineq_weight = 2
load('inequality.sage')

multiplier = 0x5DEECE66DL
addend = 0xBL
q = 2^48
R = Integers(q)

P.<x> = R[]
coeffs = [[0, 1]]
seed = x
for _ in range(15):
  for _ in range(0o3777):
    seed = multiplier*seed + addend
  seed = multiplier*seed + addend
  coeffs.append(seed.coefficients())
c, a = zip(*coeffs)

def derive_seed(x):
  x = R(x)
  ai = R(multiplier)^-1
  x = ai*(x - addend)
  for _ in range(0o3777):
    x = ai*(x - addend)
  return x.lift() ^^ multiplier

c = vector(ZZ, c)
a = vector(ZZ, a)

M = diagonal_matrix(ZZ, [q]*len(c), sparse=False).delete_rows([0]).stack(a)
lb = vector(ZZ, [16444267 << 24]*len(c)) - c
ub = vector(ZZ, [q - 1]*len(c)) - c

out, sol = solve(M, lb, ub)
x = sol[-1]
print(out)
print(sol)
print(derive_seed(x))

Yet Another PRNG

import itertools

M = 2^64 - 59
m1 = 2^32 - 107
m2 = 2^32 - 5
m3 = 2^32 - 209
a1 = [4256, 307568, 162667]
a2 = [593111, 526598, 630723]
a3 = [383732, 73391, 955684]

P = PolynomialRing(ZZ, 'x0,x1,x2,y0,y1,y2,z0,z1,z2,'+','.join([f'k{i}' for i in range(9*3 + 12*2)]))
gens = list(P.gens())

x = gens[0:3]
y = gens[3:6]
z = gens[6:9]
k = gens[9:]

for _ in range(9):
  x.append(sum(map(prod, zip(x[-3:], a1))) - k.pop(0)*m1)
  y.append(sum(map(prod, zip(y[-3:], a2))) - k.pop(0)*m2)
  z.append(sum(map(prod, zip(z[-3:], a3))) - k.pop(0)*m3)

def solve(o):
  diffs = []
  for i in range(len(o)):
    diffs.append(x[i] - z[i] - o[i]*204.inverse_mod(m3) % m3 + k[2*i]*m3)
    diffs.append(y[i] - z[i] - o[i]*102.inverse_mod(m1) % m1 + k[2*i+1]*m1)

  M, _ = Sequence(diffs).coefficient_matrix()
  B = matrix(M.right_kernel().basis())

  factors = [2^32]*9 + [2^20]*(len(o) - 3)*3 + [1]*len(o)*2 + [1]

  B = B.change_ring(QQ)
  for i, factor in enumerate(factors):
    B.rescale_col(i, 1/factor)

  B = B.LLL()

  for i, factor in enumerate(factors):
    B.rescale_col(i, factor)
  B = B.change_ring(ZZ)

  for v in B:
    v *= v[-1]
    if v[-1] == 1 and all(map(lambda v: v in (0, 1), v[-len(o)*2-1:-1])):
      return v[0:3] % m1, v[3:6] % m2, v[6:9] % m3

def brute(out):
  valid = set()
  for k1 in (0, 1, -1, -2):
    o = out - k1*M
    for k2 in itertools.product((0, -1), repeat=3):
      yz = o*102.inverse_mod(m1) % m1 + k2[0]*m1
      yx = o*204.inverse_mod(m2) % m2 + k2[2]*m2
      xz = o*204.inverse_mod(m3) % m3 + k2[1]*m3
      if yz - yx == xz:
        valid.add(o)
  return valid

hint = bytes.fromhex('67f19d3da8af1480f39ac04f7e9134b2dc4ad094475b696224389c9ef29b8a2aff8933bd3fefa6e0d03827ab2816ba0fd9c0e2d73e01aa6f184acd9c58122616f9621fb8313a62efb27fb3d3aa385b89435630d0704f0dceec00fef703d54fca')
valid = [brute(int.from_bytes(hint[i:i+8], byteorder='big')) for i in range(0, 48, 8)]
for guess in itertools.product(*valid):
  print(solve(guess))

Yet Another RSA

load('coppersmith.sage')

N = 144...997
e = 370...289

R = Integers(e)
P.<k, s> = PolynomialRing(R)

bounds = (2^400, 2^512)
f = k*(N^2 - N + N*2*s + (2*s)^2 + 2*s + 1) + 1
(k, s), = small_roots(f, bounds, m=3, d=4)
print(f.change_ring(ZZ)(k.lift(), s.lift()) // e)