LLL algorithm in Python
Polynomial-time basis reduction giving a 2^((n-1)/2)-approximate shortest vector.
- File
lll_algorithm.py, 101 lines- Level
- I0 not checked yet
- Toolchain
- Python 3.14, pytest 9
- Test framework
- pytest
- Run the tests
python3 -m pytest
# SPDX-License-Identifier: MIT
# Part of https://algorithm.tamnd.com/lll-algorithm/. The authors are listed in prov/82SX9W.jsonl.
from fractions import Fraction
from math import floor
from typing import Sequence, Union
Number = Union[int, Fraction]
Basis = list[list[Number]]
def _dot(u: Sequence[Number], v: Sequence[Number]) -> Number:
return sum((a * b for a, b in zip(u, v)), 0)
def _nearest(x: Fraction) -> int:
"""The integer nearest to x; a half is rounded up, so -1/2 <= x - result < 1/2."""
return floor(x + Fraction(1, 2))
def lll_reduce_counted(
basis: Sequence[Sequence[Number]], delta: Number = Fraction(3, 4)
) -> tuple[Basis, int]:
"""LLL-reduce a lattice basis with exact rational arithmetic.
Input: n >= 1 linearly independent rows of the same length m, with int or
Fraction entries, and a rational delta with 1/4 < delta < 1. Output: a pair
(reduced, swaps). The rows of `reduced` are a basis of the same lattice that
is size-reduced (|mu_ij| <= 1/2) and satisfies the Lovasz condition
delta * B_(k-1) <= B_k + mu_(k,k-1)^2 * B_(k-1) for k = 2..n, where B_i is the
squared length of the Gram-Schmidt vector b_i*. `swaps` counts line L6.
Cost: for an integer basis with squared lengths at most X, O(n^2 log X)
swaps and O(n^3 m log X) arithmetic operations on rationals whose sizes are
O(n log X) binary digits; one call to the Gram-Schmidt step costs O(n^2 m).
Raises ValueError for dependent rows, rows of different lengths or a bad delta.
"""
delta = Fraction(delta)
if not Fraction(1, 4) < delta < 1:
raise ValueError("delta must satisfy 1/4 < delta < 1")
b: Basis = [list(row) for row in basis]
n = len(b)
if n == 0:
return [], 0
m = len(b[0])
if m == 0 or any(len(row) != m for row in b):
raise ValueError("rows must have the same positive length")
# L1: Gram-Schmidt data, mu[i][j] for j < i and B[i] = ||b_i*||^2.
mu = [[Fraction(0)] * n for _ in range(n)]
big_b: list[Fraction] = []
star: list[list[Fraction]] = []
for i in range(n): # @line L1
v = [Fraction(x) for x in b[i]]
for j in range(i):
mu[i][j] = Fraction(_dot(b[i], star[j])) / big_b[j]
v = [x - mu[i][j] * y for x, y in zip(v, star[j])]
norm = _dot(v, v)
if norm == 0:
raise ValueError("the rows are linearly dependent")
star.append(v)
big_b.append(norm)
def reduce(k: int, j: int) -> None:
"""Subtract the nearest integer multiple of b_j from b_k (0-based)."""
q = _nearest(mu[k][j])
if q != 0:
b[k] = [x - q * y for x, y in zip(b[k], b[j])]
for i in range(j):
mu[k][i] -= q * mu[j][i]
mu[k][j] -= q
swaps = 0
k = 1 # @line L2
while k < n: # @line L3
reduce(k, k - 1) # @line L4
if delta * big_b[k - 1] > big_b[k] + mu[k][k - 1] ** 2 * big_b[k - 1]: # @line L5
# L6: exchange b_(k-1) and b_k and update the Gram-Schmidt data.
swaps += 1
u = mu[k][k - 1]
new_prev = big_b[k] + u * u * big_b[k - 1]
mu[k][k - 1] = u * big_b[k - 1] / new_prev
big_b[k] = big_b[k - 1] * big_b[k] / new_prev
big_b[k - 1] = new_prev
b[k - 1], b[k] = b[k], b[k - 1] # @line L6
for j in range(k - 1):
mu[k - 1][j], mu[k][j] = mu[k][j], mu[k - 1][j]
for i in range(k + 1, n):
t = mu[i][k]
mu[i][k] = mu[i][k - 1] - u * t
mu[i][k - 1] = t + mu[k][k - 1] * mu[i][k]
k = max(1, k - 1) # @line L7
else: # @line L8
for j in range(k - 2, -1, -1): # @line L9
reduce(k, j)
k += 1 # @line L10
return b, swaps # @line L11
def lll_reduce(basis: Sequence[Sequence[Number]], delta: Number = Fraction(3, 4)) -> Basis:
"""LLL-reduce a basis; the same as lll_reduce_counted without the swap count."""
return lll_reduce_counted(basis, delta)[0]
test_lll_algorithm.py 212 lines
# SPDX-License-Identifier: MIT
# Part of https://algorithm.tamnd.com/lll-algorithm/. The authors are listed in prov/82SX9W.jsonl.
import itertools
import random
from fractions import Fraction
from math import floor, isqrt
import pytest
from lll_algorithm import lll_reduce, lll_reduce_counted
HALF = Fraction(1, 2)
def dot(u, v):
return sum((a * b for a, b in zip(u, v)), 0)
def gram(b):
return [[Fraction(dot(x, y)) for y in b] for x in b]
def inverse(a):
"""Exact inverse of a square Fraction matrix by Gauss-Jordan elimination."""
n = len(a)
m = [list(row) + [Fraction(int(i == j)) for j in range(n)] for i, row in enumerate(a)]
for c in range(n):
p = next(r for r in range(c, n) if m[r][c] != 0)
m[c], m[p] = m[p], m[c]
piv = m[c][c]
m[c] = [x / piv for x in m[c]]
for r in range(n):
if r != c and m[r][c] != 0:
f = m[r][c]
m[r] = [x - f * y for x, y in zip(m[r], m[c])]
return [row[n:] for row in m]
def det(a):
"""Exact determinant by fraction-valued Gaussian elimination."""
m = [[Fraction(x) for x in row] for row in a]
n, d = len(m), Fraction(1)
for c in range(n):
p = next((r for r in range(c, n) if m[r][c] != 0), None)
if p is None:
return Fraction(0)
if p != c:
m[c], m[p] = m[p], m[c]
d = -d
d *= m[c][c]
for r in range(c + 1, n):
f = m[r][c] / m[c][c]
m[r] = [x - f * y for x, y in zip(m[r], m[c])]
return d
def gso(b):
"""Gram-Schmidt from scratch: returns (mu, B) with B[i] = ||b_i*||^2."""
n = len(b)
star, mu, big = [], [[Fraction(0)] * n for _ in range(n)], []
for i in range(n):
v = [Fraction(x) for x in b[i]]
for j in range(i):
mu[i][j] = Fraction(dot(b[i], star[j])) / big[j]
v = [x - mu[i][j] * y for x, y in zip(v, star[j])]
star.append(v)
big.append(dot(v, v))
return mu, big
def same_lattice(b, c):
"""True if the rows of c are a basis of L(b): c = U b with U integer, det U = +-1."""
if len(b) != len(c):
return False
gi = inverse(gram(b))
u = []
for row in c:
r = [Fraction(dot(row, x)) for x in b]
coeff = [sum((r[j] * gi[j][i] for j in range(len(b))), Fraction(0)) for i in range(len(b))]
if any(x.denominator != 1 for x in coeff):
return False
combo = [sum((coeff[i] * b[i][t] for i in range(len(b))), Fraction(0)) for t in range(len(row))]
if combo != [Fraction(x) for x in row]:
return False
u.append(coeff)
return abs(det(u)) == 1
def is_lll_reduced(b, delta):
mu, big = gso(b)
n = len(b)
if any(abs(mu[i][j]) > HALF for i in range(n) for j in range(i)):
return False
return all(delta * big[k - 1] <= big[k] + mu[k][k - 1] ** 2 * big[k - 1] for k in range(1, n))
def shortest_sq(b):
"""lambda_1^2 by enumeration of all coefficient vectors in a box that is proved large enough.
Any v = x b with ||v||^2 <= R2 = min ||b_i||^2 has |x_i| <= sqrt(R2 * (G^-1)_ii),
by Cauchy-Schwarz with the dual basis vector d_i, whose squared length is (G^-1)_ii.
"""
g = gram(b)
gi = inverse(g)
r2 = min(g[i][i] for i in range(len(b)))
box = [isqrt(floor(r2 * gi[i][i])) for i in range(len(b))]
best = r2
for x in itertools.product(*[range(-t, t + 1) for t in box]):
if any(x):
v = [sum((x[i] * b[i][t] for i in range(len(b))), 0) for t in range(len(b[0]))]
best = min(best, dot(v, v))
return best
def potential(b):
"""D = product over i < n of d_i, d_i = product of B_1..B_i (an integer for integer b)."""
_, big = gso(b)
d, acc = Fraction(1), Fraction(1)
for i in range(len(b) - 1):
acc *= big[i]
d *= acc
return d
def test_worked_example():
out, swaps = lll_reduce_counted([[1, 1, 1], [-1, 0, 2], [3, 5, 6]])
assert out == [[0, 1, 0], [1, 0, 1], [-2, 0, 1]]
assert swaps == 2
assert gso(out)[1] == [1, 2, Fraction(9, 2)]
def test_dimension_one():
assert lll_reduce([[-7, 3]]) == [[-7, 3]]
assert lll_reduce([[5]]) == [[5]]
def test_empty_basis():
assert lll_reduce([]) == []
def test_already_reduced_is_unchanged():
for b in ([[1, 0, 0], [0, 1, 0], [0, 0, 1]], [[0, 1, 0], [1, 0, 1], [-2, 0, 1]], [[2, 0], [-1, 2]]):
out, swaps = lll_reduce_counted(b)
assert out == b and swaps == 0
def test_nearly_dependent():
b = [[1, 0, 0], [1000, 1, 0], [999, 1, 1]]
out = lll_reduce(b)
assert same_lattice(b, out) and is_lll_reduced(out, Fraction(3, 4))
c = [[100001, 100000], [100000, 99999]]
out = lll_reduce(c)
assert same_lattice(c, out) and is_lll_reduced(out, Fraction(3, 4))
assert sorted(dot(v, v) for v in out) == [1, 1]
def test_fraction_input_and_delta():
b = [[Fraction(1, 2), Fraction(1, 3)], [Fraction(5, 7), 2]]
for delta in (Fraction(1, 2), Fraction(3, 4), Fraction(99, 100)):
out = lll_reduce(b, delta)
assert same_lattice(b, out) and is_lll_reduced(out, delta)
def test_errors():
with pytest.raises(ValueError):
lll_reduce([[1, 2], [2, 4]])
with pytest.raises(ValueError):
lll_reduce([[1, 2], [3]])
with pytest.raises(ValueError):
lll_reduce([[1, 0], [0, 1]], Fraction(1, 4))
with pytest.raises(ValueError):
lll_reduce([[1, 0], [0, 1]], 1)
def random_basis(rng, n, m, size):
while True:
b = [[rng.randint(-size, size) for _ in range(m)] for _ in range(n)]
if det(gram(b)) != 0:
return b
def test_random_against_oracles():
rng = random.Random(20251106)
for trial in range(150):
n = rng.randint(2, 6)
m = n + rng.randint(0, 2)
size = rng.choice([3, 20, 200])
b = random_basis(rng, n, m, size)
if trial % 5 == 0: # a skewed basis: add big multiples of the first row
b = [b[0]] + [[x + rng.randint(-50, 50) * y for x, y in zip(row, b[0])] for row in b[1:]]
delta = rng.choice([Fraction(1, 2), Fraction(3, 4), Fraction(99, 100)])
out, swaps = lll_reduce_counted(b, delta)
assert same_lattice(b, out)
assert is_lll_reduced(out, delta)
alpha = 1 / (delta - Fraction(1, 4))
n1 = dot(out[0], out[0])
assert n1 <= alpha ** (n - 1) * shortest_sq(out)
assert n1 ** n <= alpha ** (n * (n - 1) // 2) * det(gram(b))
# each swap multiplies the integer potential by less than delta
assert potential(out) <= delta ** swaps * potential(b)
assert potential(out) >= 1
def test_dimension_two_with_delta_near_one_finds_minimum():
rng = random.Random(7)
for _ in range(100):
b = random_basis(rng, 2, 2, 30)
out = lll_reduce(b, Fraction(999, 1000))
# with delta close to 1, a reduced 2D basis has ||b1||^2 <= ||b2||^2 / delta
assert dot(out[0], out[0]) * Fraction(999, 1000) <= dot(out[1], out[1])
assert dot(out[0], out[0]) <= Fraction(1000, 999) * shortest_sq(out)
What the next level needs
- the tests did not run in this build