# AeroCTF 2021 - Horcrux¶

I didn't participate in the AeroCTF but looked at the challenges afterwards. Please see nice and creative solutions by rkm0959 and Mystiz. Here I will describe an alternativ solution to the Horcrux challenge. The first part of my solution is simpler by using resultants.

However, I realized that I overcomplicated the second part a lot. I keep this writeup to keep the techniques which might be useful in other scenarios.

Here is the code:

#!/usr/bin/env python3.8

from os import urandom
from gmpy2 import next_prime
from random import randrange, getrandbits
from Crypto.Cipher import AES
from fastecdsa.curve import Curve

def bytes_to_long(data):
return int.from_bytes(data, 'big')

def generate_random_point(p):
while True:
a, x, y = (randrange(0, p) for _ in range(3))
b = (pow(y, 2, p) - pow(x, 3, p) - a * x) % p

if (4 * pow(a, 3, p) + 27 * pow(b, 2, p)) % p != 0:
break

return Curve(None, p, a, b, None, x, y).G

class DarkWizard:
def __init__(self, age):
self.power = int(next_prime(getrandbits(age)))
self.magic = generate_random_point(self.power)
self.soul = randrange(0, self.power)

def create_horcrux(self, location, weakness):
# committing a murder

# splitting the soul in half
self.soul = self.soul * pow(2, -1, self.power) % self.power

# making a horcrux
horcrux = (self.soul + murder) * self.magic

# nobody should know location and weakness of the horcrux
horcrux.x ^= location
horcrux.y ^= weakness

return horcrux

def cast_spell(self, spell_name):
spell = bytes_to_long(spell_name)

return spell %~ spell

def encrypt(key, plaintext):
cipher = AES.new(key=key, mode=AES.MODE_ECB)
padding = b'\x00' * (AES.block_size - len(plaintext) % AES.block_size)

def main():
wizard_age = 3000
horcruxes_count = 2

wizard = DarkWizard(wizard_age)
print(f'Wizard\'s power:\n{hex(wizard.power)}\n')
print(f'Wizard\'s magic:\n{wizard.magic}\n')

key = urandom(AES.key_size[0])
horcrux_length = len(key) // horcruxes_count

for i in range(horcruxes_count):
key_part = key[i * horcrux_length:(i + 1) * horcrux_length]

