Bareiss algorithm in Python
Fraction-free elimination over the integers in which every division is exact, keeping intermediate sizes polynomially bounded.
- File
bareiss_algorithm.py, 76 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/bareiss-algorithm/. The authors are listed in prov/S667TN.jsonl.
from typing import Sequence
def bareiss(a: Sequence[Sequence[int]]) -> tuple[int, list[list[int]], list[int]]:
"""Return (det, u, perm) for the n by n integer matrix a, by fraction-free elimination.
Every division is exact. u is the final working matrix: for k <= j, u[k][j] is the minor of
the row-permuted matrix in rows 0..k and columns 0..k-1, j. perm lists the original row index
of each row of u. When a is singular, det is 0 and u, perm describe the state at that point.
Cost: O(n^3) integer operations on numbers of O(n (b + log n)) bits.
"""
n = len(a)
u = [list(row) for row in a]
perm = list(range(n))
if n == 0: # @line L1
return 1, u, perm
s, prev = 1, 1 # @line L2
for k in range(n - 1): # @line L3
if u[k][k] == 0: # @line L4
p = next((i for i in range(k + 1, n) if u[i][k] != 0), None) # @line L5
if p is None: # @line L6
return 0, u, perm
u[k], u[p] = u[p], u[k] # @line L7
perm[k], perm[p] = perm[p], perm[k]
s = -s
pivot, row_k = u[k][k], u[k]
for i in range(k + 1, n): # @line L8
row_i = u[i]
f = row_i[k]
for j in range(k + 1, n): # @line L9
row_i[j] = (pivot * row_i[j] - f * row_k[j]) // prev # @line L10
row_i[k] = 0
prev = pivot # @line L11
return s * u[n - 1][n - 1], u, perm # @line L12
def bareiss_determinant(a: Sequence[Sequence[int]]) -> int:
"""Return the determinant of the square integer matrix a. Cost: O(n^3) integer operations."""
return bareiss(a)[0]
def fraction_free_echelon(a: Sequence[Sequence[int]], n: int | None = None) -> tuple[list[list[int]], int, list[int]]:
"""Return (e, r, pivots): an echelon form e of the m by n integer matrix a with integer
entries, obtained by fraction-free row operations, its rank r and its pivot columns.
The nonzero rows of e span the row space of a over Q. Every division is exact.
Cost: O(m n min(m, n)) integer operations.
"""
e = [list(row) for row in a]
m = len(e)
if n is None:
n = len(e[0]) if m else 0
r, prev, pivots = 0, 1, []
for j in range(n):
if r == m:
break
p = next((i for i in range(r, m) if e[i][j] != 0), None)
if p is None:
continue
e[r], e[p] = e[p], e[r]
pivot, row_r = e[r][j], e[r]
for i in range(r + 1, m):
row_i = e[i]
f = row_i[j]
for t in range(j + 1, n):
q, rem = divmod(pivot * row_i[t] - f * row_r[t], prev)
assert rem == 0, "inexact division"
row_i[t] = q
row_i[j] = 0
prev = pivot
pivots.append(j)
r += 1
return e, r, pivots
test_bareiss_algorithm.py 106 lines
# SPDX-License-Identifier: MIT
# Part of https://algorithm.tamnd.com/bareiss-algorithm/. The authors are listed in prov/S667TN.jsonl.
import random
from itertools import combinations, permutations
from bareiss_algorithm import bareiss, bareiss_determinant, fraction_free_echelon
def det_leibniz(a):
n = len(a)
total = 0
for p in permutations(range(n)):
inv = sum(1 for i in range(n) for j in range(i + 1, n) if p[i] > p[j])
t = 1
for i in range(n):
t *= a[p[i]][i]
total += -t if inv % 2 else t
return total
def rank_brute(a, n):
m = len(a)
for t in range(min(m, n), 0, -1):
for rows in combinations(range(m), t):
for cols in combinations(range(n), t):
if det_leibniz([[a[i][j] for j in cols] for i in rows]) != 0:
return t
return 0
def test_examples() -> None:
m = [[2, 1, 0], [1, 3, 1], [0, 1, 4]]
d, u, perm = bareiss(m)
assert d == 18 and u[1][1:] == [5, 2] and u[2][2] == 18 and perm == [0, 1, 2]
assert bareiss_determinant([[0, 1, 2, 3], [1, 0, 1, 2], [2, 1, 0, 1], [3, 2, 1, 0]]) == -12
def test_edges() -> None:
assert bareiss_determinant([]) == 1
assert bareiss_determinant([[7]]) == 7
assert bareiss_determinant([[0]]) == 0
assert bareiss_determinant([[0, 1], [1, 0]]) == -1
assert bareiss_determinant([[0, 0], [0, 5]]) == 0
assert bareiss_determinant([[1, 2], [2, 4]]) == 0
assert bareiss_determinant([[1, 2, 3], [0, 0, 0], [4, 5, 6]]) == 0
assert fraction_free_echelon([]) == ([], 0, [])
assert fraction_free_echelon([[0, 0, 3]]) == ([[0, 0, 3]], 1, [2])
def test_random_determinant() -> None:
rng = random.Random(1968)
for _ in range(800):
n = rng.randint(1, 6)
a = [[rng.randint(-4, 4) for _ in range(n)] for _ in range(n)]
if rng.random() < 0.2:
a[rng.randrange(n)] = [0] * n
if rng.random() < 0.3:
for i in range(n):
a[i][0] = 0 if i < n - 1 else a[i][0]
if n >= 2 and rng.random() < 0.15:
a[1] = [2 * x for x in a[0]]
assert bareiss_determinant(a) == det_leibniz(a)
def test_intermediate_entries_are_minors() -> None:
rng = random.Random(77)
for _ in range(400):
n = rng.randint(1, 5)
a = [[rng.randint(-5, 5) for _ in range(n)] for _ in range(n)]
d, u, perm = bareiss(a)
if d == 0:
continue
pa = [a[i] for i in perm]
for k in range(n):
for j in range(k, n):
cols = list(range(k)) + [j]
minor = det_leibniz([[pa[i][c] for c in cols] for i in range(k + 1)])
assert u[k][j] == minor
def test_random_echelon() -> None:
rng = random.Random(4096)
for _ in range(600):
m, n = rng.randint(0, 4), rng.randint(0, 5)
a = [[rng.randint(-3, 3) for _ in range(n)] for _ in range(m)]
if m >= 2 and rng.random() < 0.3:
a[-1] = [x - y for x, y in zip(a[0], a[1])]
e, r, piv = fraction_free_echelon(a, n)
assert r == rank_brute(a, n) == len(piv)
for i in range(m):
nz = [j for j in range(n) if e[i][j] != 0]
if i < r:
assert nz[0] == piv[i]
else:
assert not nz
assert piv == sorted(piv)
assert rank_brute(a + e, n) == r
def test_large_entries() -> None:
rng = random.Random(5)
for _ in range(20):
n = 5
a = [[rng.randint(-10**40, 10**40) for _ in range(n)] for _ in range(n)]
assert bareiss_determinant(a) == det_leibniz(a)
What the next level needs
- the tests did not run in this build