Problem with LLL reduction on truncated LCG schemes

Question

I am struggling to apply Freize et al. paper to break a truncated LCG.

A truncated LCG is a pseudo random generator that outputs the $n$ leading bits $y_i$ of $x_i$, where $(x_i)$ is such that $x_{i+1} = \alpha \cdot x_i + \beta\mod M$, and $\alpha$ is an integer, coprime with $M$. The first state $x_1$ is unknown.

The objective is to break the LCG, ie, from one or several consecutive outputs, find the current state $x_i$. We suppose that $\alpha$ and $M$ are known, and $\beta = 0$ (and is known too).

---

The theory

The paper states (pages 5 and 6) that their method solves systems of modular equations of the form $\displaystyle\sum_{j=1}^k a_{ij}x_j \equiv c_i \mod M$, for all $i \in \{1, \dots, k\}$.

Deriving a lattice from the vectors $a_i$, the authors apply LLL to find a small basis, of vectors $w_1, \dots, w_k$.
By multiplying the matrix $\begin{pmatrix} w_1\\ \vdots\\ w_k \end{pmatrix}$ on the right by $x$, they derive a new set of equations: $\displaystyle\sum_{j=1}^k w_{ij}x_j \equiv c_i' \mod M$, for all $i \in \{1, \dots, k\}$.

Here, i guess that the $c_i'$ are just obtained from the $c_i$ by multiplying (on the left) by the same matrix that allowed us to go from our original basis to our reduced basis.

Now, if the $c_i'$ are small enough, the equalities actually hold on $\mathbb{Z}$ and we can solve them: we have effectively solved the system.

---

In practice

Now, the authors talk a bit more about their method for a LCG, part 3 (pages 11 and 12).

We have the relations $\alpha^{k-1} x_{1} - x_k = 0$ (edit: fixed typo), so I guess that it implies that all my $c_i$ are $0$.

Now, if I want to obtain my $c_i'$, I should multiply $c$ on the left by the same matrix as the one used to change the basis of the lattice. But, given that I have $c = 0$, I should have $c' = 0$ as well, so we don't need to compute this matrix.

Furthermore, my lattice is, as descibed by the paper, $$L = \begin{pmatrix} M & 0 && \dots & 0\\ \alpha & -1 &0& \dots & 0\\ \alpha^2 & 0 & -1 & \dots & 0\\ \vdots &&&\ddots & \vdots\\ \alpha^{k-1} & 0 & 0 & \dots & -1\\ \end{pmatrix}$$ which is trivially nonsingular, and as such, the only solution to $Lx = 0$ is $x = 0$, which is obviously wrong.

Even worse, this method does not use the $y_i$, but it obviously should at some point, so I guess I'm missing something.

---

My questions

(Most probably related)

1. In my case, what are the $c_i$, and do I really obtain the $c_i'$ from them by the basis change?
2. How are the $y_i$ supposed to take part in this method?

Answer 1

The lattice basis $L$ generates every point that is a solution to the congruence $a^{i-1} x_1 - x_i \equiv 0 \pmod{m}$. By definition, the reduced basis $B = \mathsf{LLL}(L)$ also generates the same solutions; it just has shorter vectors.

Now suppose we have a set of outputs $\mathbf{y} = y_i$ of the generator, which correspond to the most significant bits of the LCG multiplied by $2^s$, and we want to find $\mathbf{z} = z_i$ such that $\mathbf{y} + \mathbf{z} = \mathbf{x}$, where $\mathbf{x} = x_i$ is the reconstructed state. The key insight here is $$ \begin{align} & L \cdot \mathbf{x} \equiv 0 \pmod{m} \\ & B \cdot \mathbf{x} \equiv 0 \pmod{m} \\ & B \cdot \mathbf{x} = m \cdot \mathbf{k} \quad \mbox{(for some unknown $\mathbf{k}$)} \\ & B \cdot (\mathbf{y} + \mathbf{z}) = m \cdot \mathbf{k} \\ & B \cdot \mathbf{z} = m \cdot \mathbf{k} - B \cdot \mathbf{y} \end{align} $$

Now, because the coefficients in the reduced basis $B$ are small, we know that all the elements of $\mathbf{k} = k_i$ are small, and we can further deduce $k_i = \lfloor (B \cdot \mathbf{y})_i / m \rceil$, or in other words, $m \cdot k_i$ is the closest multiple of $m$ in relation to each element of $B \cdot \mathbf{y}$. Once we know $k_i$, this is a simple matter of linear algebra. Here's an example with Sage:

p = previous_prime(2^32)
a = randint(2, p)
s = randint(0, p)
X = [ a^(i+1) * s % p for i in range(4) ]
Y = [ x - x % 2^16 for x in X ] # high 16 bits
Z = [ x % 2^16 for x in X ] # low 16 bits (secret)

L = matrix([
    [  p,  0,  0,  0 ],
    [a^1, -1,  0,  0 ],
    [a^2,  0, -1,  0 ],
    [a^3,  0,  0, -1 ]
])

B = L.LLL()

W1 = B * vector(Y)
W2 = vector([ round(RR(w) / p) * p - w for w in W1 ])

Z_ = list(B.solve_right(W2))
print Z_
print Z