Block Wiedemann algorithm in Python
Block version of Wiedemann's method that exploits parallel matrix-vector products for sparse matrices over finite fields.
- File
block_wiedemann_algorithm.py, 117 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/block-wiedemann-algorithm/. The authors are listed in prov/0EQBTW.jsonl.
import random
from typing import Sequence
SparseRows = Sequence[Sequence[tuple[int, int]]] # row i is a list of (column, value)
Block = list[list[int]] # n x s, stored by rows
def apply_block(rows: SparseRows, v: Block, p: int) -> Block:
"""A V modulo p for a sparse A and a dense n x s block V: s * nnz(A) multiplications."""
s = len(v[0]) if v else 0
return [[sum(a * v[j][c] for j, a in r) % p for c in range(s)] for r in rows]
def apply_vector(rows: SparseRows, u: Sequence[int], p: int) -> list[int]:
"""A u modulo p."""
return [sum(a * u[j] for j, a in r) % p for r in rows]
def minimal_generator(seq: Sequence[Block], s: int, p: int) -> list[tuple[int, list[list[int]]]]:
"""Right generators of a sequence of s x s matrices S_0, ..., S_{L-1} over F_p.
Returns pairs (d, f) with f = [f_0, ..., f_d], vectors of length s, such that
sum_j S_{i+j} f_j = 0 for 0 <= i < L - d, sorted by d. Computed as an order-L
approximant basis of the 2s x s matrix [S(x)^T; -I] by the iterative M-Basis method,
with the degree shift (0, ..., 0, 1, ..., 1): O(s^3 L^2) operations modulo p.
"""
big_l = len(seq)
m = 2 * s
g = []
for k in range(big_l):
top = [[seq[k][c][i] % p for c in range(s)] for i in range(s)]
bot = [[(p - 1) if (i == c and k == 0) else 0 for c in range(s)] for i in range(s)]
g.append(top + bot)
poly = [[[int(i == j)] for j in range(m)] for i in range(m)] # poly[i][j]: coefficients of P_ij
delta = [0] * s + [1] * s
for k in range(big_l):
disc = [] # coefficient of x^k in P G
for i in range(m):
acc = [0] * s
for j in range(m):
for t, cf in enumerate(poly[i][j]):
if cf and t <= k:
gr = g[k - t][j]
for c in range(s):
acc[c] += cf * gr[c]
disc.append([x % p for x in acc])
pivots: list[tuple[int, int]] = []
for i in sorted(range(m), key=lambda r: delta[r]):
for j, c in pivots:
if disc[i][c]:
f = disc[i][c] * pow(disc[j][c], -1, p) % p
disc[i] = [(x - f * y) % p for x, y in zip(disc[i], disc[j])]
for q in range(m):
a, b = poly[i][q], poly[j][q]
a = a + [0] * (len(b) - len(a))
poly[i][q] = [(a[t] - f * (b[t] if t < len(b) else 0)) % p for t in range(len(a))]
c = next((c for c in range(s) if disc[i][c]), None)
if c is not None:
pivots.append((i, c))
for i, _ in pivots:
poly[i] = [[0] + pol for pol in poly[i]]
delta[i] += 1
out = []
for i in sorted(range(m), key=lambda r: delta[r]):
a = poly[i][:s]
if not any(any(pol) for pol in a):
continue
d = delta[i]
f = [[a[c][d - j] if d - j < len(a[c]) else 0 for c in range(s)] for j in range(d + 1)]
out.append((d, f))
return out
def block_wiedemann_kernel(rows: SparseRows, n: int, p: int, s: int = 4,
seed: int = 0, tries: int = 10) -> list[int] | None:
"""A nonzero u with A u = 0 over F_p for a sparse n x n matrix A, or None.
Las Vegas: a returned vector is always checked, and None means that no kernel vector
was found in the given number of tries (always the case for a nonsingular A).
Cost per try: about 3n/s products of A with an n x s block, 2 n^2 s operations for the
projections, and O(s n^2) for the generator, all modulo p.
"""
rng = random.Random(seed)
for _ in range(tries):
if n == 0:
return None
x = [[rng.randrange(p) for _ in range(s)] for _ in range(n)] # @line L1
z = [[rng.randrange(p) for _ in range(s)] for _ in range(n)]
v = apply_block(rows, z, p)
big_l = 2 * (-(-n // s)) + 4 # @line L2
seq = []
for _ in range(big_l): # @line L3
seq.append([[sum(x[r][a] * v[r][c] for r in range(n)) % p for c in range(s)] for a in range(s)]) # @line L4
v = apply_block(rows, v, p)
gens = minimal_generator(seq, s, p) # @line L5
for d, f in gens: # @line L6
if d > big_l - (-(-n // s)):
break
e = next(j for j in range(d + 1) if any(f[j])) # @line L7
u = [0] * n
w = z
for j in range(e, d + 1):
for r in range(n):
u[r] = (u[r] + sum(w[r][c] * f[j][c] for c in range(s))) % p
if j < d:
w = apply_block(rows, w, p)
for _ in range(e + 2): # @line L8
if not any(u):
break
au = apply_vector(rows, u, p)
if not any(au): # @line L9
return u
u = au # @line L10
return None # @line L11
test_block_wiedemann_algorithm.py 147 lines
# SPDX-License-Identifier: MIT
# Part of https://algorithm.tamnd.com/block-wiedemann-algorithm/. The authors are listed in prov/0EQBTW.jsonl.
import random
from block_wiedemann_algorithm import apply_vector, block_wiedemann_kernel, minimal_generator
def rank_mod_p(dense, p):
m = [row[:] for row in dense]
n = len(m[0]) if m else 0
r = 0
for c in range(n):
piv = next((i for i in range(r, len(m)) if m[i][c] % p), None)
if piv is None:
continue
m[r], m[piv] = m[piv], m[r]
inv = pow(m[r][c], -1, p)
for i in range(len(m)):
if i != r and m[i][c] % p:
f = m[i][c] * inv % p
m[i] = [(x - f * y) % p for x, y in zip(m[i], m[r])]
r += 1
return r
def to_dense(rows, n):
d = [[0] * n for _ in rows]
for i, r in enumerate(rows):
for j, a in r:
d[i][j] = a
return d
def random_sparse(rng, n, p, per_row):
return [sorted((j, rng.randrange(1, p)) for j in rng.sample(range(n), min(per_row, n))) for _ in range(n)]
def make_singular(rng, rows, n, p):
rows = [list(r) for r in rows]
if n == 1:
return [[]]
i, k = rng.sample(range(n), 2)
c = rng.randrange(1, p)
di, dk = dict(rows[i]), dict(rows[k])
comb = {j: (di.get(j, 0) + c * dk.get(j, 0)) % p for j in set(di) | set(dk)}
rows[i] = sorted((j, v) for j, v in comb.items() if v)
rows[k] = rows[k]
t = rng.randrange(n)
if t not in (i, k):
rows[t] = list(rows[i]) # sometimes a second dependency
return rows
def check_kernel(rows, n, p, u):
assert len(u) == n and any(u)
assert all(0 <= x < p for x in u)
assert not any(apply_vector(rows, u, p))
def test_generator_property() -> None:
rng = random.Random(1994)
for _ in range(40):
p = rng.choice([2, 3, 101, 65521])
s = rng.randint(1, 3)
n = rng.randint(1, 10)
a = random_sparse(rng, n, p, 3)
x = [[rng.randrange(p) for _ in range(s)] for _ in range(n)]
v = [[rng.randrange(p) for _ in range(s)] for _ in range(n)]
seq = []
for _ in range(2 * (-(-n // s)) + 4):
seq.append([[sum(x[r][i] * v[r][c] for r in range(n)) % p for c in range(s)] for i in range(s)])
v = [[sum(val * v[j][c] for j, val in row) % p for c in range(s)] for row in a]
gens = minimal_generator(seq, s, p)
assert gens and [d for d, _ in gens] == sorted(d for d, _ in gens)
big_l = len(seq)
for d, f in gens:
assert any(any(fj) for fj in f)
for i in range(big_l - d):
for r in range(s):
assert sum(seq[i + j][r][c] * f[j][c] for j in range(d + 1) for c in range(s)) % p == 0
def test_edges() -> None:
assert block_wiedemann_kernel([], 0, 7) is None
u = block_wiedemann_kernel([[]], 1, 7)
assert u is not None and u[0] != 0
assert block_wiedemann_kernel([[(0, 3)]], 1, 7) is None
zero = [[] for _ in range(5)]
check_kernel(zero, 5, 2, block_wiedemann_kernel(zero, 5, 2, s=2))
a = [[(1, 1)], [(0, 1)], []] # zero row and zero column
check_kernel(a, 3, 5, block_wiedemann_kernel(a, 3, 5, s=2))
nilp = [[(i + 1, 1)] if i + 1 < 6 else [] for i in range(6)] # shift matrix, A^5 != 0
check_kernel(nilp, 6, 101, block_wiedemann_kernel(nilp, 6, 101, s=2))
def test_random_singular_and_nonsingular() -> None:
rng = random.Random(2002)
found = total = 0
for trial in range(500):
p = [2, 3, 101, 65521][trial % 4]
s = [1, 2, 4][trial % 3]
n = rng.randint(1, 16)
a = random_sparse(rng, n, p, rng.randint(1, 4))
if rng.random() < 0.6:
a = make_singular(rng, a, n, p)
singular = rank_mod_p(to_dense(a, n), p) < n
u = block_wiedemann_kernel(a, n, p, s=s, seed=trial, tries=3)
if not singular:
assert u is None
else:
total += 1
if u is not None:
check_kernel(a, n, p, u)
found += 1
elif p > 3:
raise AssertionError((p, s, n))
assert found >= 0.95 * total
def test_larger() -> None:
rng = random.Random(7)
for p, s in ((2, 8), (65521, 4)):
n = 120
a = make_singular(rng, random_sparse(rng, n, p, 5), n, p)
check_kernel(a, n, p, block_wiedemann_kernel(a, n, p, s=s, seed=1))
def test_example() -> None:
from block_wiedemann_algorithm import apply_block
p, n, s = 7, 4, 2
dense = [[1, 2, 0, 1], [0, 1, 3, 0], [1, 3, 3, 1], [2, 0, 0, 1]]
rows = [[(j, a) for j, a in enumerate(r) if a] for r in dense]
x = [[1, 0], [0, 1], [0, 0], [0, 0]]
z = [[0, 0], [0, 0], [1, 0], [0, 1]]
v = apply_block(rows, z, p)
assert v == [[0, 1], [3, 0], [3, 1], [0, 1]]
seq = []
for _ in range(8):
seq.append([[sum(x[r][a] * v[r][c] for r in range(n)) % p for c in range(s)] for a in range(s)])
v = apply_block(rows, v, p)
assert seq[:4] == [[[0, 1], [3, 0]], [[6, 2], [5, 3]], [[2, 4], [3, 4]], [[6, 5], [4, 0]]]
d, f = minimal_generator(seq, s, p)[0]
assert (d, f) == (1, [[6, 1], [5, 2]])
az = apply_block(rows, z, p)
u = [(z[r][0] * 6 + z[r][1] * 1 + az[r][0] * 5 + az[r][1] * 2) % p for r in range(n)]
assert u == [2, 1, 2, 3] and not any(apply_vector(rows, u, p))
What the next level needs
- the tests did not run in this build