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
murder = self.cast_spell(b'AVADA KEDAVRA')
# 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)
return cipher.encrypt(plaintext + padding)
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:
flag = file.read().strip()
ciphertext = encrypt(key, flag)
print(f'Ciphertext:\n{ciphertext.hex()}')
if __name__ == '__main__':
main()
Data:
Wizard's power:
0x27266a1284b761f793a529b9664693a6b1db36864a8664898d98d8010ec51afffaef2d1c79ab2078c1e0b289a1719d34b0a081ca325ba2b367017ee8e0824aaa9488409e76e923d0f0fc917ddfc1b0534e93a74246405dadfd1683f0dc31682eed0fa6b95fc235c845e16d2ef40463b7e746668dad82981fc4e05933aca410c65b36f89738f7d97502f6626c38f595338f3864638d8613fb74c16b63f3969a49ebd103ef354ed756c3cd5e67f1d2dbe5acdbc088bd6c1d503acef4ec59e4a09efac4729ca796ad25217fe74e7a0c7ef5a3e1fcd9eb9288fb89e842ef0b16642f7e84e27df4bcb623726e2c44ef46be07f9b5a5f92fe2c77d0de79fa6d46193b064207125d8935c2ff04f63e2f858e98d2518077dc58e13307f01d65ae953efd70980f3aeed320b7a6b66eb0c578dc3f05d426f412c4e3c7a9bcc68f27fe236fde41400371a39f53828824f5de3d5902cd3e7dcaee58b89c1a234188e391d412e7bc598d4f10b2bcb26aab7cd09194e80be046022ee8471
Wizard's magic:
X: 0x4f69c16e493693a9342b7b9bb0ae42d0e7425996151c631a5620ad4d7ed372d285f04df82975a379080af08322fd927cf8ea9702f4533b5981351dd12358ca61c6e34af3f902b13d2981cae7cd2416e910bfe2c90f3b64c02bf933b1d11ff8ee384f5b507d9a0b1be38b15e7db3ba700ed32122a0163ec66530911fe142e098187771f5a248627a1709800069d402c1a61a17bd053ff6541f981898df6670420b0f3ea6ebd6d01d24f6fe9e3091a7c36df0c1ad7596e8ae090c54528385422911ff8e1dc201ab56844bc92acabc495c442104de5c34a5fe6661b7213500f4ccef6e76b94c34d60772fcf0c1250e66812d85efcda4b323250f740d0d57d40697fedae0bdbaceb0582cd7c82a27610e7f9fba5ac8f84e006d034dcf481f2dc9338acabd51c28afc1b4fc1d0cd7c1cd5ad21109e8708e9458e5301868c68920ab1aefb1e9184383002e6c893f1793443f8388f3f1b6e1f4ecaecb4cf000f57f677156efdacb20d60a35b8337ee416957aadffdc889dce487
Y: 0x2393c955e5f44ab3ac052dafd83554222728393b8fd20630f3c4f122d8c86d872a7b692b3782f12a6e57352c46328f85fad2d73eaccfcbf7beb47f97da2148c4f3fc7d6751bf66569c97a64f732e7cc71767ba11b1419732035c0285ff6973fba230a3d315fba7820855208e07e6fc5bc4b2cc3868aacde3b8e9b2075bfe861b2ebccd9b6c836d85dc319290263961344d348fe8faa0aa3ebca76fe514a7b7ba313a8200727b22a8714f8bba9e5ec2c549e5bfb5857c050d1ff0471d4b01426c2ee583a7d4b6c30df5ad2d9a6902f574416f8d55ed192b29d521e5d23a54e5062400b539468dca9aa3dc558feac63c88fc696d42434ad9a83551e6860aeb6a4d84ca80713387fa8c56a1e473a82af63ec03a71202aba0ae46fb9a97132f4c92d332327e2c11b79008586b22d92d60c3155e88b5e1f9c193b363ca28f0990400afb6e8148458708b89c6023c0d5ebb746bcd754fe37a84ee4dea6baa273b2ef31a864e9586f01bae855cc0d6f6055b2546b2a918664bdb6
(On curve <None>)
Horcrux #1:
X: 0x2700ff7abe81679c770d3171c993c55da4a47cda3360e33d2b763e65517f307a39d401a256ee0634f3f1c6c244aad34422dccfe9ebd1803339972095eb2b0a57c82ee0db9ba94365cfb5270c4924c5c1ec00db2aff529aa923b113d2ca8ad6a6774fb7c655cd101ea63bc5a6ea0261f8da82d455219c7584d7de0757b2fda627c4684d3bac8f899c24178c7e0ecef6c226892b86043d3853cdb777889ce2901d8496bf0232dd000208eb2ce77e953c478551b1112ebf4b02f0086726210a50dc20ac08eb15d846084f9324b4f1f5ec73e3ef7e4a5207e04ebb866f731201e5f626084d1b61c158cfd0fdbaae8b8dd23ba599689c74a790933a3e77daa1c95fbde63b74381aff2b98f41cfebb7f1b220a4f2e8c3361734d7bd6648c720efc0a5f978917d8b84b4764e416762884f00104981e62d876d460bd8c1095cd755d8f31c5377e5ef935da77ff82e823b49d817a1f91bd4d155306173eb07efa4567d362e1cddbf0483873d5efc9286c36a58e5b3d995f3d01ca60
Y: 0x1bf9a696ef976644dd903fb148892044fe65dd16bf12966a6b6e43be4b3b52aaa348a76ed5a53bbb400a366c59c96cf88a9c6a713273a722c0c8cd6f42b7c6f5e1911d2d323780bce65be80eef4d375dd3425a9ff832e4166765c7687aaf3c6f3c1ae7c561c46c49f0075bb68d48b95be5fd2c94996a42ba255eb638ae5ce449064cc81831f2239dbb93f4944693ee42a507dc2878988928dc9ed5bb5cb3b3eba9ac414b4171e505d285edaf02d13c0748e47b17a498cb0c15883977a11cde98995a022999a555fb3f1a3b9226dc4122f3db6c86036b1b0e6ada10b23b8c89ddb590f2186191c0f1148a04a8e35c905d4053943554c58e8f0e2ccac42c2307532e2b8f55b74fccb8ab24a977c557d10bd7c8621bef3e7f326d00c3ff28a52ee85ec61bb8aef14f67ac80da5885a93f2840a5e44d9800b0352163667ba851bb15c4083955e1fba5cd9d6e7bf103a2bbb0fbc868136cf2871815762514c691e6352af18324d1ccdc21da3e1c4c6aa771aa9fccf343ad64b8
(On curve <None>)
Horcrux #2:
X: 0xa439486ae1dbd96ca0796947609757b9dc068d8ed8287f98261c0c77e0940d29e25d27e16b97bc6983be505b8295886a5e880664ecb2d9161036f5beb1dbdac08da0cc2797d574e0feac67bbcb43275ab9e9d936aa17fd9a0a13a2f95e111b5e0893c4d4c3afae3bfa4a60fc5d0673215dd876b8bc960fc765401135dc708562cbeb572fcefaef1fef51782f6e0b495c0799cf89a9f47ffec11b0e780fb41e80d75681f7c9dd5fbf5e93e1dcda8091ca6a84de5b82b3abc9194913119e429781c12c11ebc6b6e80c01683344319879dec8fe59dd2f9c5b7b6e5e211153ae3b5fd161edae59777f3e76ecdeaef518495e512119451c6a97a8f471c4349ba7df8c3e9c1ba67f7a4fef57551d58b8c87607dc830c8074eec50841a1e365be90b528e4f39b8e19e7617191ac5118bb1abc44739d65384d915cb72e226122965ebae1581f55ecb12e523b0e904b3c0a1b26cf506ca030a68fe40d34c0912272277cd0d605ab057dbb5423cca28629661ad163232e766e80949
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:
from sage.all import *
proof.arithmetic(False)
p = 0x27266a1284b761f793a529b9664693a6b1db36864a8664898d98d8010ec51afffaef2d1c79ab2078c1e0b289a1719d34b0a081ca325ba2b367017ee8e0824aaa9488409e76e923d0f0fc917ddfc1b0534e93a74246405dadfd1683f0dc31682eed0fa6b95fc235c845e16d2ef40463b7e746668dad82981fc4e05933aca410c65b36f89738f7d97502f6626c38f595338f3864638d8613fb74c16b63f3969a49ebd103ef354ed756c3cd5e67f1d2dbe5acdbc088bd6c1d503acef4ec59e4a09efac4729ca796ad25217fe74e7a0c7ef5a3e1fcd9eb9288fb89e842ef0b16642f7e84e27df4bcb623726e2c44ef46be07f9b5a5f92fe2c77d0de79fa6d46193b064207125d8935c2ff04f63e2f858e98d2518077dc58e13307f01d65ae953efd70980f3aeed320b7a6b66eb0c578dc3f05d426f412c4e3c7a9bcc68f27fe236fde41400371a39f53828824f5de3d5902cd3e7dcaee58b89c1a234188e391d412e7bc598d4f10b2bcb26aab7cd09194e80be046022ee8471
F = GF(p)
x0 = F(0x4f69c16e493693a9342b7b9bb0ae42d0e7425996151c631a5620ad4d7ed372d285f04df82975a379080af08322fd927cf8ea9702f4533b5981351dd12358ca61c6e34af3f902b13d2981cae7cd2416e910bfe2c90f3b64c02bf933b1d11ff8ee384f5b507d9a0b1be38b15e7db3ba700ed32122a0163ec66530911fe142e098187771f5a248627a1709800069d402c1a61a17bd053ff6541f981898df6670420b0f3ea6ebd6d01d24f6fe9e3091a7c36df0c1ad7596e8ae090c54528385422911ff8e1dc201ab56844bc92acabc495c442104de5c34a5fe6661b7213500f4ccef6e76b94c34d60772fcf0c1250e66812d85efcda4b323250f740d0d57d40697fedae0bdbaceb0582cd7c82a27610e7f9fba5ac8f84e006d034dcf481f2dc9338acabd51c28afc1b4fc1d0cd7c1cd5ad21109e8708e9458e5301868c68920ab1aefb1e9184383002e6c893f1793443f8388f3f1b6e1f4ecaecb4cf000f57f677156efdacb20d60a35b8337ee416957aadffdc889dce487)
y0 = F(0x2393c955e5f44ab3ac052dafd83554222728393b8fd20630f3c4f122d8c86d872a7b692b3782f12a6e57352c46328f85fad2d73eaccfcbf7beb47f97da2148c4f3fc7d6751bf66569c97a64f732e7cc71767ba11b1419732035c0285ff6973fba230a3d315fba7820855208e07e6fc5bc4b2cc3868aacde3b8e9b2075bfe861b2ebccd9b6c836d85dc319290263961344d348fe8faa0aa3ebca76fe514a7b7ba313a8200727b22a8714f8bba9e5ec2c549e5bfb5857c050d1ff0471d4b01426c2ee583a7d4b6c30df5ad2d9a6902f574416f8d55ed192b29d521e5d23a54e5062400b539468dca9aa3dc558feac63c88fc696d42434ad9a83551e6860aeb6a4d84ca80713387fa8c56a1e473a82af63ec03a71202aba0ae46fb9a97132f4c92d332327e2c11b79008586b22d92d60c3155e88b5e1f9c193b363ca28f0990400afb6e8148458708b89c6023c0d5ebb746bcd754fe37a84ee4dea6baa273b2ef31a864e9586f01bae855cc0d6f6055b2546b2a918664bdb6)
x1 = F(0x2700ff7abe81679c770d3171c993c55da4a47cda3360e33d2b763e65517f307a39d401a256ee0634f3f1c6c244aad34422dccfe9ebd1803339972095eb2b0a57c82ee0db9ba94365cfb5270c4924c5c1ec00db2aff529aa923b113d2ca8ad6a6774fb7c655cd101ea63bc5a6ea0261f8da82d455219c7584d7de0757b2fda627c4684d3bac8f899c24178c7e0ecef6c226892b86043d3853cdb777889ce2901d8496bf0232dd000208eb2ce77e953c478551b1112ebf4b02f0086726210a50dc20ac08eb15d846084f9324b4f1f5ec73e3ef7e4a5207e04ebb866f731201e5f626084d1b61c158cfd0fdbaae8b8dd23ba599689c74a790933a3e77daa1c95fbde63b74381aff2b98f41cfebb7f1b220a4f2e8c3361734d7bd6648c720efc0a5f978917d8b84b4764e416762884f00104981e62d876d460bd8c1095cd755d8f31c5377e5ef935da77ff82e823b49d817a1f91bd4d155306173eb07efa4567d362e1cddbf0483873d5efc9286c36a58e5b3d995f3d01ca60)
y1 = F(0x1bf9a696ef976644dd903fb148892044fe65dd16bf12966a6b6e43be4b3b52aaa348a76ed5a53bbb400a366c59c96cf88a9c6a713273a722c0c8cd6f42b7c6f5e1911d2d323780bce65be80eef4d375dd3425a9ff832e4166765c7687aaf3c6f3c1ae7c561c46c49f0075bb68d48b95be5fd2c94996a42ba255eb638ae5ce449064cc81831f2239dbb93f4944693ee42a507dc2878988928dc9ed5bb5cb3b3eba9ac414b4171e505d285edaf02d13c0748e47b17a498cb0c15883977a11cde98995a022999a555fb3f1a3b9226dc4122f3db6c86036b1b0e6ada10b23b8c89ddb590f2186191c0f1148a04a8e35c905d4053943554c58e8f0e2ccac42c2307532e2b8f55b74fccb8ab24a977c557d10bd7c8621bef3e7f326d00c3ff28a52ee85ec61bb8aef14f67ac80da5885a93f2840a5e44d9800b0352163667ba851bb15c4083955e1fba5cd9d6e7bf103a2bbb0fbc868136cf2871815762514c691e6352af18324d1ccdc21da3e1c4c6aa771aa9fccf343ad64b8)
x2 = F(0xa439486ae1dbd96ca0796947609757b9dc068d8ed8287f98261c0c77e0940d29e25d27e16b97bc6983be505b8295886a5e880664ecb2d9161036f5beb1dbdac08da0cc2797d574e0feac67bbcb43275ab9e9d936aa17fd9a0a13a2f95e111b5e0893c4d4c3afae3bfa4a60fc5d0673215dd876b8bc960fc765401135dc708562cbeb572fcefaef1fef51782f6e0b495c0799cf89a9f47ffec11b0e780fb41e80d75681f7c9dd5fbf5e93e1dcda8091ca6a84de5b82b3abc9194913119e429781c12c11ebc6b6e80c01683344319879dec8fe59dd2f9c5b7b6e5e211153ae3b5fd161edae59777f3e76ecdeaef518495e512119451c6a97a8f471c4349ba7df8c3e9c1ba67f7a4fef57551d58b8c87607dc830c8074eec50841a1e365be90b528e4f39b8e19e7617191ac5118bb1abc44739d65384d915cb72e226122965ebae1581f55ecb12e523b0e904b3c0a1b26cf506ca030a68fe40d34c0912272277cd0d605ab057dbb5423cca28629661ad163232e766e80949)
y2 = F(0x29cfc760529c48d68bc086a5b4403f1d4446db04abe243c99baf659b5e67cd6cdac9a658f273f4c682b9e13dda72aaa1ede42c69c98640f2fc58eabeef143a65334e4a236a6e72a157d92ab1a541c9bcf7d953b386a40d68880312ed1e900f6ba481bf134515f5a4c245a0d924db0e5105fc44fae71fc991df85219403ba5d48f68d5d91151f8411b4cbb6971bd63030b523989fa7ec94c14136ac712d4e91eca4c930d79caab7c328043c762917c4b3f868ce037cae9f19a579272c6b13ca8c19d349f5777bf6ac9c078d8472c92582daf96b30d4f7b8fdc004a36b792c133e2b6511956480892636ea91a3361afd8a3afbe3a5bf889feb5a5dc143f5a917347591634218066f2f71b36afd5257f6637152ac9a0965daa881ddc86ca8a8545b255cbb14e7738297a55c7428d9a0e79dee37f801fc1d49d205be809e0cd1b8a1b9d1d5b1b1a9ae98e7c564718cb17d10a3dd01c2914d7bb96c45a4fc942c2ed628354882464407c3208938282471aa2b8016e46c998e4)
BITS = 32
# 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)
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:
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]
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:
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
(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:
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.
# 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
20 x 20 dense matrix over Integer Ring
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.
def recover(hs, vars, solbits, padbits=20):
nbits = solbits + padbits
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
sols2.add(sol2)
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:
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
[1155305281, 1076744709, -62177711, -30259545]
from struct import pack
from Crypto.Cipher import AES
ct = bytes.fromhex("b6094505efd8e2824e8cc8698e5e68e3b2c306a2c8179fbefbfb7d2fc1a0cfb921a89f94dde88b04e655ad87d15efcc7466af652d6a330e92babceca94be63b126702153ad83d24baf7d8848bb71202771a2d80870fbd17462a715906dfc63bc")
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()
'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!
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.
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)