Throughout the writeup, I will use following variables:
- $v_1$ to represent the first IV (of around 1024 bits),
- $c_1$ to represent the first ciphertext block (of around 1024 bits),
- $b_1$ to represent the first cleartext block (of around 256 bits),
- $v_2$ to represent the second IV (of around 1024 bits),
- $c_2$ to represent the second ciphertext block (of around 1024 bits),
- $b_2$ to represent the second cleartext block (of around 256 bits)
We know that:
$$ \begin{equation} \begin{aligned} b_1 v_1 &\equiv c_1 \pmod{s} \ b_2 v_2 &\equiv c_2 \pmod{s} \end{aligned} \end{equation} $$
This means that there exist $k_1$ and $k_2$ such that:
$$ \begin{equation} \begin{aligned} b_1 v_1 &= c_1 + k_1 s \ b_2 v_2 &= c_2 + k_2 s \end{aligned} \end{equation} $$
Moreover, as $b_i$ are around 256 bits, and $v_i$, $c_i$ and $s$ are all around 1024 bits, $k_i$ must be around 256 bits each.
We can rearrange a bit the equations (multiply by $k_2$ the first one, by $k_1$ the second one):
$$ \begin{equation} \begin{aligned} k_2 b_1 v_1 &= k_2 c_1 + k_1 k_2 s \ k_1 b_2 v_2 &= k_1 c_2 + k_1 k_2 s \end{aligned} \end{equation} $$ $$ \begin{equation} \begin{aligned} k_2 b_1 v_1 - k_2 c_1 &= k_1 k_2 s \ k_1 b_2 v_2 - k_1 c_2 &= k_1 k_2 s \end{aligned} \end{equation} $$
As both equations now have an equal side, we can equate the two left sides which eliminates $s$:
$$ \begin{equation} \begin{aligned} k_2 b_1 v_1 - k_2 c_1 &= k_1 b_2 v_2 - k_1 c_2 \end{aligned} \end{equation} $$
Rearrange it:
$$ \begin{equation} \begin{aligned} k_2 b_1 v_1 - k_2 c_1 - k_1 b_2 v_2 + k_1 c_2 = 0 \end{aligned} \end{equation} $$
Here, we are left with a polynomial equation in four unknowns $b_1$, $b_2$, $k_1$ and $k_2$, which are all of size around 256 bits, which is rather small in contrast with the other variables. This reminded me of the Coppersmith method, which one can use to find small roots for modular multivariate polynomials. Our only obstacle here is that our equation is not modular.
So, we can artificially derive two modular equations from the previous one:
$$ \begin{equation} \begin{aligned} k_2 b_1 v_1 - k_2 c_1 + k_1 c_2 &\equiv 0 \pmod{v_2} \ k_2 c_1 + k_1 b_2 v_2 - k_1 c_2 &\equiv 0 \pmod{v_1} \end{aligned} \end{equation} $$
And then solve each one independently thanks to https://github.com/defund/coppersmith/blob/master/coppersmith.sage
from sage.all import *
import itertools
def small_roots(f, bounds, m=1, d=None):
"""
From https://github.com/defund/coppersmith/blob/master/coppersmith.sage
"""
if not d:
d = f.degree()
if isinstance(f, Polynomial):
x, = polygens(f.base_ring(), f.variable_name(), 1)
f = f(x)
R = f.base_ring()
N = R.cardinality()
f /= f.coefficients().pop(0)
f = f.change_ring(ZZ)
G = Sequence([], f.parent())
for i in range(m+1):
base = N^(m-i) * f^i
for shifts in itertools.product(range(d), repeat=f.nvariables()):
g = base * prod(map(power, f.variables(), shifts))
G.append(g)
B, monomials = G.coefficient_matrix()
monomials = vector(monomials)
factors = [monomial(*bounds) for monomial in monomials]
for i, factor in enumerate(factors):
B.rescale_col(i, factor)
B = B.dense_matrix().LLL()
B = B.change_ring(QQ)
for i, factor in enumerate(factors):
B.rescale_col(i, 1/factor)
H = Sequence([], f.parent().change_ring(QQ))
for h in filter(None, B*monomials):
H.append(h)
I = H.ideal()
if I.dimension() == -1:
H.pop()
elif I.dimension() == 0:
roots = []
for root in I.variety(ring=ZZ):
root = tuple(R(root[var]) for var in f.variables())
roots.append(root)
return roots
return []
v1 = 77334118706903017905851926195263987587841484120181851399779256448325060189155911977501089325280674983095989544156755173494545303262624697596415130561761174465953634200352305362018548024598549719830340102724845313751014603734472416792237506113777451564809242091589971475634784093257635837072723122869129885932
c1 = 257104261975234631627222375337907931839828921690943038556732428659431453965165130399603542051107563669989455426276044027652689911749817113376633983554884957476620329368636979168444887111773860109363809455024744449960370327325680088240165843341500932156129947860719407510938514240940526429199918706135135939861
v2 = 131234136650380837706057619662923628669346245381102353435547099562754303149511448099676487675818787057104119774585188433414155519402100381035873028169247038406113043693796360128663673745038554651890218055784588014206393800740604944579404795940741495363516367584858695705692053428372220019692899450781991366569
c2 = 112577117961628933308321194754772054386995688274339955628781905548246145645609147849567245280247535022404103263635474946325986631356316768058370788262302973339634423193842290331471740062740390145631798508179711672577501190748346477783370465157133473854847696508585568209618409989217067141141209077484850543791
bounds = (2**256, 2**256, 2**256)
shift = 2**254
Rv1 = Integers(v1)
Rv2 = Integers(v2)
Pv1.<b2_v1, k1_v1, k2_v1> = PolynomialRing(Rv1)
Pv2.<b1_v2, k1_v2, k2_v2> = PolynomialRing(Rv2)
f = (k2_v1+shift)*c1 + (k1_v1+shift)*(b2_v1+shift)*v2 - (k1_v1+shift)*c2
g = (k2_v2+shift)*(b1_v2+shift)*v1 - (k2_v2+shift)*c1 + (k1_v2+shift)*c2
roots_v1 = small_roots(f, bounds)
roots_v2 = small_roots(g, bounds)
The variables need to be “shifted” because else, all zeroes would be a valid solution. We thus need to move away from this solution, so that the smallest one is the one we are interested in.
Once the roots of the polynomials are recovered, we just need to format things properly to print the flag:
def center(value, mod):
v = value % mod
if v > mod//2:
v -= mod
return v
real_b2, real_k1_v1, real_k2_v1 = [center(int(u), v1)+int(shift) for u in roots_v1[0]]
real_b1, real_k1_v2, real_k2_v2 = [int(u)+int(shift) for u in roots_v2[0]]
flag = int(real_b1).to_bytes(32) + int(real_b2).to_bytes(32)
print(flag)