How to use LLL to solve linear systems of equations

Solving a single equation

Lets start with an easy one. Find (small) solutions to the equation

The first step is to write this in matrix form:

Solving the equation is the same as finding a linear combination of the vectors $[3], [5], [-2]$ that equals $[0]$, and finding a linear combination of vectors that gets you to a small vector is exactly what LLL is good at!

If we construct the lattice

Then the LLL algorithm will give us the shortest vector, and if the equation is solvable over $ZZ$ then we will see 0 in output.

sage: Matrix([
....:     [3],
....:     [5],
....:     [-2]
....: ]).LLL()
[ 0]
[ 0]
[-1]

This is actually useless information, as we want to know what are the values of x,y, not that we can get to 0 with a linear combination of $3$,$5$ and $-2$.

First intuition

Lets take a look at the lattice generated by the vectors $[3,1,0]$, $[5,0,1]$ and $[-2,0,0]$

In this case, if we find a linear combination of the 3 vectors, we also keep the coefficients used directly in the result vector.

So, applying LLL we see the linear combination that gets us to the smallest vector. Here there is a problem, as now the coefficients of the linear combination interfere with the original smallest vector we want to find (the solution to the equation), but we will get to it in a moment.

For now, lets see what happens if we try to apply LLL to this lattice:

sage: Matrix([
....:     [3, 1, 0, 0],
....:     [5, 0, 1, 0],
....:     [-2, 0, 0, 1]
....: ]).LLL()
[ 0 -1  1  1]
[-1 -1  0 -1]
[ 2  0  0 -1]

Lets analyze the output:

• The first vector (row) has a 0 in the first position, meaning that the result of our equation is 0. The coefficients are $a,b,c = -1,1,1$. So in this case we have $$3x + 5y = 2 \Rightarrow 3(-1) + 5(1) -2(1) = 0$$
• in the second row we see that the coefficients are $a,b,c=-1,0,-1$ $$3(-1) + 5(0) -2(-1) = -1 \Rightarrow -3 + 2 = -1$$
• same thing for the third row

We’ve found a solution with small $x,y$ values! from the output of LLL we then see that $x=-1$ and $y=1$ is a solution to $3x + 5y = 2$.

Actually this method doesn’t always work, as LLL doesn’t know that we are looking to solve an equation, it just finds the shortest vector.

Counter example

Lets try to solve

With our first lattice:

sage: Matrix([
....:     [33, 1, 0, 0],
....:     [7, 0, 1, 0],
....:     [-11, 0, 0, 1]
....: ]).LLL()
[ 0 -1  0 -3]
[ 1  1 -3  1]
[-3  0 -2 -1]

In the first row we have found a linear combination that gets us to 0, but the coefficients are $a,b,c = -1,0,-3$. We are only interested in linear combination where $c=1$ as this means that the constant term is used exactly once.

Using weights

Adding weight on the constant term

We have to tell LLL to use the constant term -11 exactly once. Remember that LLL finds a short vector. If we rescale the last column (representing the coefficient $c$) by a weight $W$ then LLL will pay a price (in terms of smallness of vector) of $W$ for every use of the constant term. If the weight is big enough then LLL will refuse to use it more than once.

sage: Matrix([
....:     [33 , 1, 0, 0],
....:     [7  , 0, 1, 0],
....:     [-11, 0, 0, 1000]
....: ]).LLL()
[  -2    1   -5    0]
[   5    1   -4    0]
[   1    1   -3 1000]

Seing 1000 in the last column means that the constant term -11 has been used only once (and this is exactly what we are looking for). Unfortunately for us, LLL didn’t manage to find a solution to the equation with $c=1$, as in the vector $[1,1,-3,1000]$ we have

.

Note: It found small x and y that gets to a small result using 11 exactly once, but we have to find a solution!

Adding weight on the equation-column

The problem is that LLL finds a short basis, and in this case some basis could be any linear combination of the vectors, like:

• $\begin{bmatrix} 1 & 1 &-3& 1000\end{bmatrix}$ or
• $\begin{bmatrix} 0 & -2 & 11 & 1000\end{bmatrix}$ (the one that we actually want to find)

