Gröbner walk in Python
Moves a Gröbner basis across the cones of weight vectors to change monomial order.
- File
groebner_walk.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/groebner-walk/. The authors are listed in prov/2Y8G4M.jsonl.
"""The Gröbner walk of Collart, Kalkbrener and Mall over Q.
Polynomials in Q[x_1, ..., x_n] are dicts from exponent tuples to Fractions.
A monomial order is a matrix order: a list of integer rows, x^a > x^b when the
vector (r . a) for the rows r is lexicographically larger. The first row is the
weight vector of the order and must be nonnegative.
groebner_walk(G, Ms, Mt) takes the reduced Gröbner basis G of an ideal I for the
order Ms and returns (the reduced Gröbner basis of I for Mt, the list of weight
vectors visited). Ideals of any dimension are allowed.
Cost model: operations in Q. Each step costs one Gröbner basis computation for
an ideal of initial forms, one division and one interreduction; the number of
steps is the number of cones of the Gröbner fan met by the segment from the
first row of Ms to the first row of Mt. No useful worst-case bound is known.
"""
from fractions import Fraction
from math import gcd
from typing import Callable, Dict, List, Sequence, Tuple
Mono = Tuple[int, ...]
Poly = Dict[Mono, Fraction]
Matrix = List[List[int]]
def order_matrix(name: str, n: int) -> Matrix:
"""The matrices of lex, grlex and grevlex with x_1 > ... > x_n."""
unit = [[1 if j == i else 0 for j in range(n)] for i in range(n)]
if name == "lex":
return unit
if name == "grlex":
return [[1] * n] + unit[: n - 1]
if name == "grevlex":
return [[1 if j < n - i else 0 for j in range(n)] for i in range(n)]
raise ValueError(name)
def dot(w: Sequence[int], a: Mono) -> int:
return sum(x * y for x, y in zip(w, a))
def key_of(M: Matrix) -> Callable[[Mono], Tuple[int, ...]]:
return lambda a: tuple(dot(r, a) for r in M)
def lead(f: Poly, M: Matrix) -> Mono:
return max(f, key=key_of(M))
def divides(a: Mono, b: Mono) -> bool:
return all(x <= y for x, y in zip(a, b))
def add_multiple(f: Poly, g: Poly, c: Fraction, u: Mono) -> None:
"""f <- f - c x^u g, in place."""
for e, v in g.items():
e2 = tuple(x + y for x, y in zip(e, u))
w = f.get(e2, Fraction(0)) - c * v
if w:
f[e2] = w
else:
f.pop(e2, None)
def divide(f: Poly, F: List[Poly], M: Matrix) -> Tuple[List[Poly], Poly]:
"""Multivariate division of f by the list F for M: quotients and remainder."""
key = key_of(M)
leads = [max(g, key=key) for g in F]
q: List[Poly] = [{} for _ in F]
p, r = dict(f), {}
while p:
m = max(p, key=key)
for i, g in enumerate(F):
if divides(leads[i], m):
c = p[m] / g[leads[i]]
u = tuple(y - x for x, y in zip(leads[i], m))
q[i][u] = q[i].get(u, Fraction(0)) + c
add_multiple(p, g, c, u)
break
else:
r[m] = p.pop(m)
return q, r
def reduced_basis(F: List[Poly], M: Matrix) -> List[Poly]:
"""Reduced Gröbner basis for M by Buchberger's algorithm with the product criterion."""
key = key_of(M)
G = [dict(f) for f in F if f]
pairs = [(i, j) for j in range(len(G)) for i in range(j)]
while pairs:
pairs.sort(key=lambda p: -sum(max(x, y) for x, y in zip(lead(G[p[0]], M), lead(G[p[1]], M))))
i, j = pairs.pop()
a, b = lead(G[i], M), lead(G[j], M)
if all(x == 0 or y == 0 for x, y in zip(a, b)):
continue
l = tuple(max(x, y) for x, y in zip(a, b))
s: Poly = {}
add_multiple(s, G[i], -1 / G[i][a], tuple(x - y for x, y in zip(l, a)))
add_multiple(s, G[j], 1 / G[j][b], tuple(x - y for x, y in zip(l, b)))
h = divide(s, G, M)[1]
if h:
pairs += [(k, len(G)) for k in range(len(G))]
G.append(h)
return interreduce(G, M)
def interreduce(G: List[Poly], M: Matrix) -> List[Poly]:
"""Turn a Gröbner basis for M into the reduced one, sorted by leading monomial."""
key = key_of(M)
minimal: List[Poly] = []
for g in sorted(G, key=lambda g: key(lead(g, M))):
if not any(divides(lead(h, M), lead(g, M)) for h in minimal):
minimal.append(g)
out: List[Poly] = []
for i, g in enumerate(minimal):
m = lead(g, M)
tail = {e: v / g[m] for e, v in g.items() if e != m}
r = divide(tail, minimal[:i] + minimal[i + 1:], M)[1]
r[m] = Fraction(1)
out.append(r)
return sorted(out, key=lambda g: key(lead(g, M)))
def initial_form(g: Poly, w: Sequence[int]) -> Poly:
top = max(dot(w, e) for e in g)
return {e: c for e, c in g.items() if dot(w, e) == top}
def primitive(v: Sequence[Fraction]) -> List[int]:
"""Positive multiple of a nonnegative rational vector with coprime integer entries."""
den = 1
for x in v:
den = den * Fraction(x).denominator // gcd(den, Fraction(x).denominator)
ints = [int(Fraction(x) * den) for x in v]
g = 0
for x in ints:
g = gcd(g, x)
return [x // g for x in ints] if g else ints
def groebner_walk(G: List[Poly], Ms: Matrix, Mt: Matrix) -> Tuple[List[Poly], List[List[int]]]:
"""Convert the reduced Gröbner basis G for Ms into the reduced basis for Mt.
Precondition: G is the reduced Gröbner basis of I for Ms, and the first rows
of Ms and Mt are nonnegative. Cost: one small Gröbner basis per cone crossed.
"""
w, tau = list(Ms[0]), list(Mt[0]) # @line L1
M_old = Ms
G = [dict(g) for g in G]
path = [list(w)]
while True: # @line L2
M_new = [list(w)] + [list(r) for r in Mt] # @line L3
In = [initial_form(g, w) for g in G]
H = reduced_basis(In, M_new) # @line L4
lifted: List[Poly] = []
for h in H: # @line L5
q, r = divide(h, In, M_old)
assert not r
hp: Poly = {}
for qi, g in zip(q, G):
for u, c in qi.items():
add_multiple(hp, g, -c, u)
lifted.append(hp)
G = interreduce(lifted, M_new) # @line L6
if w == tau: # @line L7
return G, path
u = Fraction(1) # @line L8
for g in G:
a = lead(g, M_new)
for b in g:
diff = tuple(x - y for x, y in zip(a, b))
wt, tt = dot(w, diff), dot(tau, diff)
if tt < 0:
u = min(u, Fraction(wt, wt - tt))
w = primitive([(1 - u) * x + u * y for x, y in zip(w, tau)]) # @line L9
if u == 1:
w = list(tau)
M_old = M_new
path.append(list(w))
test_groebner_walk.py 237 lines
# SPDX-License-Identifier: MIT
# Part of https://algorithm.tamnd.com/groebner-walk/. The authors are listed in prov/2Y8G4M.jsonl.
"""Tests for the Gröbner walk against Buchberger's algorithm run directly in the target order.
The oracle uses its own order keys (not matrices), its own division and its own
reduced basis. Outputs are also checked by the defining property (every
S-polynomial reduces to zero) and by ideal membership of the input.
"""
import random
from fractions import Fraction
from typing import Dict, List, Tuple
from groebner_walk import groebner_walk, order_matrix
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 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 walk(F: List[Poly], src: str, dst: str, n: int):
G = buchberger_reduced(F, src)
return groebner_walk(G, order_matrix(src, n), order_matrix(dst, n))
def P(*terms) -> Poly:
return {tuple(e): Fraction(c) for c, e in terms}
def test_order_matrices_agree_with_keys():
rng = random.Random(1)
from groebner_walk import key_of
for n in (1, 2, 3, 4):
for name in KEYS:
k = key_of(order_matrix(name, n))
for _ in range(200):
a = tuple(rng.randint(0, 4) for _ in range(n))
b = tuple(rng.randint(0, 4) for _ in range(n))
assert (k(a) > k(b)) == (KEYS[name](a) > KEYS[name](b))
def test_worked_example():
# circle and hyperbola, grevlex to lex
F = [P((1, (2, 0)), (1, (0, 2)), (-4, (0, 0))), P((1, (1, 1)), (-1, (0, 0)))]
G, path = walk(F, "grevlex", "lex", 2)
assert canon(G) == canon(buchberger_reduced(F, "lex"))
assert canon(G) == canon([P((1, (1, 0)), (1, (0, 3)), (-4, (0, 1))),
P((1, (0, 4)), (-4, (0, 2)), (1, (0, 0)))])
assert path == [[1, 1], [3, 1], [1, 0]]
def test_twisted_cubic_positive_dimensional():
# ideal of the twisted cubic (t, t^2, t^3) in variables x > y > z
F = [P((1, (0, 1, 0)), (-1, (2, 0, 0))), P((1, (0, 0, 1)), (-1, (3, 0, 0)))]
for src in KEYS:
for dst in KEYS:
G, _ = walk(F, src, dst, 3)
assert canon(G) == canon(buchberger_reduced(F, dst))
def test_edge_cases():
for n in (1, 2, 3):
# unit ideal
G, _ = walk([P((3, (0,) * n))], "grevlex", "lex", n)
assert canon(G) == canon([P((1, (0,) * n))])
# zero ideal
G, _ = walk([], "grevlex", "lex", n)
assert G == []
# a single polynomial
f = P((2, (1,) + (0,) * (n - 1)), (5, (0,) * (n - 1) + (2,)), (1, (0,) * n))
G, _ = walk([f], "grevlex", "lex", n)
assert canon(G) == canon(buchberger_reduced([f], "lex"))
# same source and target
F = [P((1, (2, 1)), (-1, (0, 1))), P((1, (1, 2)), (1, (1, 0)))]
G, path = walk(F, "lex", "lex", 2)
assert canon(G) == canon(buchberger_reduced(F, "lex"))
assert path == [[1, 0]]
def test_many_cones():
# x^5 - y, y^5 - x - 1: lex basis of degree 25, several walls on the path
F = [P((1, (5, 0)), (-1, (0, 1))), P((1, (0, 5)), (-1, (1, 0)), (-1, (0, 0)))]
G, path = walk(F, "grevlex", "lex", 2)
assert canon(G) == canon(buchberger_reduced(F, "lex"))
assert len(path) >= 3
def random_system(rng: random.Random, n: int, k: int, d: int) -> List[Poly]:
out = []
while len(out) < k:
f = random_poly(rng, n, d, rng.randint(3, 8))
top = tuple(rng.choice(range(d + 1)) for _ in range(n))
if sum(top) <= d and sum(top) > 0:
f[top] = f.get(top, Fraction(0)) + 1
f = {m: c for m, c in f.items() if c}
if f:
out.append(f)
return out
def weight_order(w: Tuple[int, ...]):
"""Register the weight order w refined by lex in the oracle, and return its matrix."""
name = "w" + ",".join(map(str, w))
KEYS[name] = lambda a, w=w: (sum(x * y for x, y in zip(w, a)),) + tuple(a)
n = len(w)
return name, [list(w)] + [[1 if j == i else 0 for j in range(n)] for i in range(n)]
def test_weight_matrix_agrees_with_key():
rng = random.Random(3)
from groebner_walk import key_of
for _ in range(50):
w = tuple(rng.randint(0, 4) for _ in range(3))
name, M = weight_order(w)
for _ in range(40):
a = tuple(rng.randint(0, 4) for _ in range(3))
b = tuple(rng.randint(0, 4) for _ in range(3))
assert (key_of(M)(a) > key_of(M)(b)) == (KEYS[name](a) > KEYS[name](b))
def test_random_against_direct_buchberger():
rng = random.Random(31337)
walked = 0
for trial in range(540):
n = rng.choice((2, 2, 3))
d = 3 if n == 2 else 2
k = rng.randint(1, n)
F = random_system(rng, n, k, d)
src = rng.choice(("grevlex", "grlex", "lex"))
Ms = order_matrix(src, n)
kind = trial % 3
if kind == 0:
dst, Mt = "lex", order_matrix("lex", n)
elif kind == 1:
dst, Mt = "grevlex", order_matrix("grevlex", n)
else:
dst, Mt = weight_order(tuple(rng.randint(1, 5) for _ in range(n)))
G, path = groebner_walk(buchberger_reduced(F, src), Ms, Mt)
expected = buchberger_reduced(F, dst)
assert canon(G) == canon(expected)
assert is_groebner(G, dst)
assert all(not remainder(f, G, dst) for f in F)
assert path[0] == Ms[0] and path[-1] == Mt[0]
walked += len(path) >= 3
assert walked >= 80 # 95 with this seed: at least one wall inside the path
What the next level needs
- the tests did not run in this build