Valiant's algorithm in Python
Reduces CFG recognition to Boolean matrix multiplication, giving sub-cubic time.
- File
valiant_algorithm.py, 252 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/valiant-algorithm/. The authors are listed in prov/44TRAH.jsonl.
"""Valiant's algorithm: context-free recognition by Boolean matrix multiplication.
A grammar is a list of productions (A, rhs), where rhs is a tuple of symbols.
The nonterminals are the left sides and the start symbol; every other symbol
is a terminal. VALIANT needs Chomsky normal form: A -> B C with B and C
nonterminals, A -> a with a terminal, and S -> () only for a start symbol S
that is on no right side. to_cnf converts any grammar.
The table T is kept as one Boolean matrix per nonterminal A: rows[A][i] is an
int whose bit j is set when A derives w[i:j]. Positions are 0..n, padded to
N, a power of two. The recursion follows Okhotin's formulation; its depth is
log2 N, and every product goes through one pluggable function bmm.
"""
def bmm_rows(x, y):
"""Boolean product of x (r rows, each an int over k columns) and y (k rows,
each an int over c columns): row i of the result is the OR of the rows
y[t] for the set bits t of x[i]."""
out = []
for row in x:
acc, t = 0, 0
while row:
if row & 1:
acc |= y[t]
row >>= 1
t += 1
out.append(acc)
return out
def is_cnf(prods, start):
nts = {a for a, _ in prods} | {start}
for a, rhs in prods:
if len(rhs) == 2 and rhs[0] in nts and rhs[1] in nts:
if start in rhs and (start, ()) in prods:
return False
elif len(rhs) == 1 and rhs[0] not in nts:
pass
elif len(rhs) == 0 and a == start:
pass
else:
return False
if (start, ()) in prods and any(start in rhs for _, rhs in prods):
return False
return True
class Valiant:
"""Fills the recognition table of w for a grammar in CNF."""
def __init__(self, prods, start, w, bmm=None):
prods = [(a, tuple(r)) for a, r in prods]
if not is_cnf(prods, start):
raise ValueError("the grammar is not in Chomsky normal form")
self.start, self.w, self.n = start, list(w), len(w)
self.bmm = bmm or bmm_rows
self.calls = 0 # number of calls of bmm
self.cells = 0 # number of cells finished at line L9
self.products = {} # block size -> number of block products (lines L13 to L15)
self.nts = sorted({a for a, _ in prods} | {start})
self.pairs = sorted({rhs for _, rhs in prods if len(rhs) == 2})
self.lhs = {bc: sorted({a for a, r in prods if r == bc}) for bc in self.pairs}
self.by_term = {}
for a, rhs in prods:
if len(rhs) == 1:
self.by_term.setdefault(rhs[0], set()).add(a)
self.accepts_empty = (start, ()) in prods
size = 1
while size < self.n + 1:
size *= 2
self.N = size
self.rows = {a: [0] * size for a in self.nts} # T, one matrix per nonterminal
self.prod = {bc: [0] * size for bc in self.pairs} # P, one matrix per pair (B, C)
for i, c in enumerate(self.w): # @line L1
for a in self.by_term.get(c, ()):
self.rows[a][i] |= 1 << (i + 1)
self.compute(0, size) # @line L2
def compute(self, l, m):
"""Fill T for all cells l <= i < j < m."""
if m - l >= 2: # @line L5
h = (l + m) // 2
self.compute(l, h) # @line L6
self.compute(h, m)
self.complete(l, h, h, m) # @line L7
def _block(self, mat, r0, r1, c0, c1):
mask = (1 << (c1 - c0)) - 1
return [(mat[i] >> c0) & mask for i in range(r0, r1)]
def _add_product(self, r0, r1, k0, k1, c0, c1):
"""P[r0:r1, c0:c1] |= T[r0:r1, k0:k1] x T[k0:k1, c0:c1], for every pair (B, C)."""
self.products[r1 - r0] = self.products.get(r1 - r0, 0) + 1
for b, c in self.pairs:
x = self._block(self.rows[b], r0, r1, k0, k1)
if not any(x):
continue
y = self._block(self.rows[c], k0, k1, c0, c1)
if not any(y):
continue
self.calls += 1
z = self.bmm(x, y)
target = self.prod[(b, c)]
for t, row in enumerate(z):
target[r0 + t] |= row << c0
def complete(self, l, m, l2, m2):
"""Fill T for rows l..m-1 and columns l2..m2-1, where m <= l2. Needs: T is
done inside [l, m) and inside [l2, m2), and P holds every split k in [m, l2)."""
if m - l == 1: # @line L8
i, j = l, l2
if j > i + 1: # @line L9
bit = 1 << j
for bc in self.pairs:
if self.prod[bc][i] & bit:
for a in self.lhs[bc]:
self.rows[a][i] |= bit
self.cells += 1
return # @line L10
a = (l + m) // 2 # @line L11
b = (l2 + m2) // 2
self.complete(a, m, l2, b) # @line L12
self._add_product(l, a, a, m, l2, b) # @line L13
self.complete(l, a, l2, b)
self._add_product(a, m, l2, b, b, m2) # @line L14
self.complete(a, m, b, m2)
self._add_product(l, a, a, m, b, m2) # @line L15
self._add_product(l, a, l2, b, b, m2)
self.complete(l, a, b, m2) # @line L16
def derives(self, a, i, j):
"""True when nonterminal a derives w[i:j], for 0 <= i < j <= n."""
return bool((self.rows[a][i] >> j) & 1)
def table(self):
"""The set of (A, i, j) with i < j and A =>* w[i:j]."""
out = set()
for a in self.nts:
for i in range(self.n):
row = self.rows[a][i]
for j in range(i + 1, self.n + 1):
if (row >> j) & 1:
out.add((a, i, j))
return out
def accepts(self):
if self.n == 0: # @line L3
return self.accepts_empty
return self.derives(self.start, 0, self.n) # @line L4
def recognize(prods, start, w, bmm=None):
"""True when w is in L(G), for G in Chomsky normal form."""
return Valiant(prods, start, w, bmm).accepts()
def to_cnf(prods, start):
"""An equivalent grammar in Chomsky normal form, and its start symbol."""
prods = [(a, tuple(r)) for a, r in prods]
nts = {a for a, _ in prods} | {start}
used = set(nts) | {x for _, r in prods for x in r}
count = [0]
def fresh(base):
while True:
count[0] += 1
name = "<%s%d>" % (base, count[0])
if name not in used:
used.add(name)
return name
s0 = fresh("S")
rules = [(s0, (start,))] + prods
nts.add(s0)
# TERM: terminals inside long right sides get their own nonterminal
term = {}
step = []
for a, rhs in rules:
if len(rhs) >= 2:
new = []
for x in rhs:
if x not in nts:
if x not in term:
term[x] = fresh("T")
x = term[x]
new.append(x)
rhs = tuple(new)
step.append((a, rhs))
for x, t in term.items():
step.append((t, (x,)))
nts.add(t)
# BIN: right sides of length at most 2
rules = []
for a, rhs in step:
while len(rhs) > 2:
t = fresh("B")
nts.add(t)
rules.append((a, (rhs[0], t)))
a, rhs = t, rhs[1:]
rules.append((a, rhs))
# DEL: remove empty right sides, keeping S0 -> () when S0 is nullable
nullable = set()
changed = True
while changed:
changed = False
for a, rhs in rules:
if a not in nullable and all(x in nullable for x in rhs):
nullable.add(a)
changed = True
step = set()
for a, rhs in rules:
if len(rhs) == 2:
step.add((a, rhs))
if rhs[1] in nullable:
step.add((a, rhs[:1]))
if rhs[0] in nullable:
step.add((a, rhs[1:]))
elif len(rhs) == 1:
step.add((a, rhs))
# UNIT: replace A -> B by A -> beta for every non-unit B -> beta
unit = {a: {a} for a in nts}
changed = True
while changed:
changed = False
for a, rhs in step:
if len(rhs) == 1 and rhs[0] in nts:
for x in list(unit):
if a in unit[x] and rhs[0] not in unit[x]:
unit[x].add(rhs[0])
changed = True
out = set()
for x in nts:
for a, rhs in step:
if a in unit[x] and not (len(rhs) == 1 and rhs[0] in nts):
out.add((x, rhs))
# keep only productions whose nonterminals derive some nonempty string, so
# that every nonterminal on a right side has a production
gen = set()
changed = True
while changed:
changed = False
for a, rhs in out:
if a not in gen and all(x in gen or x not in nts for x in rhs):
gen.add(a)
changed = True
out = {(a, rhs) for a, rhs in out if all(x in gen or x not in nts for x in rhs)}
if s0 in nullable:
out.add((s0, ()))
return sorted(out), s0
test_valiant_algorithm.py 338 lines
# SPDX-License-Identifier: MIT
# Part of https://algorithm.tamnd.com/valiant-algorithm/. The authors are listed in prov/44TRAH.jsonl.
import itertools
import random
import pytest
from valiant_algorithm import Valiant, bmm_rows, is_cnf, recognize, to_cnf
def grammar(text):
prods = []
for rule in text.split(";"):
a, rhs = rule.split("->")
for alt in rhs.split("|"):
prods.append((a.strip(), tuple(x for x in alt.split() if x != "eps")))
return prods
def spans(prods, w):
"""Oracle for any grammar: all (X, i, j) with X =>* w[i:j], by a least fixpoint
that tries every split of every right side."""
nts = {a for a, _ in prods}
n, known = len(w), set()
def fits(rhs, i, j):
if not rhs:
return i == j
for cut in itertools.combinations_with_replacement(range(i, j + 1), len(rhs) - 1):
b = (i, *cut, j)
ok = True
for t, x in enumerate(rhs):
if x in nts:
ok = (x, b[t], b[t + 1]) in known
else:
ok = b[t + 1] == b[t] + 1 and w[b[t]] == x
if not ok:
break
if ok:
return True
return False
changed = True
while changed:
changed = False
for a, rhs in prods:
for i in range(n + 1):
for j in range(i, n + 1):
if (a, i, j) not in known and fits(rhs, i, j):
known.add((a, i, j))
changed = True
return known
def cyk_table(prods, w):
"""Plain CYK on a CNF grammar: all (A, i, j) with i < j and A =>* w[i:j]."""
n = len(w)
t = {}
for i, c in enumerate(w):
t[(i, i + 1)] = {a for a, r in prods if r == (c,)}
for length in range(2, n + 1):
for i in range(n - length + 1):
j = i + length
cell = set()
for k in range(i + 1, j):
for a, r in prods:
if len(r) == 2 and r[0] in t[(i, k)] and r[1] in t[(k, j)]:
cell.add(a)
t[(i, j)] = cell
return {(a, i, j) for (i, j), cell in t.items() for a in cell}
def random_cnf(rng):
names = ["S", "A", "B", "C"][: rng.randint(1, 4)]
prods = set()
for _ in range(rng.randint(2, 9)):
a = rng.choice(names)
if rng.random() < 0.6:
prods.add((a, (rng.choice(names), rng.choice(names))))
else:
prods.add((a, (rng.choice("ab"),)))
# a name with no production would be a terminal; drop the binary rules that use one
while True:
lhs = {a for a, _ in prods}
keep = {(a, r) for a, r in prods if len(r) == 1 or set(r) <= lhs}
if keep == prods:
return sorted(prods)
prods = keep
def random_grammar(rng):
"""Any shape: empty right sides, unit productions, left recursion, long right sides."""
names = ["S", "A", "B", "C"][: rng.randint(1, 4)]
prods = set()
for _ in range(rng.randint(2, 7)):
a = rng.choice(names)
k = rng.choice([0, 1, 1, 2, 2, 3])
prods.add((a, tuple(rng.choice(names + ["a", "b"]) for _ in range(k))))
prods.add(("S", tuple(rng.choice(names + ["a", "b"]) for _ in range(rng.randint(1, 3)))))
return sorted(prods)
def words(maxn):
for n in range(maxn + 1):
for t in itertools.product("ab", repeat=n):
yield "".join(t)
HMU = grammar("S -> A B | B C ; A -> B A | a ; B -> C C | b ; C -> A B | a")
def test_hmu_example_table():
v = Valiant(HMU, "S", "baaba")
assert v.N == 8
assert v.accepts()
# the CYK table of the page, row by row (0-based spans)
expect = {
(0, 1): "B", (1, 2): "AC", (2, 3): "AC", (3, 4): "B", (4, 5): "AC",
(0, 2): "AS", (1, 3): "B", (2, 4): "CS", (3, 5): "AS",
(0, 3): "", (1, 4): "B", (2, 5): "B",
(0, 4): "", (1, 5): "ACS",
(0, 5): "ACS",
}
got = v.table()
for (i, j), cell in expect.items():
assert {a for a in "ABCS" if (a, i, j) in got} == set(cell)
assert got == cyk_table(HMU, "baaba")
# 28 cells above the diagonal of the 8 x 8 matrix, 7 of them terminal cells
assert v.cells == 21
assert v.products == {1: 24, 2: 4}
def test_cell_order_trace():
"""The order in which line L9 finishes the cells of the 8 x 8 example."""
order = []
class Traced(Valiant):
def complete(self, l, m, l2, m2):
if m - l == 1 and l2 > l + 1:
order.append((l, l2))
super().complete(l, m, l2, m2)
Traced(HMU, "S", "baaba")
assert order == [(0, 2), (1, 3), (0, 3), (4, 6), (5, 7), (4, 7), (2, 4), (3, 5), (2, 5),
(1, 4), (0, 4), (1, 5), (0, 5), (3, 6), (2, 6), (3, 7), (2, 7),
(1, 6), (0, 6), (1, 7), (0, 7)]
def test_products_for_cell_0_5():
log = []
class Logged(Valiant):
def _add_product(self, r0, r1, k0, k1, c0, c1):
if r0 <= 0 < r1 and c0 <= 5 < c1:
log.append((r0, r1, k0, k1, c0, c1))
super()._add_product(r0, r1, k0, k1, c0, c1)
v = Logged(HMU, "S", "baaba")
# k = 2, 3 by one 2 x 2 block product, then k = 1 and k = 4 by 1 x 1 products
assert log == [(0, 2, 2, 4, 4, 6), (0, 1, 1, 2, 5, 6), (0, 1, 4, 5, 5, 6)]
pairs = {bc for bc in v.pairs if v.prod[bc][0] >> 5 & 1}
# the pairs that are right sides of productions: (B, A) and (B, C) from k = 1, (A, B) from k = 2
assert pairs == {("B", "A"), ("B", "C"), ("A", "B")}
assert v.calls == 15
def test_invariant_every_split_is_known_when_a_cell_is_finished():
rng = random.Random(7)
checked = 0
class Checked(Valiant):
def complete(self, l, m, l2, m2):
if m - l == 1 and l2 > l + 1:
i, j = l, l2
need = set()
for k in range(i + 1, j):
for b, c in self.pairs:
if self.rows[b][i] >> k & 1 and self.rows[c][k] >> j & 1:
need.add((b, c))
have = {bc for bc in self.pairs if self.prod[bc][i] >> j & 1}
assert have == need
nonlocal checked
checked += 1
super().complete(l, m, l2, m2)
for _ in range(60):
g = random_cnf(rng)
w = "".join(rng.choice("ab") for _ in range(rng.randint(2, 20)))
Checked(g, "S", w)
assert checked > 3000
def test_random_cnf_tables_against_cyk_and_fixpoint():
rng = random.Random(2024)
compared, accepted = 0, 0
for _ in range(400):
g = random_cnf(rng)
w = "".join(rng.choice("ab") for _ in range(rng.randint(0, 12)))
v = Valiant(g, "S", w)
got = v.table()
assert got == cyk_table(g, w)
assert got == {(a, i, j) for a, i, j in spans(g, w) if i < j}
compared += 1
accepted += v.accepts()
assert compared >= 400 and accepted >= 50
def test_random_general_grammars_through_cnf():
"""Empty and unit productions, left recursion, ambiguity and the empty input,
against the fixpoint oracle on the original grammar."""
rng = random.Random(99)
compared, accepted, empty_accepted = 0, 0, 0
for _ in range(250):
g = random_grammar(rng)
cnf, s0 = to_cnf(g, "S")
assert is_cnf(cnf, s0)
for w in words(5):
want = ("S", 0, len(w)) in spans(g, w)
assert recognize(cnf, s0, w) == want, (g, w)
compared += 1
accepted += want
empty_accepted += want and w == ""
assert compared >= 15000
assert accepted >= 500 and empty_accepted >= 10
def test_named_grammars():
cases = [
("E -> E + E | E * E | ( E ) | x", "E", ["x", "x+x*x", "(x+x)*x", "((x))"], ["", "x+", "()", "x x"]),
("S -> S S | ( S ) | eps", "S", ["", "()", "(()())", "()()(())"], ["(", ")(", "(()"]),
("S -> A ; A -> B ; B -> S | a", "S", ["a"], ["", "aa"]),
("S -> a S b | eps", "S", ["", "ab", "aaabbb"], ["a", "ba", "aab"]),
("S -> S a | a", "S", ["a", "aaaa"], ["", "b"]),
]
for text, start, yes, no in cases:
g = grammar(text)
cnf, s0 = to_cnf(g, start)
for w in yes:
assert recognize(cnf, s0, list(w)), (text, w)
for w in no:
assert not recognize(cnf, s0, list(w)), (text, w)
def test_pluggable_multiplication():
calls = []
def cell_by_cell(x, y):
calls.append((len(x), len(y)))
out = []
for row in x:
acc = 0
for t in range(len(y)):
if row >> t & 1:
acc |= y[t]
out.append(acc)
return out
rng = random.Random(5)
for _ in range(80):
g = random_cnf(rng)
w = "".join(rng.choice("ab") for _ in range(rng.randint(0, 15)))
v1 = Valiant(g, "S", w)
v2 = Valiant(g, "S", w, bmm=cell_by_cell)
assert v1.table() == v2.table()
assert v1.calls == v2.calls
assert calls and all(r == k for r, k in calls)
assert bmm_rows([0b011, 0b100], [0b01, 0b10, 0b11]) == [0b11, 0b11]
def test_product_counts_follow_the_recurrence():
def expected(size):
# products of each block size made by COMPUTE(0, size)
def comp(s):
out = {}
if s >= 2:
for _ in range(4):
out[s // 2] = out.get(s // 2, 0) + 1
for _ in range(4):
for t, c in comp(s // 2).items():
out[t] = out.get(t, 0) + c
return out
def compute(s):
out = {}
if s >= 2:
for part in (compute(s // 2), compute(s // 2), comp(s // 2)):
for t, c in part.items():
out[t] = out.get(t, 0) + c
return out
return compute(size)
for n in (1, 3, 7, 15, 31, 40):
v = Valiant(HMU, "S", "ab" * (n // 2) + "a" * (n % 2))
assert v.products == expected(v.N)
for t, c in v.products.items():
# the closed form of the page: (N / t) (N / (2 t) - 1) products of size t
assert c == (v.N // t) * (v.N // (2 * t) - 1)
def test_rejects_non_cnf():
with pytest.raises(ValueError):
Valiant(grammar("S -> a S b | eps"), "S", "ab")
with pytest.raises(ValueError):
Valiant(grammar("S -> A ; A -> a"), "S", "a")
with pytest.raises(ValueError):
Valiant(grammar("S -> S S | eps"), "S", "")
def test_empty_input():
assert recognize([("S", ()), ("S", ("a",))], "S", "")
assert not recognize([("S", ("a",))], "S", "")
v = Valiant([("S", ("a",))], "S", "")
assert v.N == 1 and v.table() == set()
def test_long_input():
# Dyck words in CNF: S -> S S | L R | L X, X -> S R
g = grammar("S -> S S | L R | L X ; X -> S R ; L -> ( ; R -> )")
rng = random.Random(11)
w = []
depth = 0
while len(w) < 250:
if depth and (rng.random() < 0.5 or depth > 40):
w.append(")")
depth -= 1
else:
w.append("(")
depth += 1
w += [")"] * depth
v = Valiant(g, "S", w)
assert v.accepts()
assert v.N == 512
assert v.table() == cyk_table(g, w)
assert not recognize(g, "S", w[:-1])
deep = ["("] * 150 + [")"] * 150
assert recognize(g, "S", deep)
What the next level needs
- the tests did not run in this build