Faddeev-LeVerrier algorithm in Python
Computes the characteristic polynomial through a recurrence on matrix powers, dividing by integers.
- File
faddeev_leverrier_algorithm.py, 62 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/faddeev-leverrier-algorithm/. The authors are listed in prov/6P3FW5.jsonl.
from fractions import Fraction
from typing import Any, Sequence
def _divide(x: Any, k: int) -> Any:
"""Divide x by the integer k. For an int the quotient must be exact."""
if isinstance(x, int):
q, r = divmod(x, k)
assert r == 0, "inexact division by k"
return q
return x / k
def faddeev_leverrier(a: Sequence[Sequence[Any]]) -> tuple[list[Any], list[list[Any]]]:
"""Return (c, m): c lists det(xI - a) highest degree first, and m is the matrix M_n.
Works over a field of characteristic 0 or greater than n, and over the integers, where
every division by k is exact. adj(a) = (-1)^(n-1) M_n, and a M_n = -c_0 I.
Cost: n matrix products, 2n^4 + O(n^3) field operations.
"""
n = len(a)
c = [1] + [0] * n # @line L1
nm = [[0] * n for _ in range(n)]
m = nm
for k in range(1, n + 1): # @line L2
m = [[nm[i][j] + (c[k - 1] if i == j else 0) for j in range(n)] for i in range(n)] # @line L3
nm = [[sum((a[i][t] * m[t][j] for t in range(1, n)), a[i][0] * m[0][j]) for j in range(n)] for i in range(n)] # @line L4
c[k] = -_divide(sum((nm[i][i] for i in range(1, n)), nm[0][0]), k) # @line L5
return c, m # @line L6
def characteristic_polynomial(a: Sequence[Sequence[Any]]) -> list[Any]:
"""Coefficients of det(xI - a), highest degree first."""
return faddeev_leverrier(a)[0]
def determinant(a: Sequence[Sequence[Any]]) -> Any:
"""det a = (-1)^n c_0."""
c0 = faddeev_leverrier(a)[0][-1]
return c0 if len(a) % 2 == 0 else -c0
def adjugate(a: Sequence[Sequence[Any]]) -> list[list[Any]]:
"""adj a = (-1)^(n-1) M_n, the matrix of cofactors transposed."""
n = len(a)
if n == 0:
return []
m = faddeev_leverrier(a)[1]
return m if n % 2 == 1 else [[-x for x in row] for row in m]
def inverse(a: Sequence[Sequence[Any]]) -> list[list[Fraction]] | None:
"""Return the inverse of an integer or rational matrix as Fractions, or None if singular."""
n = len(a)
q = [[Fraction(x) for x in row] for row in a]
c, m = faddeev_leverrier(q)
if c[-1] == 0:
return None
return [[-x / c[-1] for x in row] for row in m] if n else []
test_faddeev_leverrier_algorithm.py 140 lines
# SPDX-License-Identifier: MIT
# Part of https://algorithm.tamnd.com/faddeev-leverrier-algorithm/. The authors are listed in prov/6P3FW5.jsonl.
import random
from fractions import Fraction
from itertools import combinations, permutations
import pytest
from faddeev_leverrier_algorithm import adjugate, characteristic_polynomial, determinant, faddeev_leverrier, inverse
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):
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
def adj_cofactors(a):
n = len(a)
if n == 1:
return [[1]]
return [[(-1) ** (i + k) * det_leibniz([[a[r][c] for c in range(n) if c != i] for r in range(n) if r != k])
for k in range(n)] for i in range(n)]
class ModP:
P = 101
def __init__(self, v: int) -> None:
self.v = v % ModP.P
def __add__(self, o):
return ModP(self.v + (o.v if isinstance(o, ModP) else o))
__radd__ = __add__
def __mul__(self, o):
return ModP(self.v * (o.v if isinstance(o, ModP) else o))
__rmul__ = __mul__
def __neg__(self):
return ModP(-self.v)
def __truediv__(self, k: int):
if k % ModP.P == 0:
raise ZeroDivisionError("k is not invertible modulo P")
return ModP(self.v * pow(k, -1, ModP.P))
def __eq__(self, o):
return self.v == (o.v if isinstance(o, ModP) else o % ModP.P)
def test_example() -> None:
m = [[2, 1, 0], [1, 3, 1], [0, 1, 4]]
c, mn = faddeev_leverrier(m)
assert c == [1, -9, 24, -18]
assert mn == [[11, -4, 1], [-4, 8, -2], [1, -2, 5]]
assert inverse(m) == [[Fraction(x, 18) for x in row] for row in mn]
def test_edges() -> None:
assert faddeev_leverrier([]) == ([1], [])
assert determinant([]) == 1 and inverse([]) == []
assert characteristic_polynomial([[5]]) == [1, -5]
assert adjugate([[5]]) == [[1]]
assert inverse([[0]]) is None
assert characteristic_polynomial([[0, 1], [0, 0]]) == [1, 0, 0]
assert inverse([[1, 2], [2, 4]]) is None
def test_random_integer() -> None:
rng = random.Random(1840)
for _ in range(600):
n = rng.randint(1, 5)
a = [[rng.randint(-5, 5) for _ in range(n)] for _ in range(n)]
if rng.random() < 0.2:
a[rng.randrange(n)] = [0] * n
if rng.random() < 0.15:
for i in range(n):
a[i][rng.randrange(n)] = 0
assert characteristic_polynomial(a) == charpoly_minors(a)
assert determinant(a) == det_leibniz(a)
assert adjugate(a) == adj_cofactors(a)
inv = inverse(a)
if det_leibniz(a) == 0:
assert inv is None
else:
prod = [[sum(a[i][t] * inv[t][j] for t in range(n)) for j in range(n)] for i in range(n)]
assert prod == [[int(i == j) for j in range(n)] for i in range(n)]
def test_rational_entries() -> None:
rng = random.Random(3)
for _ in range(100):
n = rng.randint(1, 4)
a = [[Fraction(rng.randint(-6, 6), rng.randint(1, 5)) for _ in range(n)] for _ in range(n)]
assert characteristic_polynomial(a) == charpoly_minors(a)
def test_mod_p_large_characteristic() -> None:
rng = random.Random(4)
for _ in range(150):
n = rng.randint(1, 5)
a = [[rng.randrange(101) for _ in range(n)] for _ in range(n)]
got = characteristic_polynomial([[ModP(x) for x in row] for row in a])
assert all(g == e for g, e in zip(got, charpoly_minors(a)))
def test_small_characteristic_fails() -> None:
ModP.P = 2
try:
with pytest.raises(ZeroDivisionError):
characteristic_polynomial([[ModP(1), ModP(1)], [ModP(1), ModP(0)]])
finally:
ModP.P = 101
def test_large_entries() -> None:
rng = random.Random(5)
for _ in range(10):
n = 5
a = [[rng.randint(-10**20, 10**20) for _ in range(n)] for _ in range(n)]
assert characteristic_polynomial(a) == charpoly_minors(a)
What the next level needs
- the tests did not run in this build