horcrux_location = bytes_to_long(key_part[:horcrux_length // 2])
horcrux_weakness = bytes_to_long(key_part[horcrux_length // 2:])

horcrux = wizard.create_horcrux(horcrux_location, horcrux_weakness)
print(f'Horcrux #{i + 1}:\n{horcrux}\n')

with open('flag.txt', 'rb') as file:

ciphertext = encrypt(key, flag)
print(f'Ciphertext:\n{ciphertext.hex()}')

if __name__ == '__main__':
main()


Data:

Wizard's power:

Wizard's magic:
(On curve <None>)

Horcrux #1:
(On curve <None>)

Horcrux #2:
Y: 0x29cfc760529c48d68bc086a5b4403f1d4446db04abe243c99baf659b5e67cd6cdac9a658f273f4c682b9e13dda72aaa1ede42c69c98640f2fc58eabeef143a65334e4a236a6e72a157d92ab1a541c9bcf7d953b386a40d68880312ed1e900f6ba481bf134515f5a4c245a0d924db0e5105fc44fae71fc991df85219403ba5d48f68d5d91151f8411b4cbb6971bd63030b523989fa7ec94c14136ac712d4e91eca4c930d79caab7c328043c762917c4b3f868ce037cae9f19a579272c6b13ca8c19d349f5777bf6ac9c078d8472c92582daf96b30d4f7b8fdc004a36b792c133e2b6511956480892636ea91a3361afd8a3afbe3a5bf889feb5a5dc143f5a917347591634218066f2f71b36afd5257f6637152ac9a0965daa881ddc86ca8a8545b255cbb14e7738297a55c7428d9a0e79dee37f801fc1d49d205be809e0cd1b8a1b9d1d5b1b1a9ae98e7c564718cb17d10a3dd01c2914d7bb96c45a4fc942c2ed628354882464407c3208938282471aa2b8016e46c998e4
(On curve <None>)

Ciphertext:
b6094505efd8e2824e8cc8698e5e68e3b2c306a2c8179fbefbfb7d2fc1a0cfb921a89f94dde88b04e655ad87d15efcc7466af652d6a330e92babceca94be63b126702153ad83d24baf7d8848bb71202771a2d80870fbd17462a715906dfc63bc

What we have reduces to equations of 3 points on an unknown elliptic curve: $$y_0^2 = x_0^3 + ax_0+ b,$$ $$(y_1\oplus\beta_1)^2 = (x_1\oplus\alpha_1)^3 + a(x_1\oplus\alpha_1)+ b,$$ $$(y_2\oplus\beta_2)^2 = (x_2\oplus\alpha_2)^3 + a(x_2\oplus\alpha_2)+ b,$$ where $a,b \in \mathbb{F}_p$, $p$ is a 3000-bit known prime, $\alpha_1,\alpha_2,\beta_1,\beta_2$ are 32-bit random integers. The idea here is that XORing is equivalent to adding or subtracting a 32-bit value in the field. So we will simply replace the XOR with addition in the field (recovering the original variables is trivial).

Let's setup the data:

In [1]:
from sage.all import *
proof.arithmetic(False)

F = GF(p)

y2 = F(0x29cfc760529c48d68bc086a5b4403f1d4446db04abe243c99baf659b5e67cd6cdac9a658f273f4c682b9e13dda72aaa1ede42c69c98640f2fc58eabeef143a65334e4a236a6e72a157d92ab1a541c9bcf7d953b386a40d68880312ed1e900f6ba481bf134515f5a4c245a0d924db0e5105fc44fae71fc991df85219403ba5d48f68d5d91151f8411b4cbb6971bd63030b523989fa7ec94c14136ac712d4e91eca4c930d79caab7c328043c762917c4b3f868ce037cae9f19a579272c6b13ca8c19d349f5777bf6ac9c078d8472c92582daf96b30d4f7b8fdc004a36b792c133e2b6511956480892636ea91a3361afd8a3afbe3a5bf889feb5a5dc143f5a917347591634218066f2f71b36afd5257f6637152ac9a0965daa881ddc86ca8a8545b255cbb14e7738297a55c7428d9a0e79dee37f801fc1d49d205be809e0cd1b8a1b9d1d5b1b1a9ae98e7c564718cb17d10a3dd01c2914d7bb96c45a4fc942c2ed628354882464407c3208938282471aa2b8016e46c998e4)

BITS = 32

In [2]:
# SOLVING
R = PolynomialRing(F, names='aa1, aa2, bb1, bb2, a, b')
aa1, aa2, bb1, bb2, a, b = R.gens()
bounds = dict(aa1=2**BITS, aa2=2**BITS, bb1=2**BITS, bb2=2**BITS)

In [3]:
eq0 = x0**3 + a*x0 + b - y0**2
eq1 = (x1 + aa1)**3 + a*(x1 + aa1) + b - (y1 + bb1)**2
eq2 = (x2 + aa2)**3 + a*(x2 + aa2) + b - (y2 + bb2)**2


Note that the unknowns $a,b$ are large and so more annoying than the small $\alpha_1,\alpha_2,\beta_1,\beta_2$. Since we have 3 equations, we can get rid of $a, b$ and obtain a single equation:

In [4]:
def resultant(p1, p2, var):
p1 = p1.change_ring(QQ)
p2 = p2.change_ring(QQ)
var = var.change_ring(QQ)
r = p1.resultant(p2, var)
return r.change_ring(F)

poly = eq0
poly1 = resultant(poly, eq1, b)
poly2 = resultant(poly, eq2, b)
poly = resultant(poly1, poly2, a)
poly /= poly.coefficients()[0]
print(poly.monomials())
len(poly.monomials())

[aa1^3*aa2, aa1*aa2^3, aa1^3, aa1^2*aa2, aa1*aa2^2, aa2^3, aa2*bb1^2, aa1*bb2^2, aa1^2, aa1*aa2, aa2^2, aa2*bb1, bb1^2, aa1*bb2, bb2^2, aa1, aa2, bb1, bb2, 1]

Out[4]:
20

This is a degree-4 equation over $\mathbb{F}_p$ with 4 variables and very small solutions. Using full Coppersmith-like would be rather hard here. For such extreme imbalance - max term example is $|\alpha_1^3\alpha_2| \le 2^{32\cdot 4}= 2^{128}$ versus $p \approx 2^{3000}$ - we can try simpler techniques. Linearization is one of them. Indeed, we can consider this equation as linear equation with 19 variables of varying sizes. The total size of unknowns is:

In [5]:
bits = 0
for mono in poly.monomials():
mono_bits = RR(log(mono.change_ring(ZZ).subs(**bounds), 2))
print(mono, "%.2f" % mono_bits )
bits += mono_bits
bits, "total"

aa1^3*aa2 128.00
aa1*aa2^3 128.00
aa1^3 96.00
aa1^2*aa2 96.00
aa1*aa2^2 96.00
aa2^3 96.00
aa2*bb1^2 96.00
aa1*bb2^2 96.00
aa1^2 64.00
aa1*aa2 64.00
aa2^2 64.00
aa2*bb1 64.00
bb1^2 64.00
aa1*bb2 64.00
bb2^2 64.00
aa1 32.00
aa2 32.00
bb1 32.00
bb2 32.00
1 0.00

Out[5]:
(1408.00000000000, 'total')

The rule of thumb for a modular linear equation is to compare the total bit-length of unknowns (1408) to the bit-length of the modulus (3000). We have a large margin! Let's go LLL:

In [6]:
n = len(poly.monomials())
m = matrix(ZZ, n, n)
m[0] = vector(poly.coefficients())
m[1:,1:] = p * identity_matrix(n-1)
def prmat(m):
for row in m:
print(*[{0: "0", 1: "1", p: "p"}.get(v, "x") for v in row])
prmat(m)

1 x x x x x x 1 x x x x x x x x x x x x
0 p 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 p 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 p 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 p 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 p 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 p 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 p 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 p 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 p 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 p 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 p 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 p 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 p 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 p 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 p 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 p 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 p 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 p 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 p


What we are doing is simply finding a multiple of our polynomial modulo $p$ (by a constant) such that its absolute value will be less than $p$ even with maximum values of the variables.

In [7]:
# inspired by defund's code style https://github.com/defund/coppersmith
monos = vector(poly.change_ring(ZZ).monomials())
factors = [mono(**bounds) for mono in monos]

[m.rescale_col(i, factor) for i, factor in enumerate(factors)]
m = m.LLL()
m = m.change_ring(QQ)
[m.rescale_col(i, QQ(1)/factor) for i, factor in enumerate(factors)]
m = m.change_ring(ZZ)
m

Out[7]:
20 x 20 dense matrix over Integer Ring
In [8]:
polys = []
for pol in m*monos:
maxval = sum(
(abs(int(coef)) * mono).change_ring(ZZ).subs(**bounds)
for coef, mono in pol
)
print("polynomial with max value", RR(log(maxval, 2)), "bits")
if maxval < p:
polys.append(pol)
print("got", len(polys), "polynomials")

polynomial with max value 2835.33974382588 bits
polynomial with max value 2835.30655858310 bits
polynomial with max value 2835.39889680289 bits
polynomial with max value 2835.56026947599 bits
polynomial with max value 2835.27471060592 bits
polynomial with max value 2835.39148433409 bits
polynomial with max value 2835.40453160632 bits
polynomial with max value 2835.29853369271 bits
polynomial with max value 2835.67243387373 bits
polynomial with max value 2835.54980064626 bits
polynomial with max value 2835.39233844327 bits
polynomial with max value 2835.45721853969 bits
polynomial with max value 2835.60736646230 bits
polynomial with max value 2997.92022101675 bits
polynomial with max value 3061.29094246777 bits
polynomial with max value 3061.29094246777 bits
polynomial with max value 3091.96901437270 bits
polynomial with max value 3093.29094246767 bits
polynomial with max value 3093.29094246767 bits
polynomial with max value 3125.29094246767 bits
got 13 polynomials


13 good polynomials! That's more than enough. Note that we still have degree-4 equations, but over $\mathbb{Z}$ instead of $\mathbb{F}_p$. Out of a single equation! In fact, this trick can be also used if you start with an equation over integers: choose an appropriate-sized modulus and do the process to get more equations for free!

It is left to solve the polynomial system over integers. One standard way is to use GrÃ¶bner basis/resultats, but it's performance is unpredictable and unreliable. A better way is to perform bitwise recovery (modulo 2, modulo 4, modulo 8, etc.). It is especially powerful when we have more polynomials then equations.

In [9]:
def recover(hs, vars, solbits, padbits=20):
from itertools import product
sols = {(0,) * len(vars)}
polys = [h.change_ring(Zmod(2**nbits)) for h in hs]
for i in range(nbits):
print("bit", i, "/", nbits, ":", len(sols), "candidates")
sols2 = set()
mod = 2**i
polys = [h.change_ring(Zmod(2*mod)) for h in hs]
for bits in product(range(2), repeat=len(vars)):
for sol in sols:
sol2 = tuple(ss + bit*mod for ss, bit in zip(sol, bits))
if any(poly(*sol2) for poly in polys):
continue
sols = sols2
if not sols:
print("fail", i)
return
print("sols?", i, len(sols))
# TBD: automate adding pad bits to determine right sols by smallness
for sol in sols:
# fix signs
sol = [v if v < 2**(nbits-1) else (v-2**nbits) for v in sol]
# too large solution
if any(abs(v) >= 2**solbits for v in sol):
continue
# wrong solution
if any(poly(*sol) for poly in hs):
continue
yield sol


In fact, we can see that 4 equations (against 4 variables) are sufficient here:

In [10]:
R = PolynomialRing(ZZ, names='aa1, aa2, bb1, bb2')
polys = [R(pol) for pol in polys]
sols = list(recover(polys[:4], R.gens(), BITS))
sols[0]

bit 0 / 52 : 1 candidates
bit 1 / 52 : 2 candidates
bit 2 / 52 : 2 candidates
bit 3 / 52 : 2 candidates
bit 4 / 52 : 2 candidates
bit 5 / 52 : 2 candidates
bit 6 / 52 : 2 candidates
bit 7 / 52 : 2 candidates
bit 8 / 52 : 2 candidates
bit 9 / 52 : 2 candidates
bit 10 / 52 : 2 candidates
bit 11 / 52 : 2 candidates
bit 12 / 52 : 2 candidates
bit 13 / 52 : 2 candidates
bit 14 / 52 : 2 candidates
bit 15 / 52 : 2 candidates
bit 16 / 52 : 2 candidates
bit 17 / 52 : 2 candidates
bit 18 / 52 : 2 candidates
bit 19 / 52 : 2 candidates
bit 20 / 52 : 2 candidates
bit 21 / 52 : 2 candidates
bit 22 / 52 : 2 candidates
bit 23 / 52 : 2 candidates
bit 24 / 52 : 2 candidates
bit 25 / 52 : 2 candidates
bit 26 / 52 : 2 candidates
bit 27 / 52 : 2 candidates
bit 28 / 52 : 2 candidates
bit 29 / 52 : 2 candidates
bit 30 / 52 : 2 candidates
bit 31 / 52 : 2 candidates
bit 32 / 52 : 2 candidates
bit 33 / 52 : 2 candidates
bit 34 / 52 : 2 candidates
bit 35 / 52 : 2 candidates
bit 36 / 52 : 2 candidates
bit 37 / 52 : 2 candidates
bit 38 / 52 : 2 candidates
bit 39 / 52 : 2 candidates
bit 40 / 52 : 2 candidates
bit 41 / 52 : 2 candidates
bit 42 / 52 : 2 candidates
bit 43 / 52 : 2 candidates
bit 44 / 52 : 2 candidates
bit 45 / 52 : 2 candidates
bit 46 / 52 : 2 candidates
bit 47 / 52 : 2 candidates
bit 48 / 52 : 2 candidates
bit 49 / 52 : 2 candidates
bit 50 / 52 : 2 candidates
bit 51 / 52 : 2 candidates
sols? 51 2

Out[10]:
[1155305281, 1076744709, -62177711, -30259545]
In [11]:
from struct import pack
from Crypto.Cipher import AES
sol_aa1, sol_aa2, sol_bb1, sol_bb2 = sols[0]
parts = [
(int(x1)+int(sol_aa1))^int(x1),
(int(y1)+int(sol_bb1))^int(y1),
(int(x2)+int(sol_aa2))^int(x2),
(int(y2)+int(sol_bb2))^int(y2),
]
key = pack(">4I", *parts)
AES.new(key, mode=AES.MODE_ECB).decrypt(ct).decode()

Out[11]:
'Aero{"Voldemort," said Riddle softly, "is my past, present, and future, Harry Potter...."}\x00\x00\x00\x00\x00\x00'

Bonus: this approach works for 64-bit secrets!

In [ ]:
sola = F(31337) / F(1000)
solb = F(31333337) / F(1000000)

BITS = 64
solaa1 = randrange(2**BITS)
solbb1 = randrange(2**BITS)
solaa2 = randrange(2**BITS)
solbb2 = randrange(2**BITS)

E = EllipticCurve(GF(p), [sola, solb])
x0, y0 = E.random_element().xy()
x1, y1 = E.random_element().xy()
x2, y2 = E.random_element().xy()
x1 -= solaa1
x2 -= solaa2
y1 -= solbb1
y2 -= solbb2
solaa1, solaa2, solbb1, solbb2


If you clone this notebook and run with SageMath kernel (python mode, todo so replace sage.repl.ipython_kernel with ipykernel_launcher in the kernel.json file), rerun everything from the # SOLVING cell up to recovery of aa1, aa2, bb1, bb2.

Bonus 2: here is the simple Mystiz-style solution instead of the overcomplicated part.

In [ ]:
poly /= poly.constant_coefficient() # this is crucial
n = len(poly.monomials())
factors = [
int(mono.subs(**bounds))
for mono in poly.monomials() if mono != 1
]
fmat = diagonal_matrix(QQ, factors)
m = matrix(QQ, n, n)
m[0:n-1,:1] = matrix(ZZ, n-1, 1, poly.coefficients()[:-1])
m[n-1,0] = p
m[:n-1,1:] = ~fmat
prmat(m)
m = m.LLL()
m[:,1:] *= fmat
m = m.change_ring(ZZ)

for row in m:
if row[0] == 1:
row = -row
if row[0] == -1:
sol = list(row[-4:])
if poly(*(sol + [0, 0])) == 0:
print("sol", sol)