How to use LLL to solve linear systems of equations

Table of Contents

Solving a single equation

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

$$3x + 5y = 2$$

The first step is to write this in matrix form:

$$ \begin{bmatrix}x & y & 1\end{bmatrix}\begin{bmatrix}3 \\ 5 \\ -2\end{bmatrix} = 0$$

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

$$\begin{bmatrix}3 \\ 5 \\ -2\end{bmatrix}$$

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]$

$$ \begin{bmatrix} 3 & 1 & 0 & 0\\ 5 & 0 & 1 & 0\\ -2 & 0 & 0 & 1 \end{bmatrix} $$

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

$$ \begin{bmatrix} a & b & c \end{bmatrix} \cdot \begin{bmatrix} 3 & 1 & 0 & 0 \\ 5 & 0 & 1 & 0 \\ -2 & 0 & 0 & 1 \end{bmatrix} = \begin{bmatrix} 3a+5b-2c & a & b & c \end{bmatrix} $$

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

$$33x+7y = 11$$

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.

$$ \begin{bmatrix} 33 & 1 & 0 & 0\\ 7 & 0 & 1 & 0\\ -11 & 0 & 0 & W \end{bmatrix} $$
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

$$33 -3 \cdot 7 - 11 = 1$$

.

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

$$ \begin{bmatrix} 33\cdot W & 1 & 0 & 0\\ 7 \cdot W& 0 & 1 & 0\\ -11 \cdot W& 0 & 0 & W \end{bmatrix} $$
....: 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.

$$3x + 5y -2 = 0 \mod 3 \Rightarrow 3x + 5y -2 = k\cdot 3 $$

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

$$ \begin{bmatrix} 3 & 1 & 0 & 0\\ 5 & 0 & 1 & 0\\ -2 & 0 & 0 & 1 \end{bmatrix} $$

We just add another vector with the modulus

$$ \begin{bmatrix} 3 \cdot W& 1 & 0 & 0\\ 5 \cdot W& 0 & 1 & 0\\ -2\cdot W & 0 & 0 & W \\ 3 \cdot W& 0 & 0 & 0 \end{bmatrix} $$

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.

$$ 13x - 17y + 0z = -33 \\ 18x + 0y - 3z = 63 $$

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

$$ \begin{bmatrix} x & y & z & 1 \end{bmatrix} \begin{bmatrix} 13 & 18 \\ -17 & 0 \\ 0 & -3 \\ -33 & 63 \end{bmatrix} = \begin{bmatrix} 0 & 0 \end{bmatrix} $$

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

$$ \begin{bmatrix} 13 & 18 & 1 & 0 & 0 & 0\\ -17 & 0 & 0 & 1 & 0 & 0\\ 0 & -3 & 0 & 0 & 1 & 0\\ -33 & 63 & 0 & 0 & 0 & 1 \end{bmatrix} $$

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

$$ \begin{bmatrix} 13 \cdot W & 18 \cdot W & 1 & 0 & 0 & 0\\ -17 \cdot W & 0 & 0 & 1 & 0 & 0\\ 0 & -3 \cdot W & 0 & 0 & 1 & 0\\ -33 \cdot W & 63 \cdot W & 0 & 0 & 0 & W \end{bmatrix} $$
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

$$ s \cdot k = h + r \cdot \text{flag} \mod n $$

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:

$$ \begin{bmatrix} k & flag & 1 & \_ \end{bmatrix} \cdot \begin{bmatrix} s \cdot W& 1 & 0 & 0 \\ r \cdot W& 0 & 1 & 0 \\ -h \cdot W & 0 & 0 & W \\ n \cdot W& 0 & 0 & 0 \end{bmatrix} $$
  • 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

$$ 137x + 69y + 42z = 0 \mod 1337 $$

The lattice representing it is

$$ \begin{bmatrix} 137 & 1 & 0 & 0\\ 69 & 0 & 1 & 0\\ 42 & 0 & 0 & 1 \\ 1337 & 0 & 0 & 0 \\ \end{bmatrix} $$

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

$$ \begin{bmatrix} 137 \cdot 2^{100} & 1 & 0 & 0 \\ 69 \cdot 2^{100}& 0 & 1 & 0\\ 42 \cdot 2^{100}& 0 & 0 & 1\\ 1337 \cdot 2^{100}& 0 & 0 & 0 \end{bmatrix} $$

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.

$$ \begin{bmatrix} 137 \cdot 2^{100} & \frac{1}{100} & 0 & 0\\ 69 \cdot 2^{100}& 0 & 1 & 0\\ 42 \cdot 2^{100}& 0 & 0 & 1\\ 1337 \cdot 2^{100}& 0 & 0 & 0 \end{bmatrix} $$
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

$$ \begin{bmatrix} flag & o0 & o1 & 1 & \_ & \_ \end{bmatrix} \cdot \begin{bmatrix} coef_0 & coef_1 & 1 & 0 & 0 & 0 \\ 137 & 0 & 0 & 1 & 0 & 0 \\ 0 & 137 & 0 & 0 & 1 & 0 \\ out_0 & out_1 & 0 & 0 & 0 & 1 \\ p & 0 & 0 & 0 & 0 & 0 \\ 0 & p & 0 & 0 & 0 & 0 \end{bmatrix} $$

this representing our system of equations

$$ \begin{cases} out_0 = \text{flag} \cdot \text{coef}_0 + 137 \cdot o_0 \mod p \\ out_1 = \text{flag} \cdot \text{coef}_1 + 137 \cdot o_1 \mod p \end{cases} $$
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.

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:

$$ \begin{bmatrix} coef_0 \cdot W & coef_1 \cdot W & 0 & 0 & 0 \\ 137 \cdot W & 0 & 1 & 0 & 0 \\ 0 & 137 \cdot W & 0 & 1 & 0 \\ out_0 \cdot W & out_1 \cdot W & 0 & 0 & W \\ p \cdot W & 0 & 0 & 0 & 0 \\ 0 & p \cdot W & 0 & 0 & 0 \end{bmatrix} $$
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

$$ \begin{bmatrix} 0 & 0 & 1 & 1 & 1 & W \end{bmatrix} $$
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!

Sagemath tips & tricks

Instead of manually constructing the matrix that represent the system we can do:

sage: P = PolynomialRing(ZZ, 'x,y,z')
sage: P.<x,y,z> = PolynomialRing(ZZ)
sage: mod = 13
sage: eqs = [
....:     13*x - 17*y + 33,
....:     18*x -3*z - 63
....: ]
sage: A,mons = Sequence(eqs).coefficients_monomials(sparse=False)
sage: A
[ 13 -17   0  33]
[ 18   0  -3 -63]
sage: mons
(x, y, z, 1)
sage: L = block_matrix([ # assuming i want to solve mod p
....:     [A.T,1],
....:     [mod,0]
....: ])
sage: L
[ 13  18|  1   0   0   0]
[-17   0|  0   1   0   0]
[  0  -3|  0   0   1   0]
[ 33 -63|  0   0   0   1]
[-------+---------------]
[ 13   0|  0   0   0   0]
[  0  13|  0   0   0   0]
sage: ww = diagonal_matrix([1000]*len(eqs) + [1]*3 + [1000])
sage: L *= ww
sage: L = L.LLL()
sage: L /= ww
sage: L
[ 0  0 -2  0  1  0]
[ 0  0 -3  0 -5  0]
[ 0  0  0 13  0  0]
[-1  0  0 -3  0  0]
[ 0  0 -1  5 -1  1] # this is the solution, x,y,z = -1,5,-1
[ 0  1 -1  0 -2  0]