Chapter 49: LLL — A Step-by-Step Tutorial with Four Applications#

The Lenstra–Lenstra–Lovász (LLL) algorithm (1982) is the workhorse of applied lattice theory. In a single \(O(n^4 \log B)\)-time procedure it converts an arbitrary basis of an integer lattice into a reduced basis whose first vector is provably no more than an exponential factor longer than the actual shortest lattice vector. That guarantee is, by the standards of NP-hard problems, almost laughably weak — and yet, applied carefully, LLL has broken or shaped almost every public-key cryptosystem from the late 1970s onwards: the Merkle–Hellman knapsack, Coppersmith’s RSA-with-stereotyped-messages attack, the Boneh–Venkatesan attack on the Hidden Number Problem (and hence biased ECDSA), and most recently it sits underneath every lattice-based cryptanalytic estimator for ML-KEM and ML-DSA.

This chapter is the deep dive behind W2 of the post-quantum course. We build LLL bottom-up:

  • Part I — Mathematics. Lattices, Gram–Schmidt orthogonalisation, and the Hadamard ratio — a single number that tells us how good a basis is.

  • Part II — Building blocks. The two pieces that compose LLL: size reduction and the Lovász swap condition.

  • Part III — Full algorithm. Pseudocode, decomposed Python, a worked 2D example, and the famous Hermite-factor bound.

  • Part IV — Four applications. Each is a complete, runnable end-to-end demo with the cryptanalytic punch line at the end.

Where this fits in the course

Read this chapter alongside W2 (the LLL lecture). It complements the upstream cryptanalysis-course Chapter 40 (lattice problems) by going deeper into the algorithm itself and then turning around to show four direct applications. The W1 knapsack key you generated at the end of Chapter 46 is exactly the input to Application A below.

Part I — Mathematical foundations#

49.1 Lattices and bases#

A lattice \(\mathcal L \subset \mathbb R^n\) is the set of all integer linear combinations of \(n\) linearly independent vectors \(b_1, \ldots, b_n\):

\[ \mathcal L \;=\; \mathcal L(B) \;=\; \Bigl\{ \sum_{i=1}^n c_i b_i \;:\; c_i \in \mathbb Z \Bigr\} \]

The matrix \(B \in \mathbb R^{n \times n}\) whose rows are \(b_1, \ldots, b_n\) is a basis of \(\mathcal L\). A lattice has infinitely many bases — any \(B' = U B\) with \(U \in \mathrm{GL}_n(\mathbb Z)\) (an integer matrix with \(\det U = \pm 1\)) generates the same lattice. So a lattice is an equivalence class of bases under unimodular row operations.

The determinant (or covolume) of the lattice is

\[ \det(\mathcal L) \;=\; |\det B|. \]

It is invariant under change of basis: every basis of \(\mathcal L\) has the same absolute determinant. Geometrically, \(\det \mathcal L\) is the volume of the fundamental parallelepiped spanned by any basis.

# Tiny lattice utilities, used throughout the tutorial.
import numpy as np
import math
from fractions import Fraction


def lattice_det(B):
    '''Absolute determinant (covolume) of a square integer/float basis.'''
    return abs(np.linalg.det(np.array(B, dtype=float)))


def shortest_via_enumeration_2d(B, R=10):
    '''Find the shortest non-zero vector of a 2D lattice by enumeration.

    Brute-force over c1, c2 in [-R, R].  Only honest for tiny lattices.
    '''
    B = np.array(B, dtype=float)
    best = None
    best_norm = float('inf')
    for c1 in range(-R, R + 1):
        for c2 in range(-R, R + 1):
            if c1 == 0 and c2 == 0:
                continue
            v = c1 * B[0] + c2 * B[1]
            n = float(np.linalg.norm(v))
            if n < best_norm:
                best_norm = n; best = v
    return best, best_norm


# Demo: a 2D lattice with a deliberately bad basis (almost parallel vectors).
B_bad = [[201.0, 37.0], [1648.0, 297.0]]
print('Bad basis B =', B_bad)
print('  det L =', lattice_det(B_bad))
print('  basis vector norms =',
      np.linalg.norm(B_bad, axis=1).tolist())
v_short, n_short = shortest_via_enumeration_2d(B_bad, R=20)
print('  shortest non-zero vector found by enumeration:',
      v_short.tolist(), 'with norm', round(n_short, 4))
Bad basis B = [[201.0, 37.0], [1648.0, 297.0]]
  det L = 1279.0000000000073
  basis vector norms = [204.37710243566914, 1674.5485958908448]
  shortest non-zero vector found by enumeration: [40.0, 1.0] with norm 40.0125

49.2 Gram–Schmidt orthogonalisation (GSO)#

LLL’s central trick is to measure how far a basis is from being orthogonal, then make it more orthogonal one step at a time. The measurement is the Gram–Schmidt orthogonalisation.

Definition — GSO of a basis

Given \(b_1, \ldots, b_n\), define iteratively:

\[ \mu_{i,j} \;=\; \frac{\langle b_i, b_j^* \rangle}{\langle b_j^*, b_j^* \rangle} \quad (j < i), \qquad b_i^* \;=\; b_i - \sum_{j<i} \mu_{i,j}\, b_j^*. \]

Then \(b_1^*, \ldots, b_n^*\) are pairwise orthogonal and span the same real space as \(b_1, \ldots, b_n\). The numbers \(\mu_{i,j}\) are the GSO coefficients.

The \(b_i^*\) are not lattice vectors in general (their components are typically irrational). They are an analytic tool: the squared norms \(\|b_i^*\|^2\) satisfy \(\prod_i \|b_i^*\|^2 = \det(\mathcal L)^2\) (Hadamard’s identity), so they decompose the determinant across the basis. A “good” basis has its \(\|b_i^*\|^2\) values not falling off too fast as \(i\) grows — that is the property LLL enforces.

def gram_schmidt(B):
    '''Compute the Gram-Schmidt basis B* and the coefficient matrix mu.

    Inputs
    ------
    B : array-like (n, m)
        Rows are basis vectors (n basis vectors in R^m, with m >= n).

    Returns
    -------
    Bs : (n, m) ndarray
        Orthogonalised basis vectors b_1^*, ..., b_n^*.
    mu : (n, n) ndarray
        Lower-triangular matrix of GSO coefficients (mu[i, j] = mu_{i, j}).
        The diagonal is set to 1 and the upper triangle to 0 by convention.
    '''
    B = np.array(B, dtype=float)
    n = B.shape[0]
    Bs = np.zeros_like(B)
    mu = np.zeros((n, n))                 # only the lower triangle is filled
    for i in range(n):
        Bs[i] = B[i].copy()
        for j in range(i):
            denom = float(Bs[j] @ Bs[j])
            mu[i, j] = float(B[i] @ Bs[j]) / denom if denom > 0 else 0.0
            Bs[i] = Bs[i] - mu[i, j] * Bs[j]
    return Bs, mu


# Sanity check: B*[i] is orthogonal to B*[j] for j < i.
B = [[2.0, 0.0], [3.0, 1.0]]
Bs, mu = gram_schmidt(B)
print('B  =', B)
print('B* =', Bs.tolist())
print('mu =', mu.tolist())
print('orthogonality check  <B*[0], B*[1]> =', float(Bs[0] @ Bs[1]))
B  = [[2.0, 0.0], [3.0, 1.0]]
B* = [[2.0, 0.0], [0.0, 1.0]]
mu = [[0.0, 0.0], [1.5, 0.0]]
orthogonality check  <B*[0], B*[1]> = 0.0

49.3 The Hadamard ratio#

How orthogonal is a basis? The standard scalar measure is the Hadamard ratio:

\[ H(B) \;=\; \left( \frac{\det(\mathcal L)}{\prod_{i=1}^n \|b_i\|} \right)^{1/n} \;\in\; (0, 1]. \]

It is \(1\) when \(B\) is perfectly orthogonal (Hadamard’s inequality is tight) and tends to \(0\) as \(B\) becomes more parallel (skewed). LLL pushes \(H(B)\) towards \(1\).

def hadamard_ratio(B):
    '''Return H(B) in (0, 1], a quality measure of the basis B.

    1.0 = perfectly orthogonal; near 0 = very skewed.
    '''
    B = np.array(B, dtype=float)
    n = B.shape[0]
    det = abs(np.linalg.det(B))
    norms_prod = np.prod(np.linalg.norm(B, axis=1))
    if norms_prod == 0: return 0.0
    return (det / norms_prod) ** (1.0 / n)


# Demo on the bad basis we built earlier.
print('Bad basis  : H =', round(hadamard_ratio(B_bad), 4))

# A perfectly orthogonal basis of the same lattice would have H = 1.
# Hint: the shortest vector found above plus another short independent vector
# would be near-orthogonal.  We will reach this end via LLL in §49.6.
Bad basis  : H = 0.0611

Part II — The two LLL building blocks#

49.4 Size reduction#

Once we have the GSO, we can reduce the size of a basis vector \(b_i\) without leaving the lattice by subtracting integer multiples of earlier basis vectors. Concretely:

Size-reduction step (single \(i\), single \(j < i\))

If \(|\mu_{i,j}| > 1/2\), replace

\[ b_i \;\leftarrow\; b_i - \lfloor \mu_{i,j} \rceil \, b_j, \]

where \(\lfloor \cdot \rceil\) is rounding to the nearest integer. After this, the new \(\mu_{i,j}\) has \(|\mu_{i,j}| \le 1/2\). The lattice is unchanged. The other GSO coefficients \(\mu_{i, k}\) for \(k < j\) also update accordingly.

A basis is size-reduced when every \(|\mu_{i,j}| \le 1/2\) for \(j < i\). Geometrically, every basis vector lies inside the “Voronoi-ish” region around its earlier orthogonal projection, so no vector overshoots.

def size_reduce(B, mu):
    '''Apply size reduction in place.  Mutates B and mu and returns them.

    Walks i from 1 to n-1; for each i, subtracts integer multiples of
    earlier basis vectors so that |mu_{i, j}| <= 1/2 for every j < i.
    The lattice is unchanged: only the choice of basis shifts.
    '''
    n = B.shape[0]
    for i in range(1, n):
        for j in range(i - 1, -1, -1):
            q = round(mu[i, j])
            if q != 0:
                B[i]  = B[i]  - q * B[j]
                # mu[j, j] is conventionally 0 in our representation, so the
                # vector update touches only positions strictly less than j;
                # the position-j entry decreases by exactly q (the +1 implicit
                # diagonal of b_j against itself).
                mu[i, :j] = mu[i, :j] - q * mu[j, :j]
                mu[i, j]  = mu[i, j]  - q
    return B, mu


# Demo: take the bad basis, size-reduce, observe shorter b_2.
B = np.array(B_bad, dtype=float)
Bs, mu = gram_schmidt(B)
print('before  norms =', np.round(np.linalg.norm(B, axis=1), 3).tolist(),
      '  H =', round(hadamard_ratio(B), 3))
B, mu = size_reduce(B, mu)
print('after   norms =', np.round(np.linalg.norm(B, axis=1), 3).tolist(),
      '  H =', round(hadamard_ratio(B), 3))
before  norms = [204.377, 1674.549]   H = 0.061
after   norms = [204.377, 40.012]   H = 0.395

49.5 The Lovász condition#

Size reduction alone is not enough — it cannot swap two vectors. To make any progress on the length of \(b_1\), LLL also needs a swap test:

Definition — Lovász condition (parameter \(\delta \in (1/4, 1)\))

The pair \((b_{k-1}, b_k)\) satisfies the \(\delta\)-Lovász condition iff

\[ \|b_k^*\|^2 \;\ge\; \bigl(\delta - \mu_{k, k-1}^2\bigr) \, \|b_{k-1}^*\|^2. \]

If the condition fails, swap \(b_{k-1}\) and \(b_k\) (and update the GSO on the fly). If it holds, leave them and move on.

The traditional textbook value is \(\delta = 3/4\), which gives the standard worst-case approximation factor \(2^{(n-1)/2}\). Larger \(\delta\) (\(\to 1\)) yields better-quality output at higher running cost.

Why this is the right condition

After size reduction, \(\|b_k\|^2 = \|b_k^*\|^2 + \mu_{k,k-1}^2 \|b_{k-1}^*\|^2\). If we swap \(b_{k-1}\) and \(b_k\), the new \(b_{k-1}^*\) becomes \(b_k\) itself (because nothing earlier than \(k-1\) was orthogonalised against it). The Lovász condition is exactly the condition that the new \(\|b_{k-1}^*\|\) is no more than \(1/\sqrt{\delta}\) times the old one — i.e. the swap shrinks the leading orthogonal vector by at least a constant factor. That gives the polynomial-time termination guarantee.

Part III — The full algorithm#

49.6 Pseudocode#

LLL(B, delta):
    n = number of basis vectors
    (B*, mu) = gram_schmidt(B)
    k = 1
    while k < n:
        size_reduce(B, mu)              # only on row k, suffices
        if Lovasz(B*, mu, k, delta):
            k = k + 1
        else:
            swap(B[k-1], B[k])
            (B*, mu) = gram_schmidt(B)  # refresh
            k = max(k - 1, 1)
    return B

We now translate this directly into Python.

def lll(B, delta=0.75, verbose=False):
    '''The Lenstra-Lenstra-Lovasz algorithm.

    Internally we keep B as a list-of-lists of **Python ints** (so it is
    exact even when entries exceed 2^53), and use Python floats only for
    the GSO scratch arrays.  That is enough for exact size reduction and
    a reliable Lovasz-condition comparison out to ~64-bit basis entries.
    For very-large-entry lattices a Fraction-based mu would be required.

    Parameters
    ----------
    B : array-like (n, m)
        Initial basis with n basis vectors in R^m (rows).
    delta : float in (1/4, 1)
        Quality / runtime trade-off.  3/4 is the textbook value.

    Returns
    -------
    B_red : list of lists of int
        A delta-LLL-reduced basis of the SAME lattice.
    swaps : int
        Number of swap operations performed.
    '''
    # Keep B as ints; if input was float-array, round.
    B = [[int(round(x)) for x in row] for row in (
            B.tolist() if isinstance(B, np.ndarray) else B)]
    n_rows = len(B)

    def gso():
        '''Recompute GSO Bs (float) and mu (float) from the current int B.'''
        Bs = [[float(x) for x in row] for row in B]
        mu = [[0.0] * n_rows for _ in range(n_rows)]
        for i in range(n_rows):
            for j in range(i):
                denom = sum(s * s for s in Bs[j])
                if denom > 0:
                    mu[i][j] = sum(B[i][c] * Bs[j][c] for c in range(len(B[i]))) / denom
                for c in range(len(Bs[i])):
                    Bs[i][c] -= mu[i][j] * Bs[j][c]
        return Bs, mu

    Bs, mu = gso()
    k = 1
    swaps = 0
    while k < n_rows:
        # 1. Size-reduce row k against all earlier rows (exact integer ops on B).
        for j in range(k - 1, -1, -1):
            q = round(mu[k][j])
            if q != 0:
                for c in range(len(B[k])):
                    B[k][c] -= q * B[j][c]
                for c in range(j):
                    mu[k][c] -= q * mu[j][c]
                mu[k][j] -= q

        # 2. Lovasz condition.
        lhs = sum(s * s for s in Bs[k])
        rhs = (delta - mu[k][k-1] ** 2) * sum(s * s for s in Bs[k-1])
        if lhs >= rhs:
            k += 1
        else:
            B[k], B[k-1] = B[k-1], B[k]
            swaps += 1
            Bs, mu = gso()
            k = max(k - 1, 1)
            if verbose:
                print(f'  swap, total swaps so far = {swaps}')
    return B, swaps


# Sanity: run on the bad 2D lattice.
B_red, swaps = lll(B_bad, delta=0.75, verbose=True)
B_red_np = np.array(B_red, dtype=float)
print('Reduced basis:', B_red)
print('  norms =', np.round(np.linalg.norm(B_red_np, axis=1), 3).tolist())
print('  H     =', round(hadamard_ratio(B_red_np), 3))
print('  swaps =', swaps)
  swap, total swaps so far = 1
  swap, total swaps so far = 2
Reduced basis: [[1, 32], [40, 1]]
  norms = [32.016, 40.012]
  H     = 0.999
  swaps = 2

49.7 Worked example — a deliberately bad 2D lattice#

We picked the basis \(b_1 = (201, 37)\), \(b_2 = (1648, 297)\) to be almost parallel. The Hadamard ratio of the input basis was \(\approx 0.001\) — terrible. After LLL we expect (a) much shorter basis vectors, (b) Hadamard ratio close to 1, and (c) the first output vector to coincide (up to sign) with the actual shortest non-zero lattice vector.

# Quantitative comparison.
B_in = np.array(B_bad, dtype=float)
B_red, _ = lll(B_in, delta=0.75)
B_red_np = np.array(B_red, dtype=float)
# Brute-force enumeration with a wide window so we definitely see the optimum.
v_short, _ = shortest_via_enumeration_2d(B_in, R=60)

print('       basis 1 vector    norm     basis 2 vector    norm     H')
print('  in: ', B_in[0].tolist(), '%7.2f' % np.linalg.norm(B_in[0]),
      ' ', B_in[1].tolist(), '%7.2f' % np.linalg.norm(B_in[1]),
      ' %.3f' % hadamard_ratio(B_in))
print('  out:', B_red[0], '%7.2f' % np.linalg.norm(B_red_np[0]),
      ' ', B_red[1], '%7.2f' % np.linalg.norm(B_red_np[1]),
      ' %.3f' % hadamard_ratio(B_red_np))
print()
print('shortest non-zero vector by enumeration:',
      v_short.tolist(), 'with norm', round(np.linalg.norm(v_short), 3))
print('LLL b_1 norm vs enumeration shortest:',
      f'{np.linalg.norm(B_red_np[0]):.3f}  vs  {np.linalg.norm(v_short):.3f}',
      '  -> LLL beat enumeration!' if np.linalg.norm(B_red_np[0]) <
                                     np.linalg.norm(v_short) - 1e-3 else '')
       basis 1 vector    norm     basis 2 vector    norm     H
  in:  [201.0, 37.0]  204.38   [1648.0, 297.0] 1674.55  0.061
  out: [1, 32]   32.02   [40, 1]   40.01  0.999

shortest non-zero vector by enumeration: [-1.0, -32.0] with norm 32.016
LLL b_1 norm vs enumeration shortest: 32.016  vs  32.016 

49.8 The Hermite-factor guarantee#

Theorem — LLL approximation factor

Let \(B\) be the output of LLL with parameter \(\delta = 3/4\). Then

\[ \|b_1\| \;\le\; 2^{(n-1)/2} \cdot \lambda_1(\mathcal L) \]

where \(\lambda_1(\mathcal L)\) is the actual shortest-vector length. The constant \(2^{(n-1)/2}\) is the worst-case Hermite factor.

In practice, on random lattices, LLL outputs a first vector within a factor \(\delta_{LLL}^n\) of \(\lambda_1\) where \(\delta_{LLL} \approx 1.0219\) empirically. That’s the root-Hermite factor. So for \(n = 100\) LLL typically achieves a factor \(\approx 8\), not \(2^{50}\) — much better than the worst case.

The running time is \(O(n^4 \log B_{\max})\) classical, where \(B_{\max}\) is the largest input entry.

Part IV — Four direct applications#

We now turn LLL on four cryptanalytically interesting problems. Each application is complete: setup, lattice construction, Python, and verification.

Application A — Breaking the Merkle–Hellman knapsack#

This is the bridge from W1 to W2 of the course. Given a low-density Merkle–Hellman public key \(a_1, \ldots, a_n\) and a ciphertext \(S\), we want to recover the plaintext bit-vector \(\mathbf x \in \{0, 1\}^n\) such that \(\sum x_i a_i = S\).

Lagarias–Odlyzko (1985) lattice

Build the \((n+1) \times (n+1)\) lattice with rows

\[ b_i \;=\; (\underbrace{0, \ldots, 0, 2}_{\text{position } i}, 0, \ldots, 0,\; 2 a_i) \quad (i = 1, \ldots, n), \qquad b_{n+1} \;=\; (1, 1, \ldots, 1,\; 2 S). \]

Any \(\mathbf c \in \{0, 1\}^n\) with \(\sum c_i a_i = S\) corresponds to a short lattice vector \(v = -\sum c_i b_i + b_{n+1}\) whose first \(n\) coordinates are \(\pm 1\) and last coordinate is \(0\). LLL finds it for densities \(d \lesssim 0.94\) (Coster–Joux–LaMacchia–Odlyzko–Schnorr–Stern 1992).

# A. Build a low-density Merkle-Hellman keypair, encrypt, then break with LLL.
import random as _r

def merkle_hellman_keygen(n=10, target_bits=40, seed=0):
    rng = _r.Random(seed)
    w = [rng.randrange(2, 10)]
    for _ in range(n - 1):
        w.append(sum(w) + rng.randrange(1, 10))
    m = max(sum(w) + rng.randrange(1, 100), 1 << target_bits)
    while True:
        r = rng.randrange(2, m)
        if math.gcd(r, m) == 1: break
    a = [(r * wi) % m for wi in w]
    return a, w, m, r

def encrypt(a, bits):
    return sum(b * ai for b, ai in zip(bits, a))

def lo_lattice(a, S):
    '''The (n+1)x(n+1) Lagarias-Odlyzko lattice for the knapsack (a, S).'''
    n = len(a)
    L = [[0] * (n + 1) for _ in range(n + 1)]
    for i in range(n):
        L[i][i] = 2
        L[i][n] = 2 * a[i]
    for j in range(n):
        L[n][j] = 1
    L[n][n] = 2 * S
    return L

def bits_from_short_vec(v, n):
    '''Recover (0/1)-bits from the +/-1 entries of a short LO-lattice vector.'''
    out = []
    for x in v[:n]:
        xi = int(round(x))
        if xi == 1:   out.append(0)
        elif xi == -1: out.append(1)
        else: return None
    return out


# Demo.
n = 10
a, w_secret, m_secret, r_secret = merkle_hellman_keygen(n=n, target_bits=40, seed=42)
plain = [_r.Random(7).randrange(2) for _ in range(n)]
S = encrypt(a, plain)
print(f'n = {n},  plaintext = {plain}')
print(f'density d = n / log2(max a) = {n / max(a).bit_length():.3f}')
print(f'ciphertext S = {S}')

L = np.array(lo_lattice(a, S), dtype=float)
L_red, swaps = lll(L, delta=0.75)
print(f'\nLLL on the {n+1}-dim LO lattice  ({swaps} swaps)')

recovered = None
for row in L_red:
    if abs(row[-1]) > 1e-6: continue            # need last coord = 0
    bits = bits_from_short_vec(row, n)
    if bits is None: bits = bits_from_short_vec(-row, n)
    if bits is not None: recovered = bits; break

print(f'recovered  = {recovered}')
print(f'truth      = {plain}')
assert recovered == plain
print('  -> Lagarias-Odlyzko via LLL succeeded.')
n = 10,  plaintext = [1, 1, 1, 1, 1, 1, 1, 1, 1, 1]
density d = n / log2(max a) = 0.250
ciphertext S = 5587552117795

LLL on the 11-dim LO lattice  (53 swaps)
recovered  = [1, 1, 1, 1, 1, 1, 1, 1, 1, 1]
truth      = [1, 1, 1, 1, 1, 1, 1, 1, 1, 1]
  -> Lagarias-Odlyzko via LLL succeeded.

Application B — Coppersmith’s small-root attack on RSA#

In 1996 Don Coppersmith showed that an integer polynomial \(f \in \mathbb Z[x]\) of degree \(d\) has all its small roots modulo \(N\) recoverable in polynomial time, provided the roots satisfy \(|x_0| \lesssim N^{1/d}\). The construction uses LLL on a lattice of polynomial-coefficient vectors.

The cleanest cryptographic application: stereotyped RSA messages. Suppose Alice sends \(c = m^e \bmod N\) with \(e = 3\) (low public exponent), and Eve already knows the high-order bits of \(m\). Write \(m = m_0 + x\) with \(|x| < 2^{B}\) small. Then \(f(x) = (m_0 + x)^3 - c\) has \(x\) as a small root mod \(N\), and Coppersmith’s theorem recovers it.

We give a textbook 1-polynomial Coppersmith demo (no Howgrave-Graham shifting). It works whenever \(|x| < N^{1/3}\), i.e. roughly when 1/3 of the message is unknown.

# B. Coppersmith small-root attack on RSA, e=3.
def coppersmith_small_root(f_coeffs, N, X):
    '''Recover an integer root |x0| < X of f(x) = 0 (mod N), where
    f_coeffs gives the coefficients in INCREASING degree (f_0, f_1, ..., f_d).

    Single-polynomial Coppersmith.  Builds the lattice
        rows  N * X^i e_i           for i = 0 .. d-1   (the "Hensel" rows)
              X^i * f(xX) coefficient row scaled
    Runs LLL, reads off a small integer-coefficient polynomial vanishing
    at x0, finds its real roots, and returns the integer one.
    '''
    f = list(f_coeffs)
    d = len(f) - 1                              # polynomial degree

    # Build (d+1)x(d+1) lattice.  Row i (for i < d) is X^i * N * e_i.
    # Last row is the coefficients of f(xX).
    L = [[0] * (d + 1) for _ in range(d + 1)]
    for i in range(d):
        L[i][i] = N * (X ** i)
    for i in range(d + 1):
        L[d][i] = f[i] * (X ** i)
    L = np.array(L, dtype=object)               # big-int safe

    L_red, _ = lll(L.astype(float), delta=0.75)

    # Read out the smallest-norm reduced row -> g(x) = sum (row[i] / X^i) x^i.
    # Pick the LLL-reduced row with smallest Euclidean norm.
    norms = [float(np.linalg.norm(r)) for r in L_red]
    g_scaled = L_red[int(np.argmin(norms))]
    g = [int(g_scaled[i]) // (X ** i) if X**i else 0 for i in range(d + 1)]
    # Numpy roots: pass coefficients in DECREASING order.
    g_np = np.array(g[::-1], dtype=float)
    if abs(g_np[0]) < 1e-9: g_np = g_np[1:]
    roots = np.roots(g_np) if len(g_np) > 1 else []
    candidates = [int(round(r.real)) for r in roots if abs(r.imag) < 1e-3]
    for x0 in candidates:
        if abs(x0) <= X:
            # check we hit the original f mod N
            val = sum(c * (x0 ** i) for i, c in enumerate(f)) % N
            if val == 0:
                return x0
    return None


# Demo: RSA with e=3 and a stereotyped message.
#
# Single-polynomial Coppersmith (no Howgrave-Graham shifts) has a TIGHT bound:
#       X  <  N^{2 / (d (d+1))}
# For d = 3 that is X < N^{1/6}.  We pick a 53-bit toy N so that 8-bit unknowns
# fit comfortably inside the bound (N^{1/6} ~= 460 > 256 = 2^8).  Real attacks
# extend this with shifted polynomials to reach X < N^{1/d - eps}.
p, q = 100000007, 100000037
N = p * q                            # ~10^16 (53-bit) toy modulus
e = 3
m_true = 1234567890                  # the secret message
c = pow(m_true, e, N)

# Eve knows the high bits of m: m = m_high + x with |x| < 2^8.
unknown_bits  = 8
unknown_bound = 1 << unknown_bits
m_high = (m_true >> unknown_bits) << unknown_bits   # zero out low bits
print(f'modulus N  = {N}      (~{N.bit_length()}-bit, N^(1/6) ~ {round(N ** (1/6))})')
print(f'true m     = {m_true}')
print(f'known high = {m_high}  (low {unknown_bits} bits unknown, |x| < {unknown_bound})')
print(f'ciphertext = {c}')

# f(x) = (m_high + x)^3 - c  (mod N) has x_true = m_true - m_high as small root.
# Expand: a + b*x + c1*x^2 + d*x^3 where a = m_high^3 - c, b = 3*m_high^2,
# c1 = 3*m_high, d = 1.
f0 = (m_high ** 3 - c) % N
f1 = (3 * m_high ** 2) % N
f2 = (3 * m_high) % N
f3 = 1
x0 = coppersmith_small_root([f0, f1, f2, f3], N, unknown_bound)
print(f'recovered low x = {x0}')
print(f'truth     low x = {m_true - m_high}')
assert x0 == m_true - m_high
print('  -> Coppersmith via LLL recovered the unknown low bits.')
modulus N  = 10000004400000259      (~54-bit, N^(1/6) ~ 464)
true m     = 1234567890
known high = 1234567680  (low 8 bits unknown, |x| < 256)
ciphertext = 1866831500483285
recovered low x = 210
truth     low x = 210
  -> Coppersmith via LLL recovered the unknown low bits.

Application C — Integer relation detection#

Given approximate real numbers \(\alpha_1, \ldots, \alpha_n\), an integer relation is a non-zero \(\mathbf c \in \mathbb Z^n\) with \(\sum c_i \alpha_i = 0\) (or \(\approx 0\) at floating-point precision). Finding such a relation is how a CAS guesses closed forms — the PSLQ algorithm is famous for this, but LLL works almost as well.

Construct the lattice with rows

\[ b_i \;=\; (\underbrace{0, \ldots, 0, 1}_{\text{position } i}, 0, \ldots, 0,\; \lceil M \alpha_i \rceil) \]

where \(M\) is a large precision multiplier. A short vector \(v = \sum c_i b_i\) in this lattice has the form \((c_1, \ldots, c_n,\; M \sum c_i \alpha_i)\). If a true relation exists, the last coordinate is small (hence \(v\) is short), and LLL finds it.

Demo: rediscover that \(\alpha = \sqrt 2\) satisfies \(\alpha^2 - 2 = 0\), given only the decimal expansion of \(\alpha\).

# C. Integer relation detection: discover that sqrt(2) satisfies x^2 = 2.
def integer_relation(alphas, M=10**12):
    '''Find a small integer relation among the real numbers in `alphas`.

    Returns (c_1, ..., c_n) such that sum c_i * alphas[i] is approximately 0,
    or None if LLL finds no plausibly-short relation.
    '''
    n = len(alphas)
    L = [[0] * (n + 1) for _ in range(n)]
    for i in range(n):
        L[i][i] = 1
        L[i][n] = int(round(M * alphas[i]))
    L = np.array(L, dtype=float)
    L_red, _ = lll(L, delta=0.75)
    # The shortest reduced row should have small last coord.
    candidates = sorted(L_red, key=lambda r: abs(r[-1]))
    best = candidates[0]
    return [int(x) for x in best[:n]]


# alpha = sqrt(2);  ask for a relation among (alpha^2, alpha, 1).
alpha = math.sqrt(2)
relation = integer_relation([alpha ** 2, alpha, 1], M=10 ** 12)
print(f'alpha = sqrt(2) ~= {alpha}')
print(f'integer relation among (alpha^2, alpha, 1):  {relation}')
print(f'check:  {relation[0]} * alpha^2 + {relation[1]} * alpha + {relation[2]} '
      f'= {relation[0] * alpha ** 2 + relation[1] * alpha + relation[2]:.2e}')
print()

# Same trick for the golden ratio:  phi^2 - phi - 1 = 0.
phi = (1 + math.sqrt(5)) / 2
relation = integer_relation([phi ** 2, phi, 1], M=10 ** 12)
print(f'phi = (1+sqrt 5)/2 ~= {phi}')
print(f'integer relation among (phi^2, phi, 1):       {relation}')
print(f'check:  {relation[0]} * phi^2 + {relation[1]} * phi + {relation[2]} '
      f'= {relation[0] * phi ** 2 + relation[1] * phi + relation[2]:.2e}')
alpha = sqrt(2) ~= 1.4142135623730951
integer relation among (alpha^2, alpha, 1):  [-1, 0, 2]
check:  -1 * alpha^2 + 0 * alpha + 2 = -4.44e-16

phi = (1+sqrt 5)/2 ~= 1.618033988749895
integer relation among (phi^2, phi, 1):       [-1, 1, 1]
check:  -1 * phi^2 + 1 * phi + 1 = 0.00e+00

Application D — Hidden Number Problem (biased ECDSA nonces)#

The most cryptanalytically dramatic LLL application: recovering an ECDSA secret key from biased nonces. This is exactly the bug behind the PlayStation 3 firmware-signing key compromise (Sony reused nonces, but a partial-bias variant breaks the same way).

Setup

ECDSA signs \(h_i = H(m_i)\) as \(r_i = (k_i G)_x \bmod n\) and \(s_i = k_i^{-1} (h_i + r_i \cdot d) \bmod n\). Solving for \(k_i\):

\[ k_i \;\equiv\; s_i^{-1} h_i \;+\; (s_i^{-1} r_i)\, d \pmod n. \]

If for some leak reason the most significant \(\ell\) bits of every \(k_i\) are zero — so \(k_i < B := n / 2^\ell\) — then for every signature

\[ t_i \cdot d \;-\; k_i \;\equiv\; -u_i \pmod n, \qquad t_i := s_i^{-1} r_i,\;\; u_i := s_i^{-1} h_i. \]

That is the Hidden Number Problem: recover \(d\) given many \((t_i, u_i)\) and the promise that \(t_i d + u_i\) has small residue \(k_i\) mod \(n\).

Build the \((m+1) \times (m+1)\) Boneh–Venkatesan lattice (rows = basis):

\[\begin{split} \begin{pmatrix} n & 0 & \cdots & 0 & 0 \\ 0 & n & \cdots & 0 & 0 \\ \vdots & & \ddots & & \vdots \\ 0 & 0 & \cdots & n & 0 \\ t_1 & t_2 & \cdots & t_m & B/n \end{pmatrix} \end{split}\]

A vector \((k_1, k_2, \ldots, k_m, d \cdot B / n) - (\text{shift by } \mathbf u)\) has all entries small. LLL on the lattice with target the shift recovers it (Babai’s nearest-plane reduces CVP to LLL, but for demo we use the simpler embedding).

# D. Hidden Number Problem -> recover a secret from biased modular hints.
#
# Setup (bias-only HNP).  Given an unknown d in [0, n) and many random
# multipliers t_i, suppose the residues
#       k_i := d * t_i mod n
# are guaranteed to be small (k_i < B = n / 2^ell for some leak parameter
# ell -- e.g. the top ell bits of every k_i are zero).  Recover d.
#
# In real ECDSA (with non-zero message-hash term u_i), the same lattice gets
# one more row via the Boneh-Venkatesan "embedding" trick; the analysis is
# unchanged.  We focus on the cleaner u=0 case so the lattice is explicit.

def egcd(a, b):
    if b == 0: return a, 1, 0
    g, x1, y1 = egcd(b, a % b); return g, y1, x1 - (a // b) * y1

def mod_inv(a, n):
    g, x, _ = egcd(a % n, n)
    if g != 1: raise ValueError('not invertible')
    return x % n

def is_probable_prime(N, k=20):
    if N < 2: return False
    if N % 2 == 0: return N == 2
    d_pp, r_pp = N - 1, 0
    while d_pp % 2 == 0: d_pp //= 2; r_pp += 1
    for _ in range(k):
        a = _r.randrange(2, N - 1)
        x = pow(a, d_pp, N)
        if x in (1, N - 1): continue
        for _ in range(r_pp - 1):
            x = pow(x, 2, N)
            if x == N - 1: break
        else:
            return False
    return True


# (1) Toy "group order" and secret.  We deliberately keep n at 32 bits so
#     that the LLL Gram-Schmidt scratch still fits in IEEE 754 double
#     precision -- the same algorithm scales to real ECDSA primes once you
#     swap the GSO for a Fraction-based or fpylll-based implementation.
n = (1 << 32) - 5
while not is_probable_prime(n): n += 2
print(f'toy group order n ~ 2^32  ({n.bit_length()} bits, prime)')

secret_d  = 0xCAFEBABE % n         # the hidden number
leak_bits = 8                      # top 8 bits of every k_i are zero
B_bound   = n >> leak_bits         # so 0 <= k_i < B  (about 2^24)
m_sigs    = 40                     # number of biased "signatures" --
                                   # large enough that the magic vector
                                   # is the unique shortest in L_BV.
two_l     = 1 << leak_bits         # the scaling 2^ell

# (2) Generate biased samples (t_i, k_i = d * t_i mod n, k_i < B).
#     We pick t_i so that d * t_i mod n is forced small -- exactly the
#     situation that arises from "all top bits of the nonce are zero".
ts, true_ks = [], []
rng_demo = _r.Random(2026)
d_inv = mod_inv(secret_d, n)
for i in range(m_sigs):
    k = rng_demo.randrange(1, B_bound)        # small target residue
    t = (d_inv * k) % n                       # so d * t == k mod n
    ts.append(t); true_ks.append(k)

# (3) Boneh-Venkatesan lattice (dim m+1, integer-only after scaling by 2^ell).
#       row i (i < m): (n * 2^ell) e_i
#       row m:         (t_1*2^ell, t_2*2^ell, ..., t_m*2^ell, 1)
#     The vector (k_1*2^ell, ..., k_m*2^ell, d) lies in this lattice and is
#     SHORT: each entry is below n.  Worst-case Gaussian-heuristic length is
#     roughly sqrt(m+1) * n, so the SHORT vector beats it once m * ell > log2 n.
m = m_sigs
L = [[0] * (m + 1) for _ in range(m + 1)]
for i in range(m):
    L[i][i] = n * two_l
for i in range(m):
    L[m][i] = ts[i] * two_l
L[m][m] = 1

L_red, swaps = lll(L, delta=0.75)
print(f'LLL on {m+1}-dim BV lattice ({swaps} swaps)')

# (4) Read off d.  The target lattice vector has the shape
#       (k_1 * 2^ell, k_2 * 2^ell, ..., k_m * 2^ell,  d)
#     with the first m entries of magnitude < 2^ell * B = n and the last entry
#     in [-n, n].  We scan reduced rows for a candidate d that makes EVERY
#     d * t_i mod n small.  Both sign conventions are tried.
def k_residue(d_cand, t):
    k = (d_cand * t) % n
    return min(k, n - k)              # centred mod n

best_d, best_score = None, n
for row in L_red:
    last = int(round(row[m]))
    for sign in (+1, -1):
        cand = (sign * last) % n
        if cand == 0: continue
        score = max(k_residue(cand, t) for t in ts)
        if score < best_score:
            best_d, best_score = cand, score

# The true secret has score < B by construction.  Anything else is spurious.
recovered_d = best_d if best_score <= B_bound else None
print(f'best d candidate     = {hex(best_d) if best_d else None}'
      f'  (max bias = {best_score} = ~{best_score / B_bound:.2f} * B)')
print(f'recovered (passes B) = {hex(recovered_d) if recovered_d else None}')
print(f'truth                = {hex(secret_d)}')
assert recovered_d == secret_d
print('  -> Hidden number recovered from biased modular hints.')
toy group order n ~ 2^32  (32 bits, prime)
LLL on 41-dim BV lattice (416 swaps)
best d candidate     = 0xcafebabe  (max bias = 15987859 = ~0.95 * B)
recovered (passes B) = 0xcafebabe
truth                = 0xcafebabe
  -> Hidden number recovered from biased modular hints.

49.9 Summary#

We built LLL bottom-up:

  1. Mathematics. Lattices, Gram–Schmidt, the Hadamard ratio.

  2. Building blocks. Size reduction (clean up a basis vector via integer combinations of earlier ones) and the Lovász swap test (force the leading orthogonal length to shrink by a constant factor).

  3. Algorithm. A 30-line Python implementation, polynomial-time, with the worst-case Hermite factor \(2^{(n-1)/2}\) but typically \(\approx 1.022^n\) in practice.

Then four direct applications:

  • A. Lagarias–Odlyzko: break low-density Merkle–Hellman knapsacks.

  • B. Coppersmith: recover unknown bits of an RSA plaintext when the public exponent is small.

  • C. Integer relation detection: rediscover algebraic identities from numerical data (the workhorse of computer-algebra exploration).

  • D. Hidden Number Problem: recover an ECDSA secret key from biased nonces — the same family of attacks that broke Sony’s PlayStation 3 firmware-signing key.

The same algorithm underlies all four attacks. That universality — plus the fact that LLL outputs an upper bound on \(\lambda_1\), not a solution to NP-hard SVP — is why it remains, four decades on, the single most important algorithm in applied lattice cryptanalysis.

49.10 References#

  1. Lenstra, A. K., Lenstra, H. W., Lovász, L. (1982). Factoring polynomials with rational coefficients. Mathematische Annalen 261, 515–534.

  2. Lagarias, J. C., Odlyzko, A. M. (1985). Solving low-density subset sum problems. J. ACM 32(1), 229–246.

  3. Coster, M. J., Joux, A., LaMacchia, B. A., Odlyzko, A. M., Schnorr, C.-P., Stern, J. (1992). Improved low-density subset sum algorithms. Computational Complexity 2, 111–128.

  4. Coppersmith, D. (1996). Finding a small root of a univariate modular equation. EUROCRYPT ‘96, LNCS 1070.

  5. Howgrave-Graham, N. A. (1997). Finding small roots of univariate modular equations revisited. Cryptography and Coding, LNCS 1355.

  6. Boneh, D., Venkatesan, R. (1996). Hardness of computing the most significant bits of secret keys in Diffie–Hellman and related schemes. CRYPTO ‘96, LNCS 1109.

  7. Nguyen, P. Q., Shparlinski, I. E. (2002). The insecurity of the Digital Signature Algorithm with partially known nonces. Journal of Cryptology 15(3).

  8. Galbraith, S. D. (2012). Mathematics of Public Key Cryptography. Cambridge University Press. — the modern textbook on lattices in cryptography. Chapter 17 covers LLL.

  9. Nguyen, P. Q., Vallée, B., eds. (2010). The LLL Algorithm: Survey and Applications. Springer.

  10. Albrecht, M. R., Player, R., Scott, S. (2015). On the concrete hardness of Learning with Errors. J. Mathematical Cryptology 9(3). — modern lattice-attack estimator behind every PQC parameter discussion.