The size of the vectors is (let’s define it as the sum of the element as absolute value for simplicity):

• $\begin{bmatrix} 1 & 1 &-3& 1000\end{bmatrix}$ has size 1005
• $\begin{bmatrix} 0 & -2 & 11 & 1000\end{bmatrix}$ has size 1013

So the problem is that LLL will prefer the vector $\begin{bmatrix} 1 & 1 &-3& 1000\end{bmatrix}$ with a 1 on the first element against a vector $\begin{bmatrix} 0 & -2 & 11 & 1000\end{bmatrix}$ with a 0 on the first element, as the overall size is smaller.

We can fix this, by rescaling the first column again by a big number W! Multiplying the first column by 1000 (or 1337, or 0xdeadbeef, any big number) will result in vectors like:

• $\begin{bmatrix} 1000 & 1 &-3& 1000\end{bmatrix}$ has size 2004
• $\begin{bmatrix} 0 & -2 & 11 & 1000\end{bmatrix}$ has size 1013

Now the one with a 0 in the first element is the smallest one!

With $1000$ rescaled on the first column LLL will prefer to chose big x and y to get 0 on equation-column over getting small x,y but having a non-zero value there.

Lets implement this with a lattice

....: Matrix([
....:     [33  * 1000, 1, 0, 0],
....:     [7   * 1000, 0, 1, 0],
....:     [-11 * 1000, 0, 0, 1000]
....: ]).LLL()
[   0    7  -33    0]
[1000    3  -14    0]
[   0   -2   11 1000]

We are interested in rows where we find a solution to the equation (0 as first element) and -11 is used exactly once (1000 as last element) The vector $\begin{bmatrix} 0 & -2 & 11 & 1000 \end{bmatrix}$ as these properties. We’ve found a linear combination that uses the last vector (the one with the constant term) exactly once and solves our equation.

The smallest possible solution is $x=-2$ and $y=11$.

Solving a single equation with a modulus

Adding a modulus to the equation just means adding another variable.

Also, we have no constraints on finding a solution that minimizes the use of the modulus (k*3), so there is no point in adding its corresponding 1 column to the lattice. From our original lattice

We just add another vector with the modulus

And with the same process as before, we get to our solution

sage: Matrix([
....:     [3 *1000, 1, 0, 0],
....:     [5 *1000, 0, 1, 0],
....:     [-2*1000, 0, 0, 1000],
....:     [3 *1000, 0, 0, 0]
....: ]).LLL()
[    0    -1     0     0]
[    0     0     3     0]
[-1000     0     1     0]
[    0     0     1  1000]

The vector $\begin{bmatrix} 0 & 0 & 1 & 1000 \end{bmatrix}$ is the one interesting that solves our equation $3x + 5y -2 = 0 \mod 3$ with $x=0$ and $y=1$.

Solving a system of equations

Now it should be easy to see how to solve a system of equations.

As before, the first step is writing it as a matrix:

Adding the 1s to save the result of the linear combination we get the lattice

And putting the weights to force a solution that uses the constant terms exactly once

sage: W = 10000
....: Matrix([
....:     [13*W , 18*W, 1, 0, 0, 0],
....:     [-17*W, 0   , 0, 1, 0, 0],
....:     [0    , -3*W, 0, 0, 1, 0],
....:     [33*W, -63*W, 0, 0, 0, W]
....: ]).LLL()
[     0      0     17     13    102      0]
[-10000      0     -4     -3    -24      0]
[     0      0      4      5      3  10000]
[     0 -30000      0      0      1      0]

Now we have to look for the row that has $\begin{bmatrix} 0 & 0 & a & b & c & W \end{bmatrix}$, as this is the one that solves our system of equations. In our case from the third row we extract the solution $x=4$, $y=5$ and $z=3$.

Advanced usage of weights

Usually when solving a equations each unknown variable has its own target weight. For example in CTF (or in HNP in biased nonce ecdsa) we might have an equation like

Where (in ecdsa case) $s,h,r$ are known, while $k$ and $\text{flag}$ are unknown. Also, $\text{flag}$ might be 25characters long (200bit) while $k$ 100bit.

The classic lattice that we learned to use doesn’t embed this information:

• This just says to LLL: Find a linear combination that gets the first element to 0 while using the last one exactly once.
• But we want to say : Find a linear combination that gets the first element to 0 while using the last one exactly once and the first vector ($k$) is taken about $2^{100}$ times while the second one ($\text{flag}$) about $2^{200}$.

To make an easy example, consider an equation like

The lattice representing it is

As always we have to rescale the equation-column by a big number to force it equal to 0 (solve the equation!)

Now, assume that i don’t want to find a small solution, but a particular one, where $x$ is about 100 times bigger than $y$ and $z$.

If we divide the x-column (second one) by 100, we are essentially saying to LLL that taking x 10 times is the same as taking y or z once, so naturally x will be about 10 times bigger than y and z.

sage: L = Matrix(QQ, [
....:     [137, 1, 0, 0],
....:     [69, 0, 1, 0],
....:     [42, 0, 0, 1],
....:     [1337, 0, 0, 0]
....: ]).LLL()
....: ww = diagonal_matrix([1<<100,1/100,1,1],sparse=False)
....: L *= ww
....: L = L.LLL()
....: L /= ww # rescale back just to display it better
....: L
[  0  67   2   1]
[  0 147   0  -2]
[  0 117  -1   2]
[  1 -68   0  -1]

LLL gave as solutions where x is about 100 times x and z, like $\begin{bmatrix} 117 & -1 & 2 \end{bmatrix}$

Toy problem / CTF / example

Let’s solve a more complex problem as we would do in a CTF:

import random
flag =  ???????????????? # 16 digits flag
p =     137137137137137137137137
outs = []
coeffs = []
for _ in range(2):
    randcoef = random.randint(0, p-1)
    offset = random.randint(0, 137137137137)
    coeffs.append(randcoef)
    outs.append((flag * randcoef + offset * 137) % p)
# sage: outs
# [24584806730324539451798, 74414422836370613031242]
# sage: coeffs
# [19748313136215130712968, 49984170649680941695978]

To solve this, we start by constructing the normal lattice, and then we rescale the columns

this representing our system of equations

P.<flag_sym, o0,o1> = PolynomialRing(ZZ)
offsets_sym = [o0,o1]
eqs = []
for i in range(2):
    eqs.append(outs[i] - (flag_sym * coeffs[i] + 137 * offsets_sym[i]))
A,mons = Sequence(eqs).coefficients_monomials(sparse=False)
L = block_matrix(QQ, [
    [A.T,1],
    [p,0]
])
# sage: mons
# (flag_sym, o0, o1, 1)
# sage: L.change_ring(Zmod(102))
# [80  2| 1  0  0  0]
# [67  0| 0  1  0  0]
# [ 0 67| 0  0  1  0]
# [86 38| 0  0  0  1]
# [-----+-----------]
# [43  0| 0  0  0  0]
# [ 0 43| 0  0  0  0]

Note: The code above is just a fancy and fast way to construct the lattice described above, we could have done it manually. To understand how this works go on the section An "automatic" way to LLL with sage (i don’t know how to link in markdown D: )

Now to apply the easy weights:

• The first 2 columns (equations) must result in a 0 (aren’t we solving equations after all?), so a big weight W is needed
• The last column (constant term) also must have a big weight, as we want to use the constant term exactly once

Trying with LLL on the rescaled matrix (below) we get:

W = 1<<100 # a 'big enough' number
ws = diagonal_matrix([W]*2 + [1]*3 + [W], sparse=False)
L *= ws
L = L.LLL()
L /= ws # rescale back just to display it better
L
# [0 0 -149866720202818 -156376078508744  157499495736488 1]

We see that the equation is solved and that we used the constant term exactly once. BUT we have multiple problems here:

• the flag is negative
• the offset terms $o_1$ and $o_2$ are bigger than the maximum possible generated by the random function (156376078508744 > 137137137137)

Clearly LLL didn’t give as the solution that we were trying to find.

We have to do some adjustments and tell LLL that:

• the flag is in the range (1,10000000000000000)
• the offset is in the range (0, 137137137137)

To do this, we could say that flag is about 72919 (10000000000000000//137137137137) bigger than offset, and so we can divide the flag column by that number. This will give us the correct result (try it!).

A better way is to make both the flag and the offset have the same ‘penalty’, by rescaling both of them to 1.

• The flag is in (1,10000000000000000)? Then divide the column by 10000000000000000.
• The offset is in (0,137137137137)? Then divide the column by 137137137137.

Ideally, when LLL finds the flag, it will be with a vector very close to

flag_expected = 10000000000000000
offset_expected = 137137137137
W = 1<<100 # a 'big enough' number
ws = diagonal_matrix([W]*2 + [1/flag_expected] + [1/offset_expected]*2 + [W], sparse=False)
L *= ws
L = L.LLL()
L /= ws # rescale back just to display it better
L
# [0 0 1337133713371337 102933766746 107880418507 1]

And 1337133713371337 doesn’t looks like a random number, so it is surely the flag!

An “automatic” way to LLL with sage

Note: This is a method for costructing lattices that i stole from genni that he stole from blupper that he either stole from cryptohack or crafted, so gg to them

At this point you should have a nice understanding of this problem, and after some time you will find out that crafting by hand the same (or very similar) lattice is a pain in the ass.

Fortunatly sage has some cool function that we can use to ‘cheat’ a bit.

1. First of all, we construct the equations over Multivariate Polynomial Ring over ZZ
2. Then we convert the linear equations to matrices
3. Then we add the identity matrix and modulus
4. we set weight
5. .LLL

This is basically what i do every CTF:

# step 1
P.<x,y,z> = PolynomialRing(ZZ)
eqs = []
eqs.append(
    something*x + something*y + something*z
)
...

# step 2
# !! If you have an old version of sagemath you'll need to use .coefficient_matrix instead
A,mons = Sequence(eqs).coefficients_monomials(sparse=False)

# step 3
L = block_matrix(QQ, [
    [A.T,1],
    [mod,0] # if there is a modulus
])

# step 4 
ws = diagonal_matrix([1]*len(eqs) + [mons[0]_size, mons[1]_size, ...] + [1])
L /= ws
# step 5
L = L.LLL()
L *= ws
# now profit
print(L)

The hard part of doing this is just checking that the equations are constructed correctly, and putting the correct sizes for each of the variables.

Example, template in action in CTFs

Lets solve the first part of my challenge from TRXCTF 2025, called babyDLP. You can find the detailed writeup on this blog, so this aims just to show how to use this kind “template”

What the code does is a variant of ecdsa, basically

1. Generates a secret key d, with bounds $[0,2^{195}]$
2. Exposes an oracle that calculates $$ s = \frac{h*k_1 + d*r}{k_0} \mod \text{order} $$ Where

• k_0, k_1 are unknown, with bounds $[0, 2^{32}]$ (so very small)
• s,h,r given
• d is the secret key, the same in every equation

We have a big set system of linear equations with small unknowns. This screams LLL!

For simplicity i’ll not rewrite the full code of the challenge, but just a proof of concept of the lattice needed to recover d.

This is our simplified babyDLP challenge code

order = 35578737292020366913860455505331441912070860931484986460229
d = randint(0,order)

def oracle():
    k1, k2 = randint(0,2**32), randint(0,2**32)
    h, r = randint(0,2**195),randint(0,2**195)
    s = Zmod(order)(h*k2 + d*r)/k1
    return h,s,r

Lets follow the template above.

The first hard step is to create the equations correctly (ignoring the modulus for now)

N_EQS = 10

# our unknowns are the secret key d and for each equation we have ki and ki+1
P = PolynomialRing(ZZ,'d,' + ','.join(f'k{i}' for i in range(N_EQS*2)))

d_sym = P.gen(0)
ks = list(P.gens())[1:]

eqs = []

for i in range(N_EQS):
    h,s,r = oracle()
    eq = h*ks[2*i] + d_sym*r - ks[2*i+1]*s
    eqs.append(eq)

# if we are running this locally, and i suggest to do so when crafting the LLL matrix, check that equations are correct
for eq in eqs:
    assert eq(...) % order == 0

Now that we have build the equation, what follows is extremely immediate

# just copy paste from the template
A,mons = Sequence(eqs).coefficients_monomials(sparse=False)

# step 3
L = block_matrix(QQ, [
    [A.T,1],
    [order,0]
])

# step 4 
# sage: mons
# (d, k0, k1, k2, k3, k4, k5, k6, k7, k8, k9, k10, k11, k12, k13, k14, k15, k16, k17, k18, k19)
ws = diagonal_matrix([1]*N_EQS + [2**195] + [2**32]*2*N_EQS, sparse=False)
L /= ws
# step 5
L = L.LLL()
L *= ws
# now profit

# in L[1] we have
# sage: L[1]
# (0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 17546655683529260215849600963843405506199689886176446684070, -1352514475, -2093690758, -1120079241, -3050620498, -4032608895, -2769969530, -3562303639, -2721011978, -2354803615, -1930728980, -283212853, -2173874162, -1016071909, -1072240519, -402135576, -4253596876, -2670941126, -2910651261, -2885679131, -2585751271)

# and it matches the size we expect!!! The first 10 are zeroes, meaning we solved the equation, then one big number followed by small 32bit numbers.
# In this case the sign is flipped and d is -17546655683529260215849600963843405506199689886176446684070 % order

Equation with a double modulus

Sometimes it might happen that we want to solve something like

Where obviously $n_0$ is bigger than $n_1$ (otherwise $n_1$ would be useless).

If we try to build the normal lattice

What LLL will do is use as many times as it wants both $n_0$ and $n_1$. However, this will not give us the result we are looking for.

Why standard construction doesn’t work

Lets try to solve

from random import randint
flag = 123 # secret

secret_coeffs = [flag] + [randint(0, 200) for _ in range(3)]
def poly(x):
    return (sum([c*x**i for i,c in enumerate(secret_coeffs)]) % 123457) % 12347

points = [101,202,137,314]
outs = [
    poly(points[0]),
    poly(points[1]),
    poly(points[2]),
    poly(points[3]),
]
print(outs)

At first glance it looks like we can build (omitting W on equation-column for readability)

Note: Remember that the first 3 columns represent our equations and the identity matrix added to the up-right is needed to extract the coefficient

# Here i'm building by hand just formore clarity on what is happening
L = Matrix([
    [101^0, 202^0, 137^0, 314^0, 1, 0, 0, 0],
    [101^1, 202^1, 137^1, 314^1, 0, 1, 0, 0],
    [101^2, 202^2, 137^2, 314^2, 0, 0, 1, 0],
    [101^3, 202^3, 137^3, 314^3, 0, 0, 0, 1],
    [-outs[0], -outs[1], -outs[2], -outs[3], 0, 0, 0, 1],
    [123457, 0, 0, 0, 0, 0, 0, 0],
    [0, 123457, 0, 0, 0, 0, 0, 0],
    [0, 0, 123457, 0, 0, 0, 0, 0],
    [0, 0, 0, 123457, 0, 0, 0, 0],
    [12347, 0, 0, 0, 0, 0, 0, 0],
    [0, 12347, 0, 0, 0, 0, 0, 0],
    [0, 0, 12347, 0, 0, 0, 0, 0],
    [0, 0, 0, 12347, 0, 0, 0, 0],
])
# rescaling by weights 
W = 10000
ww = diagonal_matrix([W]*4 + [1]*3 + [W], sparse=False)
L *= ww
L = L.LLL()
L /= ww
L
# [ 0  0  0  0  0  0  0  0]
# [ 0  0  0  0  0  0  0  0]
# [ 0  0  0  0  0  0  0  0]
# [ 0  0  0  0  0  0  0  0]
# [ 0  0  0  0  0  0  0  0]
# [ 0  0  0  0  1  0  0  0]
# [ 0  0  0  0  0  0 -1  0]
# [ 0  0  0  0  0 -1  0  0]
# [ 0  0 -1  0  0  0  0  0]
# [ 0  0  0  1  0  0  0  0]
# [ 0  0  0  0  0  0  0 -1]
# [ 0  1  0  0  0  0  0  0]
# [-1  0  0  0  0  0  0  0]

We definitely didn’t get the result we were looking for.

In this lattice we would expect a row in the form of

Where each $secret_coeffs_i$ is in the range $(0,200)$

But clearly we haven’t found what we where looking for.

If we think more about our lattice, we are essentially saying: Find the shortest possible vector, where

• We mostly care about the first 4 column being small (big W on equations column)
• You must use at most 1 time the $outs_i$ (the big W)
• You can use as much as you want the first modulus
• You can use as much as you want the second modulus

Remember that LLL will compute linear combinations of different rows to minimize the output. In our case, just by focusing on these 2 vectors

We can find a linear combination that gets us

How? The GCD of the two moduli is 1, so by computing xgcd:

Once we have access to the vector $\begin{bmatrix} 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \end{bmatrix}$ we can just use it to simplify the first column of the matrix (as using this vector is literally free), getting to something like

Same process for 2nd and 3rd column.

Its obvious that at this point our reduced lattice doesn’t contain any usefull information.

Trick to make it work

We can’t give permission to LLL to use for free the two moduli. What can we do then? Add weights to them!

Remember that an equation in the form of

Can be written, by adding 2 extra variable, as

Where $k_0$ and $k_1$ will have specific bounds.

Specifically,

• $k_0$ will be about $\text{something} / n_0$
• $k_1$ will be about $n_0 / n_1$

So we now have information about the size of the coefficients of the different moduli! With this we can fix our lattice by choosing carefully our weights.

The lattice (without any weights) will extremely similar as the one before

Alternatively, to use a smaller lattice, we can still give the first modulus for free, as it is true that it will be used as much as possible to reduce the equation value, and only give bounds to the external modulus (the smaller one)

From this

To this

And so our lattice becomes

Lets construct this in sagemath:

from random import randint
flag = 123

secret_coeffs = [flag] + [randint(0, 200) for _ in range(3)]
def poly(x):
    return (sum([c*x**i for i,c in enumerate(secret_coeffs)]) % 123457) % 12347

points = [101,202,137,314]
outs = [
    poly(points[0]),
    poly(points[1]),
    poly(points[2]),
    poly(points[3]),
]
# solve
P.<flag_sym,c0,c1,c2,k0,k1,k2,k3> = PolynomialRing(ZZ)
coeffs = [flag_sym,c0,c1,c2]
ks = [k0,k1,k2,k3]
eqs = []
for i in range(len(outs)):
    eqs.append(sum([(c)*points[i]**j for j,c in enumerate(coeffs)]) - ks[i]*12347 - outs[i])

A,mons = Sequence(eqs).coefficients_monomials(sparse=False)
L = block_matrix(QQ, [
    [A.T,1],
    [123457,0]
])
# we will aim to get
# (0, 0, 0, flag, c0, c1, c2, c3, c4, k0, k1, k2, 1)
# with sizes
# (0, 0, 0, 200, 200, 200, 200, 200, 200, 123457//12347, 123457//12347, 123457//12347, 1)

ws = diagonal_matrix([1]*len(outs) + [200]*len(coeffs) + [123457//12347]*len(ks) + [1], sparse=False)
L /= ws
L = L.LLL()
L *= ws
L
# [   0    0    0    0  123  150   22  128    1    2    4    3    1]
# [   1    1    1    1    1    0    0    0    0    0    0    0    0]
# [   0    0    0    0   13    0    0    0   10   10   10   10    0]
# [   2   -1    1   -1 -349    7 -177  125   -5   19    5   -3    1]
# [  -2    1    2   -2 -254 -260  311   29  -15    5   10    0    0]
# [   0   -1   -1    4 -264  319  106   16   -4    8   -2   -2   -1]
# [  -2   -1    3    0  236 -378 -104   23   -9   -4   -1   27   -1]
# [  -1    0    0    2  -18 -326   29  117  -24  -16   30   -7   -1]
# [  -3    0    1    3   26 -253  176 -104   13   13  -18   -1    2]
# [   1    0    1   -1 -539 -495  282  202    5  -15   13    1    3]
# [  -3    0    2    0   -4 -124   49  394   11    8    6  -33   -1]
# [   1   -1   -2    1  309   71  842  317    3   -5   18  -16   -3]
# [   1    0    1    0 -237  397  343 -576    7  -17   27   -8    0]

And here we found the vector we where looking for!!