Tags: hyperelliptic-curves crypto sagemath
Rating:
## The challenge
The challenge constructs a random number generator using a Hyperelliptic curve of genus 3 (and arithmetic on it), then generates some bytes using it which it XORs with the message `Hello! The flag is: hxp{...}`.
The code is in [Sage](https://sagemath.org) and first creates a hyperelliptic curve $C$ of genus $3$ of the form
$y^2 = x^7 + x$ over a 64-bit prime field. What is interesting about hyperelliptic curves is that their set of $\mathbb{K}$-rational points (for $\mathbb{K}$ being some extension of the base field) does not form a group like it is for elliptic curves. There is however a group associated to each hyperelliptic curve, its Jacobian. The elements of the Jacobian are divisors, one might imagine them as
a formal sum $ \sum_i^m n_i P_i $ of some points $ P_i $ with $ n_i \in \mathbb{Z} $. To my best understanding when mathematicians are talking about a formal sum they are talking about doing a sum, but then instead of doing it, they ... *don't do it*. In this case it also cannot be computed, as there is no addition defined on the points of the hyperelliptic curve. The sum in the divisor can contain a number of points upto and equal to the genus of the curve, some might however be defined over extensions of the base field. The operation defined on the divisors is too complicated and unnecessary for this write-up, but more can be found in the [SageMath docs](https://doc.sagemath.org/html/en/reference/arithmetic_curves/sage/schemes/hyperelliptic_curves/jacobian_morphism.html) and in the helpful [An Introduction to Elliptic and Hyperelliptic Curve Cryptography and the NTRU Cryptosystem](https://www.esat.kuleuven.be/cosic/publications/article-520.pdf).
```python
#!/usr/bin/env sage
import struct
from random import SystemRandom
p = 10000000000000001119
R.<x> = GF(p)[]; y=x
f = y + prod(map(eval, 'yyyyyyy'))
C = HyperellipticCurve(f, 0)
J = C.jacobian()
```
The RNG is seeded with three random integers in the range $[0, p^3]$ which makes them roughly 192 bits. During initialization, the RNG also creates three fixed points on the curve (with $x$-coords $11, 22$ and $33$) and transforms them into a divisors on the Jacobian. The RNG calls the `clk` function and processes the results whenever it runs out of bytes. The generation of the first byte triggers this clocking. The clocking takes the current list of the three divisors and multiplies them with the three secret random integers in `self.es`. These three divisors are then summed and this sum is decomposed into two elements (`u, v = sum(self.clk())`).
In Sage, divisors on the Jacobian are represented using the Mumford representation. This representation consists of two polynomials $u(x)$ and \$v(x)\$. It holds that \$u(x)\$ is monic, \$u(x)\$ divides \$f(x) - h(x) v(x) - v(x)^2\$ where \$f(x)\$ and \$h(x)\$ are polynomials associated with the curve. It also holds that \$deg(v(x)) < deg(u(x)) \le g\$ where \$g\$ is the genus of the curve, in our case \$3\$. When one has a divisor on the Jacobian in Sage and indexes into it or decomposes it into two elements the two polynomials \$u(x)\$ and \$v(x)\$ are returned, which is what happens on the `u, v = sum(self.clk())` line.
The rest of the RNG is straightforward, it takes the six coefficients from the Mumford representation of the summed divisor (three from \$u(x)\$ and three from \$v(x)\$) and converts them to 8-byte integers, then concatenates them into a 48-byte string which used as the RNG output.
```python
class RNG(object):
def __init__(self):
self.es = [SystemRandom().randrange(p**3) for _ in range(3)]
self.Ds = [J(C(x, min(f(x).sqrt(0,1)))) for x in (11,22,33)]
self.q = []
def clk(self):
self.Ds = [e*D for e,D in zip(self.es, self.Ds)]
return self.Ds
def __call__(self):
if not self.q:
u,v = sum(self.clk())
rs = [u[i] for i in range(3)] + [v[i] for i in range(3)]
assert 0 not in rs and 1 not in rs
self.q = struct.pack('<'+'Q'*len(rs), *rs)
r, self.q = self.q[0], self.q[1:]
return r
def __iter__(self): return self
def __next__(self): return self()
```
The rest of the challenge uses the RNG to compute 48 random bytes and XORs the output with the message.
```python
flag = open('flag.txt').read().strip()
import re; assert re.match(r'hxp\{\w+\}', flag, re.ASCII)
text = f"Hello! The flag is: {flag}"
print(bytes(k^^m for k,m in zip(RNG(), text.encode())).hex())
```
## The solution
The first important realization is that due to the message starting with `Hello! The flag is: hxp{` the first 24 bytes of the output of the RNG are known. In the algorithm, this corresponds to the knowledge of the `u[i] for i in range(3)` in the `__call__` function. These `u[i]` represent the coefficients of the \$u(x)\$ polynomial in the Mumford representation of the divisor. What we do not know are the next 24 bytes of the RNG output, which correspond to the coefficients of the \$v(x)\$ polynomial in the Mumford representation of the divisor, so our task is to somehow recover them.
To recover the \$v(x)\$ polynomial, one needs to think about the correspondence between the points in the divisor and its Mumford representation. It holds that roots of the \$u(x)\$ polynomial are \$x\$-coordinates of the points in the formal sum in the divisor. Furthermore, for \$x_i\$ a root of \$u(x)\$, \$y_i = v(x_i)\$ is the \$y\$-coordinate of the represented point. We can use this recover the \$v(x)\$ polynomial. The \$v(x)\$ polynomial has degree smaller than three so we can represent it as \$a x^2 + b x + c\$ with \$a, b, c \in \mathbb{K}\$. Using the three roots \$x_i\$ of the \$u(x)\$ polynomial we can form three linear equations over \$\mathbb{\overline{K}}\$ of the form \$v^2(x) = x^7 + x\$. Plugging in our form of the polynomial \$v(x)\$ we get \$(a x^2 + b x + c)^2 = x^7 + x\$, for each \$x\$ from the roots of the \$u(x)\$ polynomial. For easier solving we will take the square root of the equation and get \$a x^2 + b x + c = (x^7 + x)^{1/2}\$. This is what the Sage script below does, with some additional fiddling of the ordering of the coefficients of \$v(x)\$ when they are transformed into bytes for the RNG output.
```python
#!/usr/bin/env sage
from binascii import unhexlify, hexlify
from itertools import permutations
p = 10000000000000001119
k = GF(p)
kc = k.algebraic_closure()
R.<x> = k[]; y=x
f = y + prod(map(eval, 'yyyyyyy'))
C = HyperellipticCurve(f, 0)
J = C.jacobian()
msg_prefix = b"Hello! The flag is: hxp{"
with open("output.txt") as f:
content = unhexlify(f.read().strip())
bs = []
for c, m in zip(content, msg_prefix):
# Do the XOR, obtain k
b = c^^m
print(b)
bs.append(b)
u0 = int.from_bytes(bytes(bs[:8]), byteorder="little")
u1 = int.from_bytes(bytes(bs[8:16]), byteorder="little")
u2 = int.from_bytes(bytes(bs[16:24]), byteorder="little")
print(hex(u0), hex(u1), hex(u2))
ps = x^3 + u2 * x^2 + u1 * x + u0 # TODO: this ordering might be the other way around.
aps_roots = ps.roots(ring=kc, multiplicities=False)
x0, x1, x2 = aps_roots
A = Matrix(((x0^2, x0, kc(1)), (x1^2, x1, kc(1)), (x2^2, x2, kc(1))))
Y = vector((x0^7 + x0, x1^7 + x1, x2^7 + x2))
Ys = vector((-Y[0].sqrt(), -Y[1].sqrt(), -Y[2].sqrt())) # TODO: Maybe the other sqrt?
v = A.solve_right(Ys)
print(v)
v0 = int(str(v[0])).to_bytes(8, byteorder="little")
v1 = int(str(v[1])).to_bytes(8, byteorder="little")
v2 = int(str(v[2])).to_bytes(8, byteorder="little")
for a0, a1, a2 in permutations((v0, v1, v2)):
cs = []
q = bytes(bs) + a0 + a1 + a2
for c, m in zip(content, q):
# Do the XOR, obtain k
b = c^^m
cs.append(chr(b))
print("".join(cs))
```
And there it is: `Hello! The flag is: hxp{ez_P4rT_i5_ez__tL0Cm}`.
