Berkowitz algorithm in Python
Division-free algorithm for the characteristic polynomial, valid over any commutative ring.
- File
berkowitz_algorithm.py, 49 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/berkowitz-algorithm/. The authors are listed in prov/PM8CKR.jsonl.
from typing import Any, Sequence
def berkowitz(a: Sequence[Sequence[Any]], one: Any = 1) -> list[Any]:
"""Return the coefficients of det(xI - a), highest degree first, for a square matrix a.
Uses only addition, subtraction and multiplication, so it works over every commutative ring
(ints, residues modulo a composite number, polynomials, ...). `one` is the ring's unit.
Cost: about n^4/2 ring operations.
"""
n = len(a)
if n == 0: # @line L1
return [one]
p = [one, -a[n - 1][n - 1]] # @line L2
for r in range(n - 2, -1, -1): # @line L3
s = n - r - 1 # @line L4
q = [one, -a[r][r]] # @line L5
v = list(a[r][r + 1:])
for k in range(s): # @line L6
g = v[0] * a[r + 1][r]
for t in range(1, s):
g = g + v[t] * a[r + 1 + t][r]
q.append(-g) # @line L7
if k < s - 1: # @line L8
w = []
for c in range(s):
acc = v[0] * a[r + 1][r + 1 + c]
for t in range(1, s):
acc = acc + v[t] * a[r + 1 + t][r + 1 + c]
w.append(acc)
v = w
new = []
for i in range(s + 2): # @line L9
acc = None
for k in range(max(0, i - s - 1), min(i, s) + 1):
term = q[i - k] * p[k]
acc = term if acc is None else acc + term
new.append(acc)
p = new
return p # @line L10
def berkowitz_determinant(a: Sequence[Sequence[Any]], one: Any = 1) -> Any:
"""Return det a over a commutative ring, as (-1)^n times the constant coefficient."""
c0 = berkowitz(a, one)[-1]
return c0 if len(a) % 2 == 0 else -c0
test_berkowitz_algorithm.py 108 lines
# SPDX-License-Identifier: MIT
# Part of https://algorithm.tamnd.com/berkowitz-algorithm/. The authors are listed in prov/PM8CKR.jsonl.
import random
from itertools import combinations, permutations
from berkowitz_algorithm import berkowitz, berkowitz_determinant
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 charpoly_minors(a):
"""Coefficients of det(xI - a), highest first, from sums of principal minors."""
n = len(a)
out = [1]
for k in range(1, n + 1):
e = sum(det_leibniz([[a[i][j] for j in s] for i in s]) for s in combinations(range(n), k))
out.append((-1) ** k * e)
return out
class ModN:
"""An element of Z/NZ for any N >= 2, a ring that may have zero divisors."""
def __init__(self, v: int, n: int) -> None:
self.v, self.n = v % n, n
def __add__(self, o):
return ModN(self.v + o.v, self.n)
def __sub__(self, o):
return ModN(self.v - o.v, self.n)
def __mul__(self, o):
return ModN(self.v * o.v, self.n)
def __neg__(self):
return ModN(-self.v, self.n)
def matmul(x, y):
return [[sum(x[i][k] * y[k][j] for k in range(len(y))) for j in range(len(y[0]))] for i in range(len(x))]
def test_examples() -> None:
assert berkowitz([[2, 1, 0], [1, 3, 1], [0, 1, 4]]) == [1, -9, 24, -18]
assert berkowitz([[1, 2], [3, 4]]) == [1, -5, -2]
assert berkowitz([[0, 0, -5], [1, 0, 2], [0, 1, 0]]) == [1, 0, -2, 5]
assert berkowitz_determinant([[0, 1, 2, 3], [1, 0, 1, 2], [2, 1, 0, 1], [3, 2, 1, 0]]) == -12
def test_edges() -> None:
assert berkowitz([]) == [1]
assert berkowitz_determinant([]) == 1
assert berkowitz([[7]]) == [1, -7]
assert berkowitz_determinant([[7]]) == 7
assert berkowitz([[0, 0], [0, 0]]) == [1, 0, 0]
assert berkowitz([[0, 1], [0, 0]]) == [1, 0, 0]
def test_random_integer() -> None:
rng = random.Random(1984)
for _ in range(600):
n = rng.randint(1, 5)
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.2:
for i in range(n):
a[i][0] = 0
c = berkowitz(a)
assert c == charpoly_minors(a)
assert berkowitz_determinant(a) == det_leibniz(a)
# Cayley-Hamilton: c(A) = 0
acc = [[0] * n for _ in range(n)]
for coef in c:
acc = matmul(acc, a)
for i in range(n):
acc[i][i] += coef
assert acc == [[0] * n for _ in range(n)]
def test_ring_with_zero_divisors() -> None:
rng = random.Random(6)
for modulus in (4, 6, 12, 2**64):
for _ in range(150):
n = rng.randint(1, 5)
a = [[rng.randrange(modulus) for _ in range(n)] for _ in range(n)]
got = berkowitz([[ModN(x, modulus) for x in row] for row in a], ModN(1, modulus))
assert [x.v for x in got] == [x % modulus for x in charpoly_minors(a)]
def test_large_entries() -> None:
rng = random.Random(9)
for _ in range(15):
n = 5
a = [[rng.randint(-10**25, 10**25) for _ in range(n)] for _ in range(n)]
assert berkowitz(a) == charpoly_minors(a)
What the next level needs
- the tests did not run in this build