F5 algorithm in Python
Signature-based Gröbner algorithm that avoids reducing many useless syzygy pairs.
- File
f5_algorithm.py, 183 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/f5-algorithm/. The authors are listed in prov/2QB464.jsonl.
"""A signature-based Groebner basis algorithm in the style of F5, over Q.
A polynomial in K[x_1, ..., x_n] is a dict from exponent tuples to nonzero
Fractions. A monomial order is one of "lex", "grlex", "grevlex" with
x_1 > x_2 > ... > x_n. Signatures are terms u e_i of the free module R^s,
stored as (i, u), and compared position over term: first by the index i,
then by u in the same monomial order. This is the incremental order of F5.
f5(F, order) returns the reduced Groebner basis of <F> and statistics
{"reductions": ..., "zero_reductions": ...}. The criteria are the F5
(syzygy) criterion with the principal syzygies LM(g) e_i of the earlier
generators, the signatures of reductions to zero, and the rewrite rule that
keeps one candidate per signature and discards singular top-reducible ones.
Cost model: operations in Q. For a homogeneous regular sequence of n
polynomials in grevlex, F5 performs no reduction to zero (Faugere 2002), and
the matrix variant costs O(N^omega) operations, N the number of monomials of
degree at most the degree of regularity (Bardet, Faugere and Salvy 2015). In general the worst case is
doubly exponential in n (Mayr and Meyer 1982).
"""
from fractions import Fraction
from typing import Callable, Dict, List, Set, Tuple
Mono = Tuple[int, ...]
Poly = Dict[Mono, Fraction]
Key = Callable[[Mono], Tuple[object, ...]]
Sig = Tuple[int, Mono]
ORDERS: Dict[str, Key] = {
"lex": lambda a: tuple(a),
"grlex": lambda a: (sum(a),) + tuple(a),
"grevlex": lambda a: (sum(a),) + tuple(-e for e in reversed(a)),
}
def lm(f: Poly, key: Key) -> Mono:
return max(f, key=key)
def divides(a: Mono, b: Mono) -> bool:
return all(x <= y for x, y in zip(a, b))
def mmul(a: Mono, b: Mono) -> Mono:
return tuple(x + y for x, y in zip(a, b))
def mdiv(b: Mono, a: Mono) -> Mono:
return tuple(y - x for x, y in zip(a, b))
def shift(f: Poly, u: Mono, c: Fraction) -> Poly:
return {mmul(e, u): c * v for e, v in f.items()}
def axpy(f: Poly, c: Fraction, g: Poly) -> Poly:
h = dict(f)
for e, v in g.items():
w = h.get(e, Fraction(0)) - c * v
if w:
h[e] = w
else:
h.pop(e, None)
return h
def reduce_full(f: Poly, G: List[Poly], key: Key) -> Poly:
p, r = dict(f), {}
while p:
m = lm(p, key)
for g in G:
gm = lm(g, key)
if divides(gm, m):
p = axpy(p, p[m] / g[gm], shift(g, mdiv(m, gm), Fraction(1)))
break
else:
r[m] = p.pop(m)
return r
def reduced_basis(G: List[Poly], key: Key) -> List[Poly]:
minimal: List[Poly] = []
for g in sorted(G, key=lambda g: key(lm(g, key))):
if not any(divides(lm(h, key), lm(g, key)) for h in minimal):
minimal.append(g)
out = []
for i, g in enumerate(minimal):
m = lm(g, key)
tail = {e: v / g[m] for e, v in g.items() if e != m}
r = reduce_full(tail, minimal[:i] + minimal[i + 1:], key)
r[m] = Fraction(1)
out.append(r)
return sorted(out, key=lambda g: key(lm(g, key)))
def f5(F: List[Poly], order: str) -> Tuple[List[Poly], Dict[str, int]]:
"""Reduced Groebner basis of <F>, with counts of reductions and reductions to zero.
Precondition: F is a list of polynomials over Q in the same number of variables.
Cost: doubly exponential in the number of variables in the worst case.
"""
key = ORDERS[order]
def skey(s: Sig) -> Tuple[object, ...]:
return (s[0],) + key(s[1])
def smul(u: Mono, s: Sig) -> Sig:
return (s[0], mmul(u, s[1]))
F = [{e: Fraction(c) for e, c in f.items() if c} for f in F]
nonzero = [f for f in F if f]
if not nonzero:
return [], {"reductions": 0, "zero_reductions": 0}
n = len(next(iter(nonzero[0])))
one: Mono = tuple([0] * n)
B: List[Tuple[Sig, Poly]] = [] # @line L1
syz: List[Sig] = []
queue: List[Tuple[Sig, Poly]] = [((i, one), f) for i, f in enumerate(F)]
handled: Set[Sig] = set()
stats = {"reductions": 0, "zero_reductions": 0}
def regular_reduce(sigma: Sig, p: Poly) -> Poly:
p, r = dict(p), {}
while p:
m = lm(p, key)
for s, g in B:
gm = lm(g, key)
if divides(gm, m):
t = mdiv(m, gm)
if skey(smul(t, s)) < skey(sigma):
p = axpy(p, p[m] / g[gm], shift(g, t, Fraction(1)))
break
else:
r[m] = p.pop(m)
return r
def singular_reducible(sigma: Sig, r: Poly) -> bool:
m = lm(r, key)
for s, g in B:
gm = lm(g, key)
if divides(gm, m) and smul(mdiv(m, gm), s) == sigma:
return True
return False
while queue: # @line L2
queue.sort(key=lambda q: skey(q[0]), reverse=True)
sigma, p = queue.pop() # @line L3
if sigma in handled:
continue
handled.add(sigma)
i, u = sigma
if u == one: # @line L4
for s, g in B:
if s[0] < i:
syz.append((i, lm(g, key)))
if any(t[0] == i and divides(t[1], u) for t in syz): # @line L5
continue
stats["reductions"] += 1
r = regular_reduce(sigma, p) # @line L6
if not r: # @line L7
stats["zero_reductions"] += 1
syz.append(sigma)
continue
if singular_reducible(sigma, r): # @line L8
continue
rm = lm(r, key)
for s, g in B: # @line L9
gm = lm(g, key)
l = tuple(max(a, b) for a, b in zip(rm, gm))
ur, ug = mdiv(l, rm), mdiv(l, gm)
sr, sg = smul(ur, sigma), smul(ug, s)
if sr == sg:
continue
spoly = axpy(shift(r, ur, 1 / r[rm]), 1 / g[gm], shift(g, ug, Fraction(1)))
top = sr if skey(sr) > skey(sg) else sg
if top not in handled:
queue.append((top, spoly))
B.append((sigma, r)) # @line L10
return reduced_basis([g for _, g in B], key), stats # @line L11
test_f5_algorithm.py 286 lines
# SPDX-License-Identifier: MIT
# Part of https://algorithm.tamnd.com/f5-algorithm/. The authors are listed in prov/2QB464.jsonl.
"""Tests for the signature-based F5 against a plain Buchberger algorithm (product criterion only).
Each output is also checked by the defining property (every S-polynomial
reduces to zero), by containing the input in its ideal, and by ideal
membership through linear algebra on a Macaulay matrix in bounded degree.
"""
import random
from fractions import Fraction
from typing import Dict, List, Tuple
from f5_algorithm import f5
Mono = Tuple[int, ...]
Poly = Dict[Mono, Fraction]
KEYS = {
"lex": lambda a: tuple(a),
"grlex": lambda a: (sum(a),) + tuple(a),
"grevlex": lambda a: (sum(a),) + tuple(-e for e in reversed(a)),
}
def lead(f: Poly, order: str) -> Mono:
return max(f, key=KEYS[order])
def sub_scaled(f: Poly, g: Poly, c: Fraction, u: Mono) -> Poly:
h = dict(f)
for e, v in g.items():
e2 = tuple(a + b for a, b in zip(e, u))
w = h.get(e2, Fraction(0)) - c * v
if w:
h[e2] = w
else:
h.pop(e2, None)
return h
def remainder(f: Poly, G: List[Poly], order: str) -> Poly:
p, r = dict(f), {}
while p:
m = lead(p, order)
for g in G:
gm = lead(g, order)
if all(x <= y for x, y in zip(gm, m)):
p = sub_scaled(p, g, p[m] / g[gm], tuple(y - x for x, y in zip(gm, m)))
break
else:
r[m] = p.pop(m)
return r
def spoly(f: Poly, g: Poly, order: str) -> Poly:
a, b = lead(f, order), lead(g, order)
l = tuple(max(x, y) for x, y in zip(a, b))
s = sub_scaled({}, f, -1 / f[a], tuple(x - y for x, y in zip(l, a)))
return sub_scaled(s, g, 1 / g[b], tuple(x - y for x, y in zip(l, b)))
def buchberger_reduced(F: List[Poly], order: str) -> List[Poly]:
G = [dict(f) for f in F if f]
pairs = [(i, j) for j in range(len(G)) for i in range(j)]
def lcm_degree(p: Tuple[int, int]) -> int:
a, b = lead(G[p[0]], order), lead(G[p[1]], order)
return sum(max(x, y) for x, y in zip(a, b))
while pairs:
pairs.sort(key=lcm_degree, reverse=True)
i, j = pairs.pop()
a, b = lead(G[i], order), lead(G[j], order)
if all(x == 0 or y == 0 for x, y in zip(a, b)):
continue
h = remainder(spoly(G[i], G[j], order), G, order)
if h:
pairs += [(k, len(G)) for k in range(len(G))]
G.append(h)
minimal: List[Poly] = []
for g in sorted(G, key=lambda g: KEYS[order](lead(g, order))):
if not any(all(x <= y for x, y in zip(lead(h, order), lead(g, order))) for h in minimal):
minimal.append(g)
out = []
for i, g in enumerate(minimal):
m = lead(g, order)
tail = {e: v / g[m] for e, v in g.items() if e != m}
r = remainder(tail, minimal[:i] + minimal[i + 1:], order)
r[m] = Fraction(1)
out.append(r)
return sorted(out, key=lambda g: KEYS[order](lead(g, order)))
def canon(G: List[Poly]) -> List[List[Tuple[Mono, Fraction]]]:
return sorted(sorted(g.items()) for g in G)
def is_groebner(G: List[Poly], order: str) -> bool:
return all(not remainder(spoly(G[i], G[j], order), G, order)
for j in range(len(G)) for i in range(j))
def monomials_upto(n: int, e: int) -> List[Mono]:
out: List[Mono] = [()]
for _ in range(n):
out = [m + (k,) for m in out for k in range(e + 1)]
return [m for m in out if sum(m) <= e]
def echelon_size(rows: List[Poly]) -> int:
piv: Dict[Mono, Poly] = {}
for r in rows:
r = dict(r)
while r and max(r) in piv:
m = max(r)
q = piv[m]
r = sub_scaled(r, q, r[m] / q[m], tuple(0 for _ in m))
if r:
piv[max(r)] = r
return len(piv)
def macaulay_rows(G: List[Poly], n: int, e: int) -> List[Poly]:
rows = []
for g in G:
dg = max(sum(m) for m in g)
if e >= dg:
for m in monomials_upto(n, e - dg):
rows.append({tuple(a + b for a, b in zip(k, m)): v for k, v in g.items()})
return rows
def random_poly(rng: random.Random, n: int, d: int, terms: int) -> Poly:
f: Poly = {}
for _ in range(terms):
m = tuple(rng.randint(0, d) for _ in range(n))
if sum(m) <= d:
f[m] = f.get(m, Fraction(0)) + rng.randint(-4, 4)
return {m: c for m, c in f.items() if c}
def random_system(rng: random.Random, n: int, k: int, dmax: int) -> List[Poly]:
"""k polynomials with 2 to 5 terms, each with a term of its top degree."""
F = []
for _ in range(k):
dd = rng.randint(1, dmax)
top = [m for m in monomials_upto(n, dd) if sum(m) == dd]
rest = monomials_upto(n, dd)
f = {rng.choice(top): Fraction(rng.choice([-2, -1, 1, 2, 3]))}
for m in rng.sample(rest, min(len(rest), rng.randint(1, 4))):
f[m] = f.get(m, Fraction(0)) + rng.randint(-3, 3)
F.append({m: c for m, c in f.items() if c})
return F
def poly(terms: List[Tuple[Mono, int]]) -> Poly:
return {m: Fraction(c) for m, c in terms if c}
def check(F: List[Poly], order: str) -> List[Poly]:
G, _ = f5(F, order)
assert canon(G) == canon(buchberger_reduced(F, order))
assert is_groebner(G, order)
for f in F:
assert not remainder(f, G, order)
return G
CYCLIC3 = [poly([((1, 0, 0), 1), ((0, 1, 0), 1), ((0, 0, 1), 1)]),
poly([((1, 1, 0), 1), ((0, 1, 1), 1), ((1, 0, 1), 1)]),
poly([((1, 1, 1), 1), ((0, 0, 0), -1)])]
def test_worked_example_circle_hyperbola() -> None:
F = [poly([((2, 0), 1), ((0, 2), 1), ((0, 0), -4)]), poly([((1, 1), 1), ((0, 0), -1)])]
G, stats = f5(F, "grevlex")
assert G == [poly([((1, 1), 1), ((0, 0), -1)]),
poly([((2, 0), 1), ((0, 2), 1), ((0, 0), -4)]),
poly([((0, 3), 1), ((1, 0), 1), ((0, 1), -4)])]
assert stats == {"reductions": 3, "zero_reductions": 0}
def test_cyclic3_all_orders() -> None:
for order in KEYS:
check(CYCLIC3, order)
G, _ = f5(CYCLIC3, "lex")
assert canon(G) == canon([poly([((1, 0, 0), 1), ((0, 1, 0), 1), ((0, 0, 1), 1)]),
poly([((0, 2, 0), 1), ((0, 1, 1), 1), ((0, 0, 2), 1)]),
poly([((0, 0, 3), 1), ((0, 0, 0), -1)])])
def test_edge_cases() -> None:
assert f5([], "lex")[0] == []
assert f5([{}], "grevlex")[0] == []
one = poly([((0, 0), 1)])
assert f5([poly([((0, 0), 7)])], "grevlex")[0] == [one]
assert f5([poly([((1, 0), 1)]), poly([((1, 0), 1), ((0, 0), 1)])], "lex")[0] == [one]
single = poly([((2, 1), 3), ((0, 1), -6)])
assert f5([single], "grlex")[0] == [poly([((2, 1), 1), ((0, 1), -2)])]
assert f5([single, {}], "grlex")[0] == [poly([((2, 1), 1), ((0, 1), -2)])]
def test_many_reduction_steps() -> None:
F = [poly([((5, 0), 1), ((0, 1), -1)]), poly([((0, 5), 1), ((1, 0), -1), ((0, 0), -1)])]
G = check(F, "lex")
assert max(max(m[1] for m in g) for g in G) == 25
def test_random_against_buchberger() -> None:
rng = random.Random(4)
for trial in range(520):
n = rng.choice([2, 2, 3])
F = random_system(rng, n, rng.randint(2, 3), 3 if n == 2 else 2)
if trial % 10 == 0:
F.append(random_poly(rng, n, 2, rng.randint(1, 3)))
check(F, rng.choice(list(KEYS)))
def test_membership_agrees_with_linear_algebra() -> None:
rng = random.Random(11)
for _ in range(150):
n = 2
F = random_system(rng, n, 2, 3)
order = rng.choice(list(KEYS))
G, _ = f5(F, order)
Gdeg = buchberger_reduced(F, "grevlex")
e = max(max(sum(m) for m in g) for g in Gdeg) + 1
rows = macaulay_rows(Gdeg, n, e)
p: Poly = {}
for r in rng.sample(rows, min(3, len(rows))):
p = sub_scaled(p, r, Fraction(-rng.randint(-3, 3)), (0,) * n)
if rng.random() < 0.5:
p = sub_scaled(p, {rng.choice(monomials_upto(n, e)): Fraction(1)}, Fraction(-1), (0,) * n)
member = echelon_size(rows + [p]) == echelon_size(rows)
assert (not remainder(p, G, order)) == member
def hilbert_prefix(G: List[Poly], order: str, n: int, top: int) -> List[int]:
"""Number of standard monomials of each degree 0..top for the leading monomials of G."""
lms = [lead(g, order) for g in G]
counts = [0] * (top + 1)
for m in monomials_upto(n, top):
if not any(all(x <= y for x, y in zip(a, m)) for a in lms):
counts[sum(m)] += 1
return counts
def regular_series(degs: List[int], n: int, top: int) -> List[int]:
"""Coefficients of prod (1 - t^d) / (1 - t)^n up to t^top."""
c = [1] + [0] * top
for d in degs:
c = [c[k] - (c[k - d] if k >= d else 0) for k in range(top + 1)]
for _ in range(n):
for k in range(1, top + 1):
c[k] += c[k - 1]
return c
def random_homogeneous(rng: random.Random, n: int, d: int) -> Poly:
top = [m for m in monomials_upto(n, d) if sum(m) == d]
f = {m: Fraction(rng.randint(-5, 5)) for m in top}
return {m: c for m, c in f.items() if c}
def test_no_reduction_to_zero_on_regular_sequences() -> None:
rng = random.Random(2002)
checked = 0
for _ in range(40):
n = rng.choice([2, 3])
degs = [rng.randint(1, 3 if n == 2 else 2) for _ in range(n)]
F = [random_homogeneous(rng, n, d) for d in degs]
if any(not f for f in F):
continue
G, stats = f5(F, "grevlex")
top = sum(d - 1 for d in degs) + 2
if hilbert_prefix(G, "grevlex", n, top) == regular_series(degs, n, top):
assert stats["zero_reductions"] == 0
checked += 1
check(F, "grevlex")
assert checked >= 30
def test_fewer_zero_reductions_than_buchberger() -> None:
G, stats = f5(CYCLIC3, "grevlex")
assert stats["zero_reductions"] == 0
What the next level needs
- the tests did not run in this build