Structured Gaussian elimination in Python
Eliminates sparse columns and rows first to shrink a large sparse system before dense solving, used in integer factorization.
- File
structured_gaussian_elimination.py, 134 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/structured-gaussian-elimination/. The authors are listed in prov/FPBFNE.jsonl.
"""Structured Gaussian elimination (LaMacchia and Odlyzko) for a sparse system A x = c over Q or Z/pZ.
Light columns are pivoted away with rules that create fill-in only in a growing set of heavy
columns; the small remaining system on the heavy columns is solved densely, then the pivots
are back-substituted.
"""
from fractions import Fraction
from typing import Optional, Union
Number = Union[int, Fraction]
Row = dict[int, Number]
class Inconsistent(ValueError):
"""Raised when A x = c has no solution."""
def _field(p: Optional[int]):
if p is None:
return Fraction, lambda a: Fraction(1) / a
return (lambda a: a % p), (lambda a: pow(a, -1, p))
def dense_solve(rows: list[Row], rhs: list[Number], cols: list[int], p: Optional[int]) -> dict[int, Number]:
"""Gauss-Jordan elimination on the given rows restricted to cols; free variables are set to 0."""
norm, inv = _field(p)
m = [[norm(r.get(j, 0)) for j in cols] + [norm(v)] for r, v in zip(rows, rhs)]
k = len(cols)
where: list[int] = []
top = 0
for col in range(k):
piv = next((i for i in range(top, len(m)) if m[i][col] != 0), None)
if piv is None:
continue
m[top], m[piv] = m[piv], m[top]
s = inv(m[top][col])
m[top] = [norm(v * s) for v in m[top]]
for i in range(len(m)):
if i != top and m[i][col] != 0:
f = m[i][col]
m[i] = [norm(a - f * b) for a, b in zip(m[i], m[top])]
where.append(col)
top += 1
if any(row[k] != 0 for row in m[top:]):
raise Inconsistent("the dense part is inconsistent")
x = {j: norm(0) for j in cols}
for i, col in enumerate(where):
x[cols[col]] = m[i][k]
return x
def structured_solve(a: list[Row], c: list[Number], n: int, p: Optional[int] = None,
stats: Optional[dict] = None) -> list[Number]:
"""Return one solution x of A x = c, where A has n columns and row i is the sparse dict a[i].
Arithmetic is exact over Q (p is None) or over Z/pZ for a prime p. Raises Inconsistent when
there is no solution. If stats is a dict, it receives the number of heavy columns and pivots.
"""
if len(a) != len(c):
raise ValueError("a and c must have the same length")
norm, inv = _field(p)
rows = [{j: norm(v) for j, v in r.items() if norm(v) != 0} for r in a]
rhs = [norm(v) for v in c]
if any(not 0 <= j < n for r in rows for j in r):
raise ValueError("column index out of range")
active = set(range(len(rows))) # @line L1
light = set(range(n))
heavy: list[int] = []
pivots: list[tuple[int, int]] = []
x: list[Number] = [norm(0)] * n
colrows: dict[int, set[int]] = {j: set() for j in range(n)}
for i, r in enumerate(rows):
for j in r:
colrows[j].add(i)
def light_entries(i: int) -> list[int]:
return [j for j in rows[i] if j in light]
while light: # @line L2
empty = next((j for j in sorted(light) if not colrows[j]), None)
if empty is not None: # @line L3
light.discard(empty) # a free variable, x_j = 0
continue
single = next((j for j in sorted(light) if len(colrows[j]) == 1), None)
if single is not None: # @line L4
r = next(iter(colrows[single]))
pivots.append((r, single))
active.discard(r)
light.discard(single)
for j in rows[r]:
colrows[j].discard(r)
continue
cand = next((i for i in sorted(active) if len(light_entries(i)) == 1), None)
if cand is not None: # @line L5
r, j = cand, light_entries(cand)[0]
pivots.append((r, j))
active.discard(r)
light.discard(j)
for k in rows[r]:
colrows[k].discard(r)
s = inv(rows[r][j])
for i in sorted(colrows[j]): # the other active rows meeting column j
f = norm(rows[i][j] * s)
for k, v in rows[r].items(): # row r lives in column j and the heavy columns only
nv = norm(rows[i].get(k, 0) - f * v)
if nv == 0:
if k in rows[i]:
del rows[i][k]
colrows[k].discard(i)
else:
if k not in rows[i]:
colrows[k].add(i)
rows[i][k] = nv
rhs[i] = norm(rhs[i] - f * rhs[r])
continue
j = max(sorted(light), key=lambda col: len(colrows[col])) # @line L6
light.discard(j)
heavy.append(j)
act = sorted(active) # @line L7
xh = dense_solve([rows[i] for i in act], [rhs[i] for i in act], sorted(heavy), p)
for j, v in xh.items():
x[j] = v
for r, j in reversed(pivots): # @line L8
s = rhs[r] - sum(v * x[k] for k, v in rows[r].items() if k != j)
x[j] = norm(s * inv(rows[r][j]))
if stats is not None:
stats["heavy"] = len(heavy)
stats["pivots"] = len(pivots)
stats["dense_rows"] = len(act)
return x # @line L9
test_structured_gaussian_elimination.py 163 lines
# SPDX-License-Identifier: MIT
# Part of https://algorithm.tamnd.com/structured-gaussian-elimination/. The authors are listed in prov/FPBFNE.jsonl.
import math
import random
from fractions import Fraction
import pytest
from structured_gaussian_elimination import Inconsistent, structured_solve
def rank(rows, n, p):
m = [[Fraction(r.get(j, 0)) if p is None else r.get(j, 0) % p for j in range(n)] for r in rows]
rk = 0
for col in range(n):
piv = next((i for i in range(rk, len(m)) if m[i][col] != 0), None)
if piv is None:
continue
m[rk], m[piv] = m[piv], m[rk]
for i in range(rk + 1, len(m)):
if m[i][col] != 0:
f = m[i][col] / m[rk][col] if p is None else m[i][col] * pow(m[rk][col], -1, p) % p
m[i] = [a - f * b if p is None else (a - f * b) % p for a, b in zip(m[i], m[rk])]
rk += 1
return rk
def consistent(rows, c, n, p):
aug = [{**r, n: v} for r, v in zip(rows, c)]
return rank(rows, n, p) == rank(aug, n + 1, p)
def residual_ok(rows, c, x, p):
for r, v in zip(rows, c):
s = sum(a * x[j] for j, a in r.items()) - v
if (s if p is None else s % p) != 0:
return False
return True
def random_sparse(rng, m, n, p, density):
rows = []
for _ in range(m):
r = {}
for j in range(n):
if rng.random() < density:
v = rng.randrange(1, p) if p else rng.choice([-3, -2, -1, 1, 2, 3])
r[j] = v
rows.append(r)
return rows
def test_worked_example():
# x0 + x1 = 3, x1 + x2 + x4 = 10, x2 + x3 = 7, x0 + x3 + x4 = 10, x2 + x4 = 8; the solution is (1, 2, 3, 4, 5)
a = [{0: 1, 1: 1}, {1: 1, 2: 1, 4: 1}, {2: 1, 3: 1}, {0: 1, 3: 1, 4: 1}, {2: 1, 4: 1}]
c = [3, 10, 7, 10, 8]
stats = {}
x = structured_solve(a, c, 5, stats=stats)
assert x == [1, 2, 3, 4, 5]
assert stats == {"heavy": 1, "pivots": 4, "dense_rows": 1}
assert structured_solve(a, c, 5, p=7) == [1, 2, 3, 4, 5]
def test_edges():
assert structured_solve([], [], 0) == []
assert structured_solve([], [], 3) == [0, 0, 0] # no equations, every x works
assert structured_solve([{}], [0], 0) == []
with pytest.raises(Inconsistent):
structured_solve([{}], [1], 2) # a zero row with a nonzero right side
assert structured_solve([{0: 5}], [10], 1) == [2]
assert structured_solve([{0: 3}], [1], 1, p=7) == [5]
x = structured_solve([{1: 2}, {1: 4}], [2, 4], 3) # zero columns 0 and 2
assert x == [0, 1, 0]
with pytest.raises(Inconsistent):
structured_solve([{0: 1, 1: 1}, {0: 2, 1: 2}], [1, 3], 2) # singular and inconsistent
with pytest.raises(ValueError):
structured_solve([{3: 1}], [1], 2)
def test_triangular_needs_no_heavy_columns():
rng = random.Random(3)
for _ in range(50):
n = rng.randrange(1, 30)
perm_r, perm_c = list(range(n)), list(range(n))
rng.shuffle(perm_r)
rng.shuffle(perm_c)
rows = [dict() for _ in range(n)]
for i in range(n):
rows[perm_r[i]][perm_c[i]] = rng.randrange(1, 101)
for j in range(i + 1, n):
if rng.random() < 0.2:
rows[perm_r[i]][perm_c[j]] = rng.randrange(101)
x0 = [rng.randrange(101) for _ in range(n)]
c = [sum(v * x0[j] for j, v in r.items()) % 101 for r in rows]
stats = {}
assert structured_solve(rows, c, n, p=101, stats=stats) == x0
assert stats["heavy"] == 0
@pytest.mark.parametrize("p", [2, 3, 7, 101])
def test_random_mod_p(p):
rng = random.Random(1000 + p)
for _ in range(250):
m, n = rng.randrange(0, 13), rng.randrange(0, 13)
rows = random_sparse(rng, m, n, p, rng.choice([0.1, 0.25, 0.5]))
if rng.random() < 0.6:
x0 = [rng.randrange(p) for _ in range(n)]
c = [sum(v * x0[j] for j, v in r.items()) % p for r in rows]
else:
c = [rng.randrange(p) for _ in range(m)]
if consistent(rows, c, n, p):
x = structured_solve(rows, c, n, p)
assert residual_ok(rows, c, x, p) and all(0 <= v < p for v in x)
else:
with pytest.raises(Inconsistent):
structured_solve(rows, c, n, p)
def test_random_over_q_and_full_rank_uniqueness():
rng = random.Random(77)
for _ in range(200):
m, n = rng.randrange(0, 11), rng.randrange(0, 11)
rows = random_sparse(rng, m, n, None, rng.choice([0.15, 0.3, 0.6]))
x0 = [Fraction(rng.randrange(-9, 10), rng.randrange(1, 5)) for _ in range(n)]
c = [sum(v * x0[j] for j, v in r.items()) for r in rows]
x = structured_solve(rows, c, n)
assert residual_ok(rows, c, x, None)
if rank(rows, n, None) == n:
assert x == x0
def test_large_entries():
rng = random.Random(9)
for _ in range(30):
n = rng.randrange(1, 9)
rows = [{j: rng.randrange(-2**64, 2**64) for j in range(n) if rng.random() < 0.5} for _ in range(n + 2)]
x0 = [rng.randrange(-2**40, 2**40) for _ in range(n)]
c = [sum(v * x0[j] for j, v in r.items()) for r in rows]
x = structured_solve(rows, c, n)
assert residual_ok(rows, c, x, None)
if rank(rows, n, None) == n:
assert x == x0
def test_sparse_reduction_shrinks_the_system():
# 310 rows with 10 entries in 300 columns, column j chosen with probability about 1/(j+1),
# like the factor base of index calculus: the dense part has far fewer columns
rng = random.Random(5)
p, n, m, w = 1_000_003, 300, 310, 10
rows = []
for _ in range(m):
r = {}
while len(r) < w:
j = int(math.exp(rng.random() * math.log(n + 1))) - 1
r[min(j, n - 1)] = rng.randrange(1, p)
rows.append(r)
x0 = [rng.randrange(p) for _ in range(n)]
c = [sum(v * x0[j] for j, v in r.items()) % p for r in rows]
stats = {}
x = structured_solve(rows, c, n, p, stats)
assert residual_ok(rows, c, x, p)
assert stats["heavy"] < n // 5
What the next level needs
- the tests did not run in this build