matrix_vass_review.py

#!/usr/bin/env python3
"""Independent finite checks for the matrix-power and reversible-rewriting review.

The matrix checks are exact integer computations.  The rewriting check uses a
small dependency-free Buchberger implementation over Q, independently of the
reachability search it is compared with.
"""

from __future__ import annotations

from collections import deque
from fractions import Fraction
from itertools import combinations, permutations, product
import json
import math


def matmul(a, b):
    n = len(a)
    return tuple(
        tuple(sum(a[i][k] * b[k][j] for k in range(n)) for j in range(n))
        for i in range(n)
    )


def eye(n):
    return tuple(tuple(int(i == j) for j in range(n)) for i in range(n))


def nonnegative(a):
    return all(x >= 0 for row in a for x in row)


def signed_permutation_matrix(perm, signs):
    n = len(perm)
    return tuple(
        tuple(signs[i] if j == perm[i] else 0 for j in range(n))
        for i in range(n)
    )


def signed_permutation_audit():
    checked = 0
    positive_hits = 0
    maximum_period = 0
    for n in range(1, 6):
        ident = eye(n)
        for perm in permutations(range(n)):
            for signs in product((-1, 1), repeat=n):
                m = signed_permutation_matrix(perm, signs)
                powers = [ident]
                p = ident
                while True:
                    p = matmul(p, m)
                    if p == ident:
                        break
                    powers.append(p)
                    assert len(powers) <= 2 * math.lcm(*range(1, n + 1))
                period = len(powers)
                good = {k for k, value in enumerate(powers) if nonnegative(value)}
                # Multiples of a good positive exponent stay good.
                for k in good - {0}:
                    assert all((r * k) % period in good for r in range(1, period + 1))
                checked += 1
                positive_hits += bool(good - {0})
                maximum_period = max(maximum_period, period)
    return {
        "signed_permutation_matrices": checked,
        "with_positive_nonnegative_power": positive_hits,
        "maximum_group_period": maximum_period,
    }


def signed_cycle_first_hit(p):
    perm = tuple((i + 1) % p for i in range(p))
    signs = (-1,) + (1,) * (p - 1)
    m = signed_permutation_matrix(perm, signs)
    power = eye(p)
    for n in range(1, 2 * p + 1):
        power = matmul(power, m)
        if nonnegative(power):
            return n
    raise AssertionError("signed cycle did not attain its predicted hit")


def matrix_lrs_audit():
    matrices = 0
    entry_recurrences = 0
    interleaved_values = 0
    quantifier_separations = 0
    for entries in product((-1, 0, 1), repeat=4):
        m = (entries[:2], entries[2:])
        powers = [eye(2)]
        for _ in range(12):
            powers.append(matmul(powers[-1], m))
        trace = m[0][0] + m[1][1]
        determinant = m[0][0] * m[1][1] - m[0][1] * m[1][0]
        for i, j in product(range(2), repeat=2):
            for n in range(11):
                assert powers[n + 2][i][j] == (
                    trace * powers[n + 1][i][j] - determinant * powers[n][i][j]
                )
                entry_recurrences += 1
        entries_by_index = [(i, j) for i in range(2) for j in range(2)]
        interleaved = [
            powers[n][i][j]
            for n in range(13)
            for i, j in entries_by_index
        ]
        for n in range(13):
            for r, (i, j) in enumerate(entries_by_index):
                assert interleaved[4 * n + r] == powers[n][i][j]
                interleaved_values += 1
        good = [nonnegative(p) for p in powers]
        if any(good[1:]) and not all(good[7:]):
            quantifier_separations += 1
        matrices += 1
    assert signed_cycle_first_hit(3) == 6
    assert signed_cycle_first_hit(5) == 10
    assert signed_cycle_first_hit(7) == 14
    assert signed_cycle_first_hit(11) == 22
    assert not nonnegative(((-1,),)) and nonnegative(((1,),))
    return {
        "small_integer_matrices": matrices,
        "cayley_hamilton_entry_checks": entry_recurrences,
        "interleaved_value_checks": interleaved_values,
        "bounded_quantifier_separations": quantifier_separations,
        "signed_cycle_primes": [3, 5, 7, 11],
    }


def monomials(nvars, degree):
    if nvars == 1:
        return [(degree,)]
    answer = []
    for first in range(degree + 1):
        for rest in monomials(nvars - 1, degree - first):
            answer.append((first,) + rest)
    return answer


def divides(a, b):
    return all(x <= y for x, y in zip(a, b))


def sub(a, b):
    return tuple(x - y for x, y in zip(a, b))


def add(a, b):
    return tuple(x + y for x, y in zip(a, b))


def reachable(source, target, rules):
    todo = deque([source])
    seen = {source}
    while todo:
        current = todo.popleft()
        if current == target:
            return True
        for left, right in rules:
            for before, after in ((left, right), (right, left)):
                if divides(before, current):
                    nxt = add(sub(current, before), after)
                    if nxt not in seen:
                        seen.add(nxt)
                        todo.append(nxt)
    return False


def poly_clean(p):
    return {m: c for m, c in p.items() if c}


def poly_add(p, q, factor=Fraction(1)):
    out = dict(p)
    for m, c in q.items():
        out[m] = out.get(m, Fraction(0)) + factor * c
    return poly_clean(out)


def poly_mul_monomial(p, monomial, coefficient=Fraction(1)):
    return {add(m, monomial): coefficient * c for m, c in p.items()}


def leading(p):
    m = max(p)
    return m, p[m]


def reduce_poly(p, basis):
    p = poly_clean(p)
    remainder = {}
    while p:
        lm, lc = leading(p)
        divisor = None
        for g in basis:
            gm, gc = leading(g)
            if divides(gm, lm):
                divisor = (g, gm, gc)
                break
        if divisor is None:
            remainder[lm] = remainder.get(lm, Fraction(0)) + lc
            del p[lm]
        else:
            g, gm, gc = divisor
            p = poly_add(
                p,
                poly_mul_monomial(g, sub(lm, gm), lc / gc),
                factor=Fraction(-1),
            )
    return poly_clean(remainder)


def s_polynomial(f, g):
    fm, fc = leading(f)
    gm, gc = leading(g)
    lcm = tuple(max(x, y) for x, y in zip(fm, gm))
    first = poly_mul_monomial(f, sub(lcm, fm), Fraction(1, 1) / fc)
    second = poly_mul_monomial(g, sub(lcm, gm), Fraction(1, 1) / gc)
    return poly_add(first, second, factor=Fraction(-1))


def groebner(generators):
    basis = [poly_clean(g) for g in generators if poly_clean(g)]
    pairs = deque(combinations(range(len(basis)), 2))
    while pairs:
        i, j = pairs.popleft()
        remainder = reduce_poly(s_polynomial(basis[i], basis[j]), basis)
        if remainder:
            new_index = len(basis)
            basis.append(remainder)
            for old in range(new_index):
                pairs.append((old, new_index))
    return basis


def binomial(left, right):
    if left == right:
        return {}
    return poly_clean({left: Fraction(-1), right: Fraction(1)})


def ideal_contains_difference(source, target, basis):
    return not reduce_poly(binomial(source, target), basis)


def reversible_rewriting_audit():
    systems = 0
    queries = 0
    nvars = 2
    rule_monomials = monomials(nvars, 2)
    rule_universe = list(combinations(rule_monomials, 2))
    for mask in range(1 << len(rule_universe)):
        rules = [rule_universe[i] for i in range(len(rule_universe)) if mask & (1 << i)]
        basis = groebner([binomial(a, b) for a, b in rules])
        systems += 1
        for degree in range(1, 6):
            states = monomials(nvars, degree)
            for source in states:
                for target in states:
                    assert reachable(source, target, rules) == ideal_contains_difference(
                        source, target, basis
                    )
                    queries += 1

    # Three-variable hostile extension with a fixed six-edge rule universe.
    nvars = 3
    degree_two = monomials(nvars, 2)
    selected_edges = [
        (degree_two[0], degree_two[1]),
        (degree_two[1], degree_two[2]),
        (degree_two[0], degree_two[3]),
        (degree_two[2], degree_two[4]),
        (degree_two[3], degree_two[5]),
        (degree_two[4], degree_two[5]),
    ]
    for mask in range(1 << len(selected_edges)):
        rules = [selected_edges[i] for i in range(len(selected_edges)) if mask & (1 << i)]
        basis = groebner([binomial(a, b) for a, b in rules])
        systems += 1
        states = monomials(nvars, 3)
        for source in states:
            for target in states:
                assert reachable(source, target, rules) == ideal_contains_difference(
                    source, target, basis
                )
                queries += 1

    # Ideals forget orientation: y-x is in <y-x>, but y cannot reach x
    # in the one-way system x -> y.
    x, y = (1, 0), (0, 1)
    directed_reverse_reachable = False
    directed_ideal_membership = ideal_contains_difference(y, x, groebner([binomial(x, y)]))
    assert directed_ideal_membership and not directed_reverse_reachable
    return {
        "degree_preserving_reversible_systems": systems,
        "reachability_ideal_comparisons": queries,
        "directed_hostile_example": {
            "reverse_reachable": directed_reverse_reachable,
            "ideal_membership": directed_ideal_membership,
        },
    }


def main():
    output = {
        "matrix_lrs": matrix_lrs_audit(),
        "signed_permutations": signed_permutation_audit(),
        "reversible_rewriting": reversible_rewriting_audit(),
        "status": "pass",
    }
    print(json.dumps(output, indent=2, sort_keys=True))


if __name__ == "__main__":
    main()