Sylvester matrix in Python
Matrix of coefficients of two polynomials whose determinant is their resultant.
- File
sylvester_matrix.py, 129 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/sylvester-matrix/. The authors are listed in prov/NEATH6.jsonl.
"""The Sylvester matrix of two univariate polynomials, stored implicitly.
A polynomial is a list of coefficients from the highest degree down:
f = [a_m, ..., a_0] has formal degree m = len(f) - 1. The leading coefficient
may be zero (a formal degree). Coefficients are elements of an integral
domain with exact division, such as int or Fraction.
Syl(f, g) is the (m + n) x (m + n) matrix whose row i < n holds the
coefficients of x^(n-1-i) f and whose row n + j holds those of x^(m-1-j) g,
from x^(m+n-1) down to x^0.
Cost (ring operations, which over Z means operations on integers whose size
grows with the input, not bit operations): entry O(1), combine O(m n),
determinant O((m + n)^3) by Bareiss elimination with exact divisions, rank
O((m + n)^3). Space: O(m + n) for the matrix itself.
"""
from __future__ import annotations
from typing import Any, Sequence
class SylvesterMatrix:
"""Implicit Sylvester matrix: only the two coefficient vectors are stored."""
def __init__(self, f: Sequence[Any], g: Sequence[Any]) -> None:
if len(f) == 0 or len(g) == 0:
raise ValueError("each polynomial needs a formal degree, at least one coefficient")
self.f = list(f)
self.g = list(g)
self.m = len(self.f) - 1
self.n = len(self.g) - 1
self.size = self.m + self.n
def check_invariant(self) -> None:
"""The dense matrix agrees with entry(), and each row is a shifted coefficient vector."""
m, n, N = self.m, self.n, self.size
rows = self.to_dense()
assert len(rows) == N and all(len(r) == N for r in rows)
for i in range(n):
assert rows[i][i:i + m + 1] == self.f
for j in range(m):
assert rows[n + j][j:j + n + 1] == self.g
def entry(self, r: int, c: int) -> Any:
"""Entry in row r, column c. Cost O(1)."""
if not (0 <= r < self.size and 0 <= c < self.size): # @line E1
raise IndexError("index out of range")
if r < self.n: # @line E2
k = c - r # @line E3
return self.f[k] if 0 <= k <= self.m else 0 * self.f[0]
k = c - (r - self.n) # @line E4
return self.g[k] if 0 <= k <= self.n else 0 * self.g[0]
def to_dense(self) -> list[list[Any]]:
"""The full matrix as a list of rows. Cost O((m + n)^2)."""
return [[self.entry(r, c) for c in range(self.size)] for r in range(self.size)]
def combine(self, s: Sequence[Any], t: Sequence[Any]) -> list[Any]:
"""Row vector (s, t) times Syl(f, g): the coefficients of s f + t g.
s has n coefficients (degree < n) and t has m (degree < m), highest first.
The result has m + n coefficients, highest first. Cost O(m n).
"""
if len(s) != self.n or len(t) != self.m: # @line C1
raise ValueError("s needs n coefficients and t needs m")
out = [0 * self.f[0]] * self.size # @line C2
for i, si in enumerate(s): # @line C3
for k, a in enumerate(self.f):
out[i + k] = out[i + k] + si * a
for j, tj in enumerate(t): # @line C4
for k, b in enumerate(self.g):
out[j + k] = out[j + k] + tj * b
return out # @line C5
def determinant(self) -> Any:
"""det Syl(f, g), the resultant Res_(m,n)(f, g), by Bareiss elimination.
Every division is exact in the coefficient domain. Cost O((m + n)^3) ring operations.
"""
N = self.size
if N == 0: # @line L1
return 1
A = self.to_dense() # @line L2
sign, prev = 1, 1 # @line L3
for k in range(N - 1): # @line L4
p = next((r for r in range(k, N) if A[r][k] != 0), None) # @line L5
if p is None: # @line L6
return 0 * A[0][0]
if p != k: # @line L7
A[k], A[p] = A[p], A[k]
sign = -sign
for i in range(k + 1, N): # @line L8
for j in range(k + 1, N):
A[i][j] = _exact_div(A[i][j] * A[k][k] - A[i][k] * A[k][j], prev) # @line L9
A[i][k] = 0 * A[i][k]
prev = A[k][k] # @line L10
return sign * A[N - 1][N - 1] # @line L11
def rank(self) -> int:
"""Rank over the fraction field, by fraction-free elimination. Cost O((m + n)^3)."""
A = self.to_dense()
N = self.size
rank, prev, row = 0, 1, 0
for col in range(N):
p = next((r for r in range(row, N) if A[r][col] != 0), None)
if p is None:
continue
A[row], A[p] = A[p], A[row]
for i in range(row + 1, N):
for j in range(col + 1, N):
A[i][j] = _exact_div(A[i][j] * A[row][col] - A[i][col] * A[row][j], prev)
A[i][col] = 0 * A[i][col]
prev = A[row][col]
row += 1
rank += 1
return rank
def _exact_div(a: Any, b: Any) -> Any:
"""a / b where b divides a exactly (int or Fraction)."""
if isinstance(a, int) and isinstance(b, int):
q, r = divmod(a, b)
if r != 0:
raise ArithmeticError("division is not exact")
return q
return a / b
test_sylvester_matrix.py 157 lines
# SPDX-License-Identifier: MIT
# Part of https://algorithm.tamnd.com/sylvester-matrix/. The authors are listed in prov/NEATH6.jsonl.
import itertools
import random
from fractions import Fraction
import pytest
from sylvester_matrix import SylvesterMatrix
def leibniz_det(A):
"""Determinant by the permutation expansion (independent oracle)."""
N = len(A)
total = 0
for perm in itertools.permutations(range(N)):
inv = sum(1 for a in range(N) for b in range(a + 1, N) if perm[a] > perm[b])
prod = 1
for r in range(N):
prod *= A[r][perm[r]]
total += -prod if inv % 2 else prod
return total
def gauss_rank(A):
"""Rank over Q by Gaussian elimination with Fractions."""
M = [[Fraction(x) for x in row] for row in A]
rank, row = 0, 0
N = len(M)
for col in range(N):
p = next((r for r in range(row, N) if M[r][col] != 0), None)
if p is None:
continue
M[row], M[p] = M[p], M[row]
for i in range(row + 1, N):
q = M[i][col] / M[row][col]
for j in range(col, N):
M[i][j] -= q * M[row][j]
row += 1
rank += 1
return rank
def from_roots(lead, roots):
"""Coefficients, highest first, of lead * prod (x - r)."""
c = [lead]
for r in roots:
c = [a - r * b for a, b in zip(c + [0], [0] + c)]
return c
def poly_gcd_degree(f, g):
"""Degree of gcd(f, g) over Q by the Euclidean algorithm."""
def strip(p):
p = [Fraction(x) for x in p]
while p and p[0] == 0:
p.pop(0)
return p
a, b = strip(f), strip(g)
while b:
while len(a) >= len(b) and a:
q = a[0] / b[0]
a = strip([x - q * y for x, y in zip(a, b + [0] * (len(a) - len(b)))])
a, b = b, a
return len(a) - 1
def naive_combine(f, g, s, t):
"""s f + t g by direct multiplication, highest first, padded to m + n."""
m, n = len(f) - 1, len(g) - 1
out = [0] * (m + n)
for i, si in enumerate(s):
for k, a in enumerate(f):
out[i + k] += si * a
for j, tj in enumerate(t):
for k, b in enumerate(g):
out[j + k] += tj * b
return out
def test_example():
f = [1, -3, 2] # x^2 - 3x + 2 = (x - 1)(x - 2)
g = [1, 0, 1] # x^2 + 1
S = SylvesterMatrix(f, g)
assert S.to_dense() == [[1, -3, 2, 0], [0, 1, -3, 2], [1, 0, 1, 0], [0, 1, 0, 1]]
assert S.determinant() == 10
assert S.rank() == 4
S.check_invariant()
h = SylvesterMatrix([1, -3, 2], [1, 1, -2]) # common root x = 1
assert h.determinant() == 0 and h.rank() == 3
def test_edges():
assert SylvesterMatrix([5], [7]).determinant() == 1 # empty matrix
assert SylvesterMatrix([3], [1, 2, 3]).determinant() == 9 # a0^n
assert SylvesterMatrix([1, 2, 3], [3]).determinant() == 9 # b0^m
assert SylvesterMatrix([0, 1, 1], [0, 1, 2]).determinant() == 0 # both formal leads zero
assert SylvesterMatrix([0, 0], [1, 0]).determinant() == 0 # f = 0 with formal degree 1
with pytest.raises(ValueError):
SylvesterMatrix([], [1])
with pytest.raises(IndexError):
SylvesterMatrix([1, 1], [1, 1]).entry(2, 0)
def test_random_against_leibniz_and_combine():
rng = random.Random(59)
for _ in range(500):
m, n = rng.randint(0, 3), rng.randint(0, 3)
f = [rng.randint(-4, 4) for _ in range(m + 1)]
g = [rng.randint(-4, 4) for _ in range(n + 1)]
S = SylvesterMatrix(f, g)
S.check_invariant()
A = S.to_dense()
assert S.determinant() == leibniz_det(A)
assert S.rank() == gauss_rank(A)
if f[0] != 0 or g[0] != 0:
# rank = m + n - deg gcd once one formal leading coefficient is nonzero
assert S.rank() == m + n - poly_gcd_degree(f, g)
s = [rng.randint(-3, 3) for _ in range(n)]
t = [rng.randint(-3, 3) for _ in range(m)]
assert S.combine(s, t) == naive_combine(f, g, s, t)
# Res(g, f) = (-1)^(m n) Res(f, g)
assert SylvesterMatrix(g, f).determinant() == (-1) ** (m * n) * S.determinant()
def test_product_formula_and_gcd_rank():
rng = random.Random(7)
for _ in range(300):
m, n = rng.randint(1, 4), rng.randint(1, 4)
ra = [rng.randint(-3, 3) for _ in range(m)]
rb = [rng.randint(-3, 3) for _ in range(n)]
a, b = rng.choice([1, -1, 2, 3]), rng.choice([1, -2, 5])
f, g = from_roots(a, ra), from_roots(b, rb)
S = SylvesterMatrix(f, g)
expected = a ** n * b ** m
for x in ra:
for y in rb:
expected *= x - y
assert S.determinant() == expected
assert (S.determinant() == 0) == (poly_gcd_degree(f, g) >= 1)
assert S.rank() == m + n - poly_gcd_degree(f, g)
def test_fractions_and_large_sizes():
rng = random.Random(11)
for _ in range(40):
m, n = rng.randint(3, 7), rng.randint(3, 7)
f = [Fraction(rng.randint(-9, 9), rng.randint(1, 5)) for _ in range(m + 1)]
g = [Fraction(rng.randint(-9, 9), rng.randint(1, 5)) for _ in range(n + 1)]
f[0] = f[0] or Fraction(1)
g[0] = g[0] or Fraction(1)
S = SylvesterMatrix(f, g)
d = S.determinant()
# Scaling: Res(c f, g) = c^n Res(f, g)
assert SylvesterMatrix([3 * x for x in f], g).determinant() == 3 ** n * d
assert (d == 0) == (poly_gcd_degree(f, g) >= 1)
assert S.rank() == m + n - poly_gcd_degree(f, g)
What the next level needs
- the tests did not run in this build