Wu's characteristic set method in Python
Decomposes polynomial systems into characteristic sets and proves geometric theorems mechanically.
- File
wu_characteristic_set_method.py, 148 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/wu-characteristic-set-method/. The authors are listed in prov/RYT59J.jsonl.
"""Wu's characteristic set method over a commutative ring R (Z, Q or Z/pZ).
A polynomial in R[x_0, ..., x_{n-1}] is a dict {exponent tuple: coefficient},
with the variable order x_0 < x_1 < ... < x_{n-1}. Only ring operations
(+, -, *) and zero tests on coefficients are used, so the code runs on int,
Fraction or any field element class. Cost model: operations in R; over Z the
coefficients can grow exponentially in the number of reductions.
charset(P) -> (C, P') a Ritt-Wu characteristic set of P
decompose(P) -> [C_1, ...] with Zero(P) = union of Zero(C_i / I_i)
wu_prove(H, g) -> (proved, remainder, nondegeneracy conditions)
"""
from __future__ import annotations
def clean(p):
return {e: c for e, c in p.items() if c != 0}
def add(p, q, s=1):
r = dict(p)
for e, c in q.items():
r[e] = r.get(e, 0) + c * s
return clean(r)
def mul(p, q):
r = {}
for e1, c1 in p.items():
for e2, c2 in q.items():
e = tuple(a + b for a, b in zip(e1, e2))
r[e] = r.get(e, 0) + c1 * c2
return clean(r)
def cls(p):
"""Class: index of the main variable, -1 for a constant."""
return max((i for e in p for i, k in enumerate(e) if k), default=-1)
def deg(p, i):
return max((e[i] for e in p), default=-1)
def rank(p):
c = cls(p)
return (c, deg(p, c) if c >= 0 else 0)
def coeff(p, i, k):
return {e[:i] + (0,) + e[i + 1:]: c for e, c in p.items() if e[i] == k}
def init(p):
c = cls(p)
return coeff(p, c, deg(p, c)) if c >= 0 else dict(p)
def shift(p, i, k):
return {e[:i] + (e[i] + k,) + e[i + 1:]: c for e, c in p.items()}
def prem(p, q):
"""Pseudo-remainder of p by q in the main variable of q."""
i = cls(q)
d, h, r = deg(q, i), init(q), dict(p)
while r and deg(r, i) >= d:
k = deg(r, i)
r = add(mul(h, r), mul(shift(coeff(r, i, k), i, k - d), q), -1)
return r
def prem_chain(p, C):
r = dict(p)
for b in reversed(C):
if not r:
break
r = prem(r, b)
return r
def reduced(p, C):
return all(deg(p, cls(b)) < deg(b, cls(b)) for b in C)
def is_constant(C):
return len(C) == 1 and cls(C[0]) < 0
def basic_set(P):
"""A basic set: an ascending set of lowest rank contained in P."""
Ps = sorted((p for p in P if p), key=rank)
if not Ps:
return []
if cls(Ps[0]) < 0:
return [Ps[0]] # a nonzero constant
B = []
for p in Ps:
if (not B or cls(p) > cls(B[-1])) and reduced(p, B):
B.append(p)
return B
def charset(P):
"""Return (C, P') with C a characteristic set of P and P' = P plus remainders."""
P = [p for p in P if p]
while True: # @line L1
B = basic_set(P) # @line L2
if is_constant(B): # @line L3
return B, P
R = []
for p in P: # @line L4
if not any(p is b for b in B):
r = prem_chain(p, B)
if r and r not in R:
R.append(r)
if not R: # @line L5
return B, P
P = P + R # @line L6
def decompose(P):
"""Wu's zero decomposition: Zero(P) is the union of Zero(C / I_C)."""
out, stack = [], [[p for p in P if p]]
while stack: # @line D1
Q = stack.pop()
C, Q2 = charset(Q) # @line D2
if is_constant(C): # @line D3
continue
out.append(C) # @line D4
for b in C: # @line D5
h = init(b)
if cls(h) >= 0:
stack.append(Q2 + [h]) # @line D6
return out
def wu_prove(H, g):
"""Prove that g vanishes on Zero(H) where all initials of C are nonzero."""
C, _ = charset(H) # @line W1
if is_constant(C):
return True, {}, [] # the hypotheses have no zeros
r = prem_chain(g, C) # @line W2
conds = [init(b) for b in C if cls(init(b)) >= 0]
return (not r), r, conds # @line W3
test_wu_characteristic_set_method.py 197 lines
# SPDX-License-Identifier: MIT
# Part of https://algorithm.tamnd.com/wu-characteristic-set-method/. The authors are listed in prov/RYT59J.jsonl.
"""Brute-force tests over Z/7Z, where every point of the affine space can be checked."""
import itertools
import random
from wu_characteristic_set_method import (add, mul, cls, deg, init, prem_chain,
basic_set, charset, decompose, wu_prove,
reduced, is_constant)
P7 = 7
class Mod:
__slots__ = ("v",)
def __init__(self, v):
self.v = v % P7
def _o(self, o):
return o.v if isinstance(o, Mod) else o
def __add__(self, o):
return Mod(self.v + self._o(o))
__radd__ = __add__
def __sub__(self, o):
return Mod(self.v - self._o(o))
def __rsub__(self, o):
return Mod(self._o(o) - self.v)
def __mul__(self, o):
return Mod(self.v * self._o(o))
__rmul__ = __mul__
def __neg__(self):
return Mod(-self.v)
def __eq__(self, o):
return self.v == self._o(o) % P7
def __ne__(self, o):
return not self == o
def __hash__(self):
return self.v
def ev(p, pt, m=None):
s = 0
for e, c in p.items():
t = c.v if isinstance(c, Mod) else c
for v, k in zip(pt, e):
t *= v ** k
s += t
return s % m if m else s
def rand_poly(rng, n, terms, maxdeg, mod=True):
p = {}
for _ in range(terms):
e = tuple(rng.randint(0, maxdeg) for _ in range(n))
if sum(e) > maxdeg + 1:
continue
c = rng.randint(1, P7 - 1)
p = add(p, {e: Mod(c) if mod else rng.choice([-2, -1, 1, 2])})
return p
def zeros(P, n, m=P7):
pts = itertools.product(range(m), repeat=n)
return {pt for pt in pts if all(ev(p, pt, m) == 0 for p in P)}
def quasi_zeros(C, n, m=P7):
h = [init(b) for b in C]
return {pt for pt in zeros(C, n, m) if all(ev(g, pt, m) != 0 for g in h)}
def check_ascending(C):
if is_constant(C):
return
for k, b in enumerate(C):
assert cls(b) >= 0
if k:
assert cls(b) > cls(C[k - 1])
assert reduced(b, C[:k])
def test_charset_properties_mod7():
rng = random.Random(11)
n = 3
for _ in range(60):
P = [rand_poly(rng, n, rng.randint(1, 3), 2) for _ in range(rng.randint(1, 3))]
P = [p for p in P if p]
C, P2 = charset(P)
check_ascending(C)
zP = zeros(P, n)
if is_constant(C):
assert not zP
continue
for p in P:
assert prem_chain(p, C) == {}
assert zP <= zeros(C, n)
assert quasi_zeros(C, n) <= zP
def test_decomposition_mod7():
rng = random.Random(12)
n = 3
sizes = []
for _ in range(40):
P = [rand_poly(rng, n, rng.randint(1, 3), 2) for _ in range(rng.randint(1, 3))]
P = [p for p in P if p]
D = decompose(P)
sizes.append(len(D))
union = set()
for C in D:
check_ascending(C)
union |= quasi_zeros(C, n)
assert union == zeros(P, n)
assert max(sizes) >= 2
def test_integers_pointwise():
rng = random.Random(13)
n = 3
for _ in range(30):
P = [rand_poly(rng, n, rng.randint(1, 3), 2, mod=False) for _ in range(2)]
P = [p for p in P if p]
C, _ = charset(P)
if is_constant(C):
continue
for p in P:
assert prem_chain(p, C) == {}
for m in (5, 11):
assert quasi_zeros(C, n, m) <= zeros(P, n, m)
assert zeros(P, n, m) <= zeros(C, n, m)
def var(n, i, k=1):
e = [0] * n
e[i] = k
return {tuple(e): 1}
def parallelogram():
# variables u1 < u2 < u3 < x1 < x2 < x3 < x4 (indices 0..6)
n = 7
u1, u2, u3, x1, x2, x3, x4 = (var(n, i) for i in range(n))
H = [add(x2, u3, -1), # C on the line y = u3
add(mul(add(x1, u1, -1), u3), mul(x2, u2), -1), # BC parallel to AD
add(mul(x4, x1), mul(x3, x2), -1), # O on AC
add(mul(x4, add(u2, u1, -1)), mul(u3, add(x3, u1, -1)), -1)] # O on BD
return n, H, (u1, u2, u3, x1, x2, x3, x4)
def test_parallelogram_theorem():
n, H, (u1, u2, u3, x1, x2, x3, x4) = parallelogram()
ok, r, conds = wu_prove(H, add(mul({(0,) * n: 2}, x3), x1, -1)) # 2 x3 = x1
assert ok and r == {}
ok2, _, _ = wu_prove(H, add(mul({(0,) * n: 2}, x4), x2, -1)) # 2 x4 = x2
assert ok2
bad, r3, _ = wu_prove(H, add(x3, x4, -1)) # false claim
assert not bad and r3
C, _ = charset(H)
assert [cls(b) for b in C] == [3, 4, 5, 6]
assert C[1] == add(x2, u3, -1)
assert C[0] == add(mul(u3, x1), mul(u3, add(u1, u2)), -1) # x1 = u1 + u2
# nondegeneracy conditions: u3, 2 u1 u3^2 and 2 u1 u3^3 (u1 + u2)
assert conds[0] == u3 and len(conds) == 3
def test_edges():
n = 2
one = {(0, 0): 1}
x, y = var(n, 0), var(n, 1)
assert charset([])[0] == []
assert decompose([]) == [[]]
assert charset([one])[0] == [one] # the unit ideal
assert decompose([one]) == []
assert decompose([x, add(x, one, -1)]) == [] # x = 0 and x = 1
C, _ = charset([mul(x, y)]) # a single generator
assert C == [mul(x, y)]
D = decompose([mul(x, y)]) # x y = 0: y = 0 (x != 0) or x = 0
assert len(D) == 2
assert charset([{}, x])[0] == [x] # the zero polynomial is dropped
assert basic_set([mul(y, y), add(y, x)]) == [add(y, x)]
ok, r, _ = wu_prove([one], x)
assert ok
# many reduction steps: y^6 - x by y^2 - x leaves x^3 - x
C = [add(mul(y, y), x, -1)]
assert prem_chain(add(var(n, 1, 6), x, -1), C) == add(var(n, 0, 3), x, -1)
What the next level needs
- the tests did not run in this build