diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..7a60b85 --- /dev/null +++ b/.gitignore @@ -0,0 +1,2 @@ +__pycache__/ +*.pyc diff --git a/README.md b/README.md index 6e3f1b5..1e66a90 100644 --- a/README.md +++ b/README.md @@ -1,2 +1,70 @@ # math0 Math explorations. + +## gfq: minimum-distance checker for linear codes over GF(q) + +`gfq/` is an open checker for the first step in [AREAS.md](AREAS.md) §1: +given a generator matrix over GF(q), compute the exact minimum distance +and print a minimum-weight codeword as a certificate that anyone can +re-check. It makes no claims about new code records. + +**Fields.** q ∈ {2, 3, 4, 5, 7, 8, 9}. Prime fields are residues mod p. +GF(4), GF(8) and GF(9) use explicit polynomial tables +(GF(2)[x]/(x²+x+1), GF(2)[x]/(x³+x+1), GF(3)[x]/(x²+1)); elements are the +integers 0..q−1 encoding polynomial coefficients in base p. There is no +field GF(6), and q=6 is rejected. Zero dependencies (Python ≥ 3.10). + +**Use.** + +``` +python3 -m gfq --q 3 examples/golay11.txt # prints [11,6,5] + certificate +python3 -m gfq --q 4 --method bz --json matrix.txt # machine-readable +cat matrix.txt | python3 -m gfq --q 3 - # stdin +pip install . # optional: provides `gfq-check` +``` + +Matrix files: one row per line, entries separated by spaces or commas +(brackets as pasted from Magma/Sage also work), compact digit strings +such as `20121` work for q ≤ 9, `#` comments are ignored, and an optional +leading `q k n` header is accepted. Exit codes: 0 exact, 2 budget +exhausted (only bounds claimed), 1 input error. The printed certificate +is checked to lie in the row space with the claimed weight before it is +reported. + +**Methods and limits.** `--method brute` enumerates all +(q^k−1)/(q−1) normalised codewords: exact, feasible only for small q^k +(`auto` uses it up to 400,000 by default). `--method bz` is a +Brouwer–Zimmermann-style enumeration: in g disjoint information sets it +enumerates message weights s = 1, 2, …, which proves +d ≥ min(upper bound, g·(s+1)) at each level, and stops when the bounds +meet (exact) or the `--max-enum` budget runs out (then it reports an +honest interval and exit code 2, never a guessed d). Cost per level is +about g·C(k,s)·(q−1)^(s−1) codewords, so it is strong when d is small +relative to k or several disjoint information sets exist, and weak for +large-k codes with large d — those need the budget raised or a stronger +implementation (bitsets/GPU, better information sets). Finite +enumerations here are labelled VERIFIED, not "proved"; expected values +in the tests are KNOWN theorems. + +**Tests.** `python3 -m unittest discover -s tests` validates the field +axioms exhaustively for every supported q, and the distance against +known codes: ternary Golay [11,6,5] and [12,6,6], Hamming codes over +GF(3)/GF(4)/GF(8)/GF(9) (including [21,18,3] and [40,36,3] cases where +brute force is infeasible), Reed–Solomon examples over +GF(5)/GF(7)/GF(8)/GF(9), and brute-vs-BZ agreement on small random codes. + +**Attacking the smallest open codetables entries next.** The census in +[data/codetables/](data/codetables/) counts 80,652 open non-binary +entries. A sane order: (1) re-run `python3 data/codetables/census.py +--csv entries.csv` and sort open entries by gap, then by n — single-gap, +small-n entries first, remembering AREAS.md's caution that large GF(7) +gaps mostly reflect weak upper bounds; (2) for a target [n,k] over GF(q), +fetch the record construction's generator matrix from codetables.de and +re-verify its d with this checker before searching (a construction that +cannot be re-checked is not a baseline); (3) search where Grassl's +quasi-cyclic index-2 sweep has not: quasi-twisted codes of index 3–4, +constacyclic shifts, and one-column extensions of record codes, using +this checker (BZ mode, budget set explicitly) to certify every candidate +before believing it; (4) only then consider submitting a verified +improvement, generator matrix included, to codes@codetables.de. No +search or submission has been done in this repository yet. diff --git a/examples/golay11.txt b/examples/golay11.txt new file mode 100644 index 0000000..b36cc76 --- /dev/null +++ b/examples/golay11.txt @@ -0,0 +1,7 @@ +# Ternary Golay code [11,6,5] over GF(3): cyclic shifts of g(x)=x^5+x^4+2x^3+x^2+2 +2 0 1 2 1 1 0 0 0 0 0 +0 2 0 1 2 1 1 0 0 0 0 +0 0 2 0 1 2 1 1 0 0 0 +0 0 0 2 0 1 2 1 1 0 0 +0 0 0 0 2 0 1 2 1 1 0 +0 0 0 0 0 2 0 1 2 1 1 diff --git a/gfq/__init__.py b/gfq/__init__.py new file mode 100644 index 0000000..c9c3e5e --- /dev/null +++ b/gfq/__init__.py @@ -0,0 +1,9 @@ +"""Open minimum-distance checker for linear codes over GF(q), q in 3..9 +(plus GF(2)). See gfq.distance for methods, gfq.cli for the CLI.""" +from .field import GF, SUPPORTED_Q +from .distance import MinDistanceResult, min_distance, brute_force_min_distance, bz_min_distance +from .matrix import parse_matrix + +__all__ = ["GF", "SUPPORTED_Q", "MinDistanceResult", "min_distance", + "brute_force_min_distance", "bz_min_distance", "parse_matrix"] +__version__ = "0.1.0" diff --git a/gfq/__main__.py b/gfq/__main__.py new file mode 100644 index 0000000..bfdcd0c --- /dev/null +++ b/gfq/__main__.py @@ -0,0 +1,4 @@ +from .cli import main + +if __name__ == "__main__": + raise SystemExit(main()) diff --git a/gfq/cli.py b/gfq/cli.py new file mode 100644 index 0000000..04b2c7d --- /dev/null +++ b/gfq/cli.py @@ -0,0 +1,85 @@ +"""Command-line minimum-distance checker. + +Usage: + python -m gfq --q 3 examples/golay11.txt + python -m gfq --q 4 --method bz --json matrix.txt + cat matrix.txt | python -m gfq --q 3 - + +Exit codes: 0 = exact distance (with certificate codeword printed); +2 = enumeration budget exhausted, only bounds are claimed; +1 = input/usage error. +""" +from __future__ import annotations + +import argparse +import json +import sys + +from .field import GF +from .matrix import parse_matrix, row_space_contains, weight +from .distance import min_distance + + +def main(argv=None) -> int: + ap = argparse.ArgumentParser( + prog="gfq-check", + description="Exact minimum distance of a linear code over GF(q) " + "from a generator matrix.") + ap.add_argument("matrix", help="generator-matrix file, or - for stdin") + ap.add_argument("--q", type=int, required=True, + help="field order: one of 2,3,4,5,7,8,9") + ap.add_argument("--method", choices=["auto", "brute", "bz"], default="auto") + ap.add_argument("--max-enum", type=int, default=2_000_000, + help="Brouwer-Zimmermann enumeration budget (default 2000000)") + ap.add_argument("--brute-limit", type=int, default=400_000, + help="auto: use brute force at or below this many codewords") + ap.add_argument("--json", action="store_true", help="machine-readable output") + args = ap.parse_args(argv) + + try: + F = GF(args.q) + text = sys.stdin.read() if args.matrix == "-" else open(args.matrix).read() + G = parse_matrix(text, F) + res = min_distance(G, F, method=args.method, + brute_limit=args.brute_limit, max_enum=args.max_enum) + except (ValueError, OSError) as e: + print(f"error: {e}", file=sys.stderr) + return 1 + + # Certificate sanity: the printed codeword really is in the row space + # and really has the claimed weight. This is a cheap independent check + # of the *upper* bound (the lower bound rests on the enumeration). + cert_ok = (row_space_contains(G, F, res.codeword) + and weight(res.codeword) == res.upper_bound) + + if args.json: + print(json.dumps({ + "q": args.q, "n": res.n, "k": res.k, + "d": res.d, "lower_bound": res.lower_bound, + "upper_bound": res.upper_bound, "complete": res.complete, + "codeword_certificate": res.codeword, + "certificate_checked": cert_ok, + "method": res.method, "enumerated": res.enumerated, + "num_info_sets": res.num_info_sets, + }, indent=2)) + else: + print(res.params()) + print(f"field: GF({args.q}) method: {res.method}" + + (f" disjoint info sets: {res.num_info_sets}" if res.method == "bz" else "") + + f" combinations enumerated: {res.enumerated}") + label = ("minimum-weight codeword" if res.complete + else "lightest codeword found (upper-bound certificate only)") + print(f"{label} (certificate, weight {res.upper_bound}): " + + " ".join(str(a) for a in res.codeword)) + print(f"certificate in row space with claimed weight: {cert_ok}") + if not res.complete: + print("status: INCOMPLETE - enumeration budget exhausted; " + "only the stated bounds are claimed, not an exact d.") + if not cert_ok: # pragma: no cover - would indicate an internal bug + print("error: internal certificate check failed", file=sys.stderr) + return 1 + return 0 if res.complete else 2 + + +if __name__ == "__main__": # pragma: no cover + raise SystemExit(main()) diff --git a/gfq/codes.py b/gfq/codes.py new file mode 100644 index 0000000..401029f --- /dev/null +++ b/gfq/codes.py @@ -0,0 +1,83 @@ +"""Classical code constructions with KNOWN parameters. + +These exist so the checker can be validated against parameters that are +theorems in the literature (Hamming: d=3; Reed-Solomon: MDS, d=n-k+1; +ternary Golay: [11,6,5] / [12,6,6]), not merely against itself. The test +suite additionally cross-checks the Brouwer-Zimmermann method against +brute force on the smaller examples, so no expected value rests on a +single code path. +""" +from __future__ import annotations + +from itertools import product + +from .field import GF +from .matrix import Matrix, nullspace_basis + + +def hamming_generator(F: GF, r: int) -> Matrix: + """Generator matrix of the q-ary Hamming code with redundancy r: + [ (q^r-1)/(q-1), (q^r-1)/(q-1) - r, 3 ] (KNOWN: d=3 for r >= 2).""" + if r < 2: + raise ValueError("Hamming codes need r >= 2") + cols = [] + for vec in product(F.elements(), repeat=r): + if all(a == 0 for a in vec): + continue + first = next(a for a in vec if a != 0) + if first == 1: # one representative per 1-dim subspace + cols.append(vec) + H = [[cols[j][i] for j in range(len(cols))] for i in range(r)] + return nullspace_basis(H, F) + + +def reed_solomon_generator(F: GF, n: int, k: int, points=None) -> Matrix: + """[n,k,n-k+1] Reed-Solomon code (KNOWN: MDS) by evaluating all + polynomials of degree < k at n distinct points of GF(q).""" + if not 1 <= k <= n <= F.q: + raise ValueError("need 1 <= k <= n <= q for this RS construction") + pts = list(points) if points is not None else F.elements()[:n] + if len(set(pts)) != n: + raise ValueError("evaluation points must be distinct") + return [[F.pow(x, i) for x in pts] for i in range(k)] + + +def repetition_generator(F: GF, n: int) -> Matrix: + """[n,1,n] repetition code.""" + return [[1] * n] + + +def single_parity_generator(F: GF, n: int) -> Matrix: + """[n,n-1,2] single-parity-check code (n >= 2).""" + if n < 2: + raise ValueError("need n >= 2") + rows = [] + for i in range(n - 1): + row = [0] * n + row[i] = 1 + row[n - 1] = F.neg(1) + rows.append(row) + return rows + + +# Ternary Golay code: cyclic [11,6,5] code over GF(3) generated by +# g(x) = x^5 + x^4 + 2x^3 + x^2 + 2, a degree-5 factor of x^11 - 1 over +# GF(3). (Its parameters are classical; in this repo the value d=5 was +# also confirmed by exhaustively enumerating all 3^6 codewords of this +# exact matrix with an independent throwaway script, and the committed +# tests re-confirm it by brute force.) Rows are the cyclic shifts +# x^i g(x), i = 0..5, coefficients low -> high. +_GOLAY_G = [2, 0, 1, 2, 1, 1] + + +def ternary_golay_generator(extended: bool = False) -> Matrix: + """[11,6,5] ternary Golay generator; with ``extended=True`` an overall + parity coordinate is appended giving the [12,6,6] extended code.""" + rows = [] + for i in range(6): + row = [0] * i + list(_GOLAY_G) + [0] * (11 - i - len(_GOLAY_G)) + if extended: + F = GF(3) + row.append(F.neg(sum(row) % 3)) + rows.append(row) + return rows diff --git a/gfq/distance.py b/gfq/distance.py new file mode 100644 index 0000000..42f8343 --- /dev/null +++ b/gfq/distance.py @@ -0,0 +1,179 @@ +"""Exact minimum distance of a linear code over GF(q). + +Two exact methods: + +``brute`` + Enumerate every codeword (up to the harmless normalisation that the + first nonzero message coefficient is 1, which does not change weights). + Cost is (q^k - 1)/(q - 1) codewords. Exact, but only feasible for + small q^k. + +``bz`` (Brouwer-Zimmermann style) + Find g pairwise disjoint information sets. In each corresponding + systematic form, enumerate all row combinations whose message weight + is <= s, for s = 1, 2, ... This gives an upper bound ``ub`` (the + lightest codeword seen). It also gives a *proved* lower bound: + any codeword c with weight < g*(s+1) has, by pigeonhole on the g + disjoint information sets, at most s nonzero coordinates in one of + them, so c would have been enumerated in that form. Hence + d >= min(ub, g*(s+1)). The loop stops when the two bounds meet, + which is an exact proof, or when an enumeration budget is exhausted, + in which case the (honest) interval [lower_bound, upper_bound] is + returned with ``complete=False`` and no claim of exactness. + +``auto`` picks brute force when its estimated cost is below +``brute_limit``, else Brouwer-Zimmermann. + +No method here is a heuristic search: a result with ``complete=True`` +is exact (the finite enumeration above). Such finite checks are +labelled VERIFIED in project reporting, not "proved" in the +human-mathematics sense, although the certificate codeword plus the +enumeration argument above let anyone re-verify them independently. +""" +from __future__ import annotations + +from dataclasses import dataclass +from itertools import combinations, product +from math import comb + +from .field import GF +from .matrix import Matrix, find_disjoint_info_sets, rref, validate_matrix, weight + + +@dataclass +class MinDistanceResult: + n: int + k: int + d: int | None # exact distance, or None if incomplete + lower_bound: int + upper_bound: int + codeword: list[int] # a minimum- (or lightest-found-) weight codeword + method: str + enumerated: int # number of (normalised) combinations examined + num_info_sets: int = 0 # g, for the BZ method + complete: bool = True + + def params(self) -> str: + if self.complete: + return f"[{self.n},{self.k},{self.d}]" + return f"[{self.n},{self.k},d in [{self.lower_bound},{self.upper_bound}]] (incomplete)" + + +def _basis(G: Matrix, F: GF) -> Matrix: + rows, _ = rref(G, F) + if not rows: + raise ValueError("generator matrix has rank 0: the code is {0} and " + "minimum distance is undefined") + return rows + + +def _combos(k: int, s: int, F: GF): + """Yield coeff vectors (tuples, length k) of weight exactly s whose + first nonzero entry is 1 (one representative per scalar multiple).""" + nz = F.nonzero() + for supp in combinations(range(k), s): + for rest in product(nz, repeat=s - 1): + coeffs = [0] * k + coeffs[supp[0]] = 1 + for pos, val in zip(supp[1:], rest): + coeffs[pos] = val + yield coeffs + + +def _combo_count(k: int, s: int, q: int) -> int: + return comb(k, s) * (q - 1) ** (s - 1) if s >= 1 else 0 + + +def _eval_combo(rows: Matrix, coeffs, F: GF, cutoff: int): + """Return (weight, codeword); weight may be reported as > cutoff early + only after the full combination is formed (cancellation forbids + earlier pruning), but counting stops once it reaches cutoff.""" + n = len(rows[0]) + out = [0] * n + add, mul = F._add, F._mul + for c, row in zip(coeffs, rows): + if c == 1: + out = [add[a][b] for a, b in zip(out, row)] + elif c: + out = [add[a][mul[c][b]] for a, b in zip(out, row)] + w = 0 + for a in out: + if a: + w += 1 + if w >= cutoff: + break + return w, out + + +def brute_force_min_distance(G: Matrix, F: GF) -> MinDistanceResult: + validate_matrix(G, F) + basis = _basis(G, F) + k, n = len(basis), len(basis[0]) + best_w, best_c, total = n + 1, None, 0 + for s in range(1, k + 1): + for coeffs in _combos(k, s, F): + total += 1 + w, c = _eval_combo(basis, coeffs, F, best_w) + if w < best_w: + best_w, best_c = w, c + if best_w == 1: + break + if best_w == 1: + break + return MinDistanceResult(n, k, best_w, best_w, best_w, best_c, + "brute", total, 0, True) + + +def bz_min_distance(G: Matrix, F: GF, max_enum: int = 2_000_000) -> MinDistanceResult: + validate_matrix(G, F) + basis = _basis(G, F) + k, n = len(basis), len(basis[0]) + + # Upper bound seed: every nonzero row of the *given* G is a codeword. + ub, best_c = n + 1, None + for row in G: + w = weight(row) + if 0 < w < ub: + ub, best_c = w, list(row) + + info_sets = find_disjoint_info_sets(basis, F) + g = len(info_sets) + enumerated = 0 + lb = 1 + s = 0 + while lb < ub: + s += 1 + if s > k: # everything enumerated in the first form: exact + lb = ub + break + level_cost = g * _combo_count(k, s, F.q) + if enumerated + level_cost > max_enum: + return MinDistanceResult(n, k, None, lb, ub, best_c, "bz", + enumerated, g, False) + for rows, _info in info_sets: + for coeffs in _combos(k, s, F): + enumerated += 1 + w, c = _eval_combo(rows, coeffs, F, ub) + if 0 < w < ub: + ub, best_c = w, c + lb = min(ub, g * (s + 1)) + return MinDistanceResult(n, k, ub, ub, ub, best_c, "bz", + enumerated, g, True) + + +def min_distance(G: Matrix, F: GF, method: str = "auto", + brute_limit: int = 400_000, + max_enum: int = 2_000_000) -> MinDistanceResult: + """Exact minimum distance. See module docstring for methods/limits.""" + if method not in ("auto", "brute", "bz"): + raise ValueError(f"unknown method {method!r}; use auto|brute|bz") + validate_matrix(G, F) + k = len(_basis(G, F)) + if method == "brute": + return brute_force_min_distance(G, F) + if method == "bz": + return bz_min_distance(G, F, max_enum=max_enum) + brute_cost = (F.q ** k - 1) // (F.q - 1) if F.q > 1 else 0 + if brute_cost <= brute_limit: + return brute_force_min_distance(G, F) + return bz_min_distance(G, F, max_enum=max_enum) diff --git a/gfq/field.py b/gfq/field.py new file mode 100644 index 0000000..69ec76b --- /dev/null +++ b/gfq/field.py @@ -0,0 +1,156 @@ +"""Finite fields GF(q) for q in {2,3,4,5,7,8,9}. + +Elements are encoded as integers 0..q-1. For a prime field GF(p) the +encoding is the residue itself. For an extension field GF(p^m) the integer +is the base-p digit string of the polynomial coefficients: +element = sum_i digit_i * x^i <-> code = sum_i digit_i * p^i. +Multiplication is polynomial multiplication modulo a fixed irreducible +polynomial (see _IRREDUCIBLE). Addition/multiplication tables are built +once at construction, which is cheap for q <= 9. + +There is no field of order 6, so q=6 (and any non-prime-power, or a prime +power other than the supported ones) raises ValueError. +""" +from __future__ import annotations + + +# Irreducible polynomials, coefficients low -> high over GF(p). +_IRREDUCIBLE = { + 4: (2, [1, 1, 1]), # x^2 + x + 1 over GF(2) + 8: (2, [1, 1, 0, 1]), # x^3 + x + 1 over GF(2) + 9: (3, [1, 0, 1]), # x^2 + 1 over GF(3) (no root mod 3) +} +_PRIMES = {2, 3, 5, 7} +SUPPORTED_Q = (2, 3, 4, 5, 7, 8, 9) + + +def _digits(code: int, p: int, m: int) -> list[int]: + out = [] + for _ in range(m): + out.append(code % p) + code //= p + return out + + +def _from_digits(digs: list[int], p: int) -> int: + code = 0 + for d in reversed(digs): + code = code * p + (d % p) + return code + + +class GF: + """The finite field with ``q`` elements; elements are ints 0..q-1.""" + + def __init__(self, q: int): + if q in _PRIMES: + self.q, self.p, self.m = q, q, 1 + self.modulus = None + elif q in _IRREDUCIBLE: + self.q = q + self.p, self.modulus = _IRREDUCIBLE[q] + self.m = len(self.modulus) - 1 + else: + raise ValueError( + f"no supported field of order {q}; supported q: {SUPPORTED_Q} " + "(q must be a prime power; there is no field GF(6))" + ) + self.zero = 0 + self.one = 1 + self._add = [[self._add_raw(a, b) for b in range(q)] for a in range(q)] + self._mul = [[self._mul_raw(a, b) for b in range(q)] for a in range(q)] + self._neg = [self._find_neg(a) for a in range(q)] + self._inv = [0] + [self._find_inv(a) for a in range(1, q)] + + # -- raw arithmetic (used only to build tables) --------------------- + def _add_raw(self, a: int, b: int) -> int: + if self.m == 1: + return (a + b) % self.p + da, db = _digits(a, self.p, self.m), _digits(b, self.p, self.m) + return _from_digits([(x + y) % self.p for x, y in zip(da, db)], self.p) + + def _mul_raw(self, a: int, b: int) -> int: + if self.m == 1: + return (a * b) % self.p + p, mod = self.p, self.modulus + da, db = _digits(a, p, self.m), _digits(b, p, self.m) + prod = [0] * (2 * self.m - 1) + for i, x in enumerate(da): + for j, y in enumerate(db): + prod[i + j] = (prod[i + j] + x * y) % p + # reduce modulo the (monic) irreducible polynomial + for deg in range(len(prod) - 1, self.m - 1, -1): + c = prod[deg] % p + if c: + shift = deg - self.m + for j in range(len(mod)): + prod[shift + j] = (prod[shift + j] - c * mod[j]) % p + return _from_digits(prod[: self.m], p) + + def _find_neg(self, a: int) -> int: + for b in range(self.q): + if self._add[a][b] == 0: + return b + raise AssertionError("no additive inverse") # pragma: no cover + + def _find_inv(self, a: int) -> int: + for b in range(self.q): + if self._mul[a][b] == 1: + return b + raise AssertionError("no multiplicative inverse") # pragma: no cover + + # -- public arithmetic ------------------------------------------------ + def _check(self, a: int) -> int: + if not isinstance(a, int) or not 0 <= a < self.q: + raise ValueError(f"{a!r} is not an element of GF({self.q}) (0..{self.q - 1})") + return a + + def add(self, a: int, b: int) -> int: + return self._add[self._check(a)][self._check(b)] + + def sub(self, a: int, b: int) -> int: + return self._add[self._check(a)][self._neg[self._check(b)]] + + def neg(self, a: int) -> int: + return self._neg[self._check(a)] + + def mul(self, a: int, b: int) -> int: + return self._mul[self._check(a)][self._check(b)] + + def inv(self, a: int) -> int: + self._check(a) + if a == 0: + raise ZeroDivisionError("0 has no inverse") + return self._inv[a] + + def div(self, a: int, b: int) -> int: + return self.mul(a, self.inv(b)) + + def pow(self, a: int, e: int) -> int: + self._check(a) + if e < 0: + return self.pow(self.inv(a), -e) + r = 1 + for _ in range(e): + r = self._mul[r][a] + return r + + def elements(self) -> list[int]: + return list(range(self.q)) + + def nonzero(self) -> list[int]: + return list(range(1, self.q)) + + def primitive_element(self) -> int: + """Smallest generator of the multiplicative group.""" + for a in range(1, self.q): + seen, cur = set(), 1 + for _ in range(self.q - 1): + seen.add(cur) + cur = self._mul[cur][a] + if len(seen) == self.q - 1: + return a + raise AssertionError("multiplicative group is cyclic") # pragma: no cover + + def __repr__(self) -> str: # pragma: no cover + return f"GF({self.q})" diff --git a/gfq/matrix.py b/gfq/matrix.py new file mode 100644 index 0000000..80f6f84 --- /dev/null +++ b/gfq/matrix.py @@ -0,0 +1,190 @@ +"""Matrices over GF(q): parsing, row reduction, systematic forms. + +A matrix is a list of rows; each row is a list of field elements (ints). +""" +from __future__ import annotations + +import re + +from .field import GF + +Matrix = list[list[int]] + + +def validate_matrix(G: Matrix, F: GF) -> tuple[int, int]: + """Return (rows, cols); raise ValueError on ragged/out-of-field input.""" + if not G or not G[0]: + raise ValueError("matrix must have at least one row and one column") + n = len(G[0]) + for row in G: + if len(row) != n: + raise ValueError("matrix rows have different lengths") + for a in row: + F._check(a) + return len(G), n + + +def parse_matrix(text: str, F: GF) -> Matrix: + """Parse a generator matrix from text. + + Accepted layouts (auto-detected per line): + * whitespace- or comma-separated integers: ``1 0 2`` / ``1,0,2`` + * bracketed rows as pasted from Magma/Sage/Python: ``[1,0,2]`` + * compact digit strings (only unambiguous because q <= 9): ``102`` + Blank lines and ``#`` comments are ignored. A leading header line of + the form ``q k n`` or ``k n`` (three/two integers on their own line, + where the remaining row count matches) is also accepted. + """ + lines = [] + for raw in text.splitlines(): + line = raw.split("#", 1)[0].strip().rstrip(";") + if line: + lines.append(line) + if not lines: + raise ValueError("no matrix data found") + + rows: Matrix = [] + for idx, line in enumerate(lines): + s = line.strip().strip("[]") + if re.fullmatch(r"[0-9]+", s) and len(s) > 1 and (rows or len(lines) > 1): + # compact row, unless it is a lone header candidate handled below + vals = [int(ch) for ch in s] + # Heuristic: a first line like "3 6 11" never reaches here (spaces), + # and a single compact line is one row. Keep it as a row. + else: + toks = [t for t in re.split(r"[\s,\[\]]+", s) if t != ""] + if not toks: + continue + if idx == 0 and len(toks) in (2, 3) and len(lines) > 1: + # Possible header "k n" / "q k n": only treat as header when + # the numbers are inconsistent with being a matrix row of the + # width determined by the second line. + second = lines[1].strip().strip("[]") + if re.fullmatch(r"[0-9]+", second): + width2 = len(second) + else: + width2 = len([t for t in re.split(r"[\s,\[\]]+", second) if t]) + hdr = [int(t) for t in toks] + # Header only when the values read as (q,k,n)/(k,n) + # consistent with the remaining rows; otherwise the line + # is a matrix row (and raggedness is reported as such). + consistent = ( + (len(hdr) == 3 and hdr[1] == len(lines) - 1 and hdr[2] == width2) + or (len(hdr) == 2 and hdr[0] == len(lines) - 1 and hdr[1] == width2)) + if consistent: + if len(hdr) == 3 and hdr[0] != F.q: + raise ValueError( + f"header says q={hdr[0]} but --q is {F.q}") + continue + try: + vals = [int(t) for t in toks] + except ValueError as e: + raise ValueError(f"cannot parse matrix line {idx + 1}: {raw!r}") from e + rows.append(vals) + validate_matrix(rows, F) + return rows + + +def rref(G: Matrix, F: GF, col_order: list[int] | None = None): + """Reduced row echelon form. + + Columns are considered as pivot candidates in ``col_order`` (default + left-to-right). Returns (rref_rows, pivot_cols) with zero rows dropped; + row i has its pivot in column pivot_cols[i]. The returned matrix is in + the original column order (not permuted). + """ + rows = [list(r) for r in G] + ncols = len(rows[0]) if rows else 0 + order = list(range(ncols)) if col_order is None else list(col_order) + pivots: list[int] = [] + r = 0 + for c in order: + if r >= len(rows): + break + piv = next((i for i in range(r, len(rows)) if rows[i][c] != 0), None) + if piv is None: + continue + rows[r], rows[piv] = rows[piv], rows[r] + inv = F.inv(rows[r][c]) + rows[r] = [F.mul(inv, a) for a in rows[r]] + for i in range(len(rows)): + if i != r and rows[i][c] != 0: + f = rows[i][c] + rows[i] = [F.sub(a, F.mul(f, b)) for a, b in zip(rows[i], rows[r])] + pivots.append(c) + r += 1 + return rows[:r], pivots + + +def rank(G: Matrix, F: GF) -> int: + return len(rref(G, F)[0]) + + +def nullspace_basis(H: Matrix, F: GF) -> Matrix: + """Basis (as rows) of {x : H x^T = 0}; used to turn a parity-check + matrix into a generator matrix.""" + R, pivots = rref(H, F) + n = len(H[0]) + free = [c for c in range(n) if c not in set(pivots)] + basis = [] + for f in free: + x = [0] * n + x[f] = 1 + for i, pc in enumerate(pivots): + x[pc] = F.neg(R[i][f]) + basis.append(x) + return basis + + +def row_space_contains(G: Matrix, F: GF, vec: list[int]) -> bool: + """True iff ``vec`` is a GF(q)-linear combination of the rows of G.""" + R, _ = rref(G, F) + return rank(R + [list(vec)], F) == len(R) + + +def systematic_form(G: Matrix, F: GF, col_order: list[int]): + """Return (rows, info_cols) where rows form [I | *] on ``info_cols``: + row i is 1 in column info_cols[i] and 0 in the other info columns. + ``col_order`` lists preferred information columns first.""" + R, pivots = rref(G, F, col_order=col_order) + return R, pivots + + +def find_disjoint_info_sets(G: Matrix, F: GF, max_sets: int | None = None): + """Greedily find pairwise disjoint information sets (pivot column sets). + + Returns a list of (systematic_rows, info_cols). Each successive set is + searched with the already-used columns deprioritised; a set counts only + if its pivots avoid every column used by earlier sets. At least one set + is always returned for a nonzero code. Finding the maximum number of + disjoint information sets is not attempted (greedy is enough for the + Brouwer-Zimmermann lower bound; fewer sets only slow the proof). + """ + n = len(G[0]) + k = rank(G, F) + used: set[int] = set() + out = [] + cap = max_sets if max_sets is not None else max(1, n // max(k, 1)) + while len(out) < cap: + order = [c for c in range(n) if c not in used] + sorted(used) + rows, pivots = systematic_form(G, F, order) + if len(pivots) < k: + break + if out and any(c in used for c in pivots): + break + out.append((rows, pivots)) + used.update(pivots) + return out + + +def weight(vec: list[int]) -> int: + return sum(1 for a in vec if a != 0) + + +def combine(rows: Matrix, coeffs: list[int], F: GF) -> list[int]: + n = len(rows[0]) + out = [0] * n + for c, row in zip(coeffs, rows): + if c: + out = [F.add(a, F.mul(c, b)) for a, b in zip(out, row)] + return out diff --git a/pyproject.toml b/pyproject.toml new file mode 100644 index 0000000..f80fc30 --- /dev/null +++ b/pyproject.toml @@ -0,0 +1,16 @@ +[build-system] +requires = ["setuptools>=61"] +build-backend = "setuptools.build_meta" + +[project] +name = "math0-gfq" +version = "0.1.0" +description = "Exact minimum-distance checker for linear codes over GF(q), q in 2..9 (prime powers)" +requires-python = ">=3.10" +dependencies = [] + +[project.scripts] +gfq-check = "gfq.cli:main" + +[tool.setuptools] +packages = ["gfq"] diff --git a/tests/test_cli.py b/tests/test_cli.py new file mode 100644 index 0000000..a1ab491 --- /dev/null +++ b/tests/test_cli.py @@ -0,0 +1,65 @@ +import io +import json +import os +import unittest +from contextlib import redirect_stdout, redirect_stderr + +from gfq.cli import main + +ROOT = os.path.join(os.path.dirname(__file__), "..") + + +class TestCLI(unittest.TestCase): + def run_cli(self, argv): + out, err = io.StringIO(), io.StringIO() + with redirect_stdout(out), redirect_stderr(err): + rc = main(argv) + return rc, out.getvalue(), err.getvalue() + + def test_golay_example(self): + rc, out, _ = self.run_cli( + ["--q", "3", os.path.join(ROOT, "examples", "golay11.txt")]) + self.assertEqual(rc, 0) + self.assertIn("[11,6,5]", out) + self.assertIn("True", out) # certificate check passed + + def test_golay_json(self): + rc, out, _ = self.run_cli( + ["--q", "3", "--json", "--method", "bz", + os.path.join(ROOT, "examples", "golay11.txt")]) + self.assertEqual(rc, 0) + data = json.loads(out) + self.assertEqual((data["n"], data["k"], data["d"]), (11, 6, 5)) + self.assertTrue(data["certificate_checked"]) + self.assertEqual(sum(1 for a in data["codeword_certificate"] if a), 5) + + def test_wrong_field_element_error(self): + import tempfile + with tempfile.NamedTemporaryFile("w", suffix=".txt", delete=False) as f: + f.write("1 4\n0 1\n") + path = f.name + try: + rc, _, err = self.run_cli(["--q", "3", path]) + finally: + os.unlink(path) + self.assertEqual(rc, 1) + self.assertIn("error", err) + + def test_incomplete_exit_code_2(self): + import tempfile + from gfq import GF, codes + G = codes.hamming_generator(GF(4), 3) + with tempfile.NamedTemporaryFile("w", suffix=".txt", delete=False) as f: + f.write("\n".join(" ".join(map(str, r)) for r in G)) + path = f.name + try: + rc, out, _ = self.run_cli( + ["--q", "4", "--method", "bz", "--max-enum", "1", path]) + finally: + os.unlink(path) + self.assertEqual(rc, 2) + self.assertIn("INCOMPLETE", out) + + +if __name__ == "__main__": + unittest.main() diff --git a/tests/test_distance.py b/tests/test_distance.py new file mode 100644 index 0000000..560e971 --- /dev/null +++ b/tests/test_distance.py @@ -0,0 +1,143 @@ +"""Validation against codes whose parameters are KNOWN theorems. + +Every expected [n,k,d] below is classical (Hamming d=3, Reed-Solomon MDS +d=n-k+1, ternary Golay [11,6,5]/[12,6,6], repetition, single-parity). +Where size permits, brute force and Brouwer-Zimmermann must agree, so no +expected value rests on one code path alone. +""" +import unittest + +from gfq import GF, min_distance, brute_force_min_distance, bz_min_distance +from gfq import codes +from gfq.matrix import parse_matrix, row_space_contains, weight + + +def check(tc, G, q, n, k, d, methods=("auto", "brute", "bz")): + F = GF(q) + tc.assertEqual(len(G[0]), n) + for m in methods: + res = min_distance(G, F, method=m) + tc.assertTrue(res.complete, f"{m}: incomplete") + tc.assertEqual((res.n, res.k, res.d), (n, k, d), m) + tc.assertEqual(weight(res.codeword), d, m) + tc.assertTrue(row_space_contains(G, F, res.codeword), m) + return res + + +class TestKnownCodes(unittest.TestCase): + def test_ternary_golay(self): + check(self, codes.ternary_golay_generator(False), 3, 11, 6, 5) + + def test_ternary_golay_extended(self): + check(self, codes.ternary_golay_generator(True), 3, 12, 6, 6) + + def test_hamming_gf3(self): + check(self, codes.hamming_generator(GF(3), 2), 3, 4, 2, 3) + + def test_hamming_gf3_r3(self): + # [13,10,3]: 3^10/2 normalised codewords, brute still feasible. + check(self, codes.hamming_generator(GF(3), 3), 3, 13, 10, 3) + + def test_hamming_gf4(self): + check(self, codes.hamming_generator(GF(4), 2), 4, 5, 3, 3) + + def test_hamming_gf4_r3_bz_only(self): + # [21,18,3]: brute force (4^18) infeasible; BZ must close exactly. + F = GF(4) + res = bz_min_distance(codes.hamming_generator(F, 3), F) + self.assertTrue(res.complete) + self.assertEqual((res.n, res.k, res.d), (21, 18, 3)) + + def test_hamming_gf8_gf9(self): + # n=(q^2-1)/(q-1)=q+1: GF(8) -> [9,7,3], GF(9) -> [10,8,3] + check(self, codes.hamming_generator(GF(8), 2), 8, 9, 7, 3, + methods=("auto", "bz")) + check(self, codes.hamming_generator(GF(9), 2), 9, 10, 8, 3, + methods=("auto", "bz")) + + def test_reed_solomon(self): + check(self, codes.reed_solomon_generator(GF(8), 7, 3), 8, 7, 3, 5) + check(self, codes.reed_solomon_generator(GF(8), 7, 4), 8, 7, 4, 4, + methods=("auto", "bz")) + check(self, codes.reed_solomon_generator(GF(9), 9, 4), 9, 9, 4, 6, + methods=("auto", "bz")) + check(self, codes.reed_solomon_generator(GF(7), 7, 3), 7, 7, 3, 5) + check(self, codes.reed_solomon_generator(GF(5), 5, 2), 5, 5, 2, 4) + check(self, codes.reed_solomon_generator(GF(3), 3, 2), 3, 3, 2, 2) + + def test_trivial_codes(self): + for q in (3, 4, 9): + F = GF(q) + check(self, codes.repetition_generator(F, 6), q, 6, 1, 6) + check(self, codes.single_parity_generator(F, 5), q, 5, 4, 2) + ident = [[1 if i == j else 0 for j in range(4)] for i in range(4)] + check(self, ident, q, 4, 4, 1) + + def test_binary_sanity(self): + # binary Hamming [7,4,3] and Golay-free sanity that GF(2) also works + check(self, codes.hamming_generator(GF(2), 3), 2, 7, 4, 3) + + +class TestBehaviour(unittest.TestCase): + def test_rank_deficient_generator(self): + F = GF(3) + G = [[1, 0, 1], [2, 0, 2], [0, 1, 1]] # rows 1,2 proportional: k=2 + res = min_distance(G, F) + self.assertEqual(res.k, 2) + self.assertEqual(res.d, 2) + + def test_zero_code_rejected(self): + with self.assertRaises(ValueError): + min_distance([[0, 0, 0]], GF(3)) + + def test_out_of_field_rejected(self): + with self.assertRaises(ValueError): + min_distance([[1, 3]], GF(3)) + + def test_bz_incomplete_is_honest(self): + # A starved budget must yield bounds, never a claimed exact d. + F = GF(4) + res = bz_min_distance(codes.hamming_generator(F, 3), F, max_enum=1) + self.assertFalse(res.complete) + self.assertIsNone(res.d) + self.assertLessEqual(res.lower_bound, 3) + self.assertGreaterEqual(res.upper_bound, 3) + + def test_brute_matches_bz_random_small(self): + # Deterministic pseudo-random small codes (fixed LCG, no `random` + # import needed): independent methods must agree exactly. + seed = 12345 + def rnd(mod): + nonlocal seed + seed = (seed * 1103515245 + 12345) & 0x7FFFFFFF + return seed % mod + for q, k, n in ((3, 4, 9), (4, 3, 8), (5, 3, 7), (9, 3, 8), (8, 3, 8)): + F = GF(q) + for _ in range(3): + G = [[rnd(q) for _ in range(n)] for _ in range(k)] + try: + a = brute_force_min_distance(G, F) + except ValueError: + continue + b = bz_min_distance(G, F) + self.assertTrue(b.complete) + self.assertEqual(a.d, b.d) + self.assertEqual(a.k, b.k) + + def test_parse_formats(self): + F = GF(3) + a = parse_matrix("1 0 2\n0 1 1\n", F) + b = parse_matrix("[1,0,2]\n[0,1,1]\n", F) + c = parse_matrix("102\n011\n", F) + d_ = parse_matrix("# hdr\n3 2 3\n1 0 2\n0 1 1\n", F) + self.assertEqual(a, b) + self.assertEqual(a, c) + self.assertEqual(a, d_) + with self.assertRaises(ValueError): + parse_matrix("1 0\n1 0 0\n", F) # ragged + with self.assertRaises(ValueError): + parse_matrix("1 5\n0 1\n", F) # 5 not in GF(3) + + +if __name__ == "__main__": + unittest.main() diff --git a/tests/test_field.py b/tests/test_field.py new file mode 100644 index 0000000..7255e8c --- /dev/null +++ b/tests/test_field.py @@ -0,0 +1,51 @@ +import unittest +from gfq import GF + + +class TestField(unittest.TestCase): + def test_supported_and_unsupported(self): + for q in (2, 3, 4, 5, 7, 8, 9): + self.assertEqual(GF(q).q, q) + for q in (1, 6, 10, 11, 16): + with self.assertRaises(ValueError): + GF(q) + + def test_field_axioms_exhaustive(self): + # VERIFIED by finite enumeration for every supported field. + for q in (2, 3, 4, 5, 7, 8, 9): + F = GF(q) + for a in range(q): + self.assertEqual(F.add(a, 0), a) + self.assertEqual(F.mul(a, 1), a) + self.assertEqual(F.add(a, F.neg(a)), 0) + if a: + self.assertEqual(F.mul(a, F.inv(a)), 1) + for b in range(q): + self.assertEqual(F.add(a, b), F.add(b, a)) + self.assertEqual(F.mul(a, b), F.mul(b, a)) + for c in range(q): + self.assertEqual(F.mul(a, F.add(b, c)), + F.add(F.mul(a, b), F.mul(a, c))) + # multiplicative group cyclic of order q-1 + g = F.primitive_element() + self.assertEqual({F.pow(g, e) for e in range(q - 1)}, set(range(1, q))) + + def test_prime_fields_are_residues(self): + F = GF(5) + self.assertEqual(F.mul(3, 4), 2) + self.assertEqual(F.inv(2), 3) + self.assertEqual(GF(3).add(2, 2), 1) + + def test_gf4_table_spot_check(self): + # GF(4)=GF(2)[x]/(x^2+x+1): element 2 is x, x*x = x+1 = 3. + F = GF(4) + self.assertEqual(F.mul(2, 2), 3) + self.assertEqual(F.mul(2, 3), 1) # x(x+1)=1 mod (x^2+x+1) + with self.assertRaises(ZeroDivisionError): + F.inv(0) + with self.assertRaises(ValueError): + F.add(4, 1) + + +if __name__ == "__main__": + unittest.main()