rectangle-consistency-check.py

#!/usr/bin/env python3
"""Exact finite corroboration; universal bounds use the accompanying proofs."""
import json
import random
from fractions import Fraction
from itertools import product
from math import comb


def code(generators):
    result = {0}
    for g in generators:
        result |= {b ^ g for b in result}
    return sorted(result)


def cell(b, n, i, j):
    return (b >> (i * n + j)) & 1


def rectangle(u, v, n):
    return sum(1 << (i * n + j) for i in range(n) for j in range(n)
               if (u >> i) & 1 and (v >> j) & 1)


def score(b, q):
    return 2 * (b & q).bit_count() - b.bit_count()


def consistent_basis(n):
    rows = [sum(1 << (i * n + j) for j in range(n)) for i in range(n)]
    cols = [sum(1 << (i * n + j) for i in range(n)) for j in range(n)]
    return rows + cols[:-1]


def predicted_l(u, v, n):
    U, V = u.bit_count(), v.bit_count()
    return max(0, U * (2 * V - n), V * (2 * U - n), n * (U + V - n))


def check(b, n, i, k, j, ell):
    return cell(b, n, i, j) ^ cell(b, n, i, ell) ^ cell(b, n, k, j) ^ cell(b, n, k, ell)


def repair(b, n, k, ell):
    return sum((cell(b, n, i, ell) ^ cell(b, n, k, j) ^ cell(b, n, k, ell))
               << (i * n + j) for i in range(n) for j in range(n))


l_results, k_results = [], []
for n in range(1, 5):
    L = code(consistent_basis(n))
    total = 0
    for u, v in product(range(1 << n), repeat=2):
        q = rectangle(u, v, n)
        actual = max(score(b, q) for b in L)
        assert actual == predicted_l(u, v, n)
        total += actual
    mu = Fraction(n * comb(n - 1, n // 2), 1 << n)
    mean = Fraction(total, 1 << (2 * n))
    assert mean == n * mu - mu * mu
    l_results.append({"n": n, "dimension": 2 * n - 1, "F_L": str(mean)})

for n in range(1, 7):
    cut_generators = [sum(((int(i == k) ^ int(j == k))) << (i * n + j)
                         for i in range(n) for j in range(n)) for k in range(n - 1)]
    K = code(cut_generators)
    total = 0
    for u, v in product(range(1 << n), repeat=2):
        q = rectangle(u, v, n)
        actual = max(score(b, q) for b in K)
        t = max(0, u.bit_count() + v.bit_count() - n)
        assert actual == 2 * ((t * t) // 4)
        total += actual
    mean = Fraction(total, 1 << (2 * n))
    assert mean == Fraction(n - 1, 8)
    k_results.append({"n": n, "dimension": n - 1, "F_K": str(mean)})

rng = random.Random(20261003)
repair_results = []
for n in (3, 4):
    rows = consistent_basis(n)[:n]
    candidates = [("row", rows), ("row_with_column_defect", rows[:-1] + [1])]
    candidates += [("generic", [rng.getrandbits(n * n) for _ in range(n)])]
    if n % 2 == 0:
        candidates.append(("half_square", [rectangle((1 << (n // 2)) - 1,
                                                     (1 << (n // 2)) - 1, n)]))
    for label, gs in candidates:
        C = code(gs)
        anchors = []
        nonzero_checks = failures = 0
        for k, ell in product(range(n), repeat=2):
            H = sum(1 << (i * n + j) for i, j in product(range(n), repeat=2)
                    if any(check(b, n, i, k, j, ell) for b in C))
            anchors.append((H.bit_count(), k, ell, H))
            nonzero_checks += H.bit_count()
            for i, j in product(range(n), repeat=2):
                failures += sum(check(b, n, i, k, j, ell) for b in C)
        rho = Fraction(failures, len(C) * n ** 4)
        assert Fraction(nonzero_checks, n ** 4) == 2 * rho
        if label == "half_square":
            assert rho == Fraction(1, 8)
        h, k, ell, H = min(anchors)
        assert h <= 2 * rho * n * n
        for b in C:
            tb = repair(b, n, k, ell)
            assert (b ^ tb) & ~H == 0
            assert all(check(tb, n, i, kk, j, ll) == 0
                       for i, kk, j, ll in product(range(n), repeat=4))
            for d in C:
                assert repair(b ^ d, n, k, ell) == tb ^ repair(d, n, k, ell)
        total = 0
        for u, v in product(range(1 << n), repeat=2):
            q = rectangle(u, v, n)
            actual = max(score(b, q) for b in C)
            assert actual <= predicted_l(u, v, n) + h
            total += actual
        if label == "row_with_column_defect":
            for kk in range(n):
                hh = next(item[0] for item in anchors if item[1:3] == (kk, 1))
                assert hh <= n - 1
        repair_results.append({"n": n, "family": label,
                               "dimension": len(C).bit_length() - 1,
                               "rho": str(rho), "common_h": h,
                               "mean_score": str(Fraction(total, 1 << (2 * n)))})

print(json.dumps({"status": "passed", "scope": "finite corroboration, not formal verification",
                  "full_consistent_spaces": l_results, "cut_codes": k_results,
                  "common_linear_repairs": repair_results}, indent=2))