verify.py

1from itertools import product23# KRONECKER45CODES = tuple(range(1, 16))6LIB3 = (3, 5, 6, 7, 9, 10, 11, 12, 13, 14)7K2 = (3, 5, 6, 9, 10, 12)8K3 = (7, 11, 13, 14, 15)91011def tile(c):12    return [((c >> i) & 1) | (((c >> (i + 2)) & 1) << 1) for i in range(2)]131415def kron(A, wa, B, wb):16    out = []17    for ra in A:18        for rb in B:19            v = 020            for j in range(wa):21                if (ra >> j) & 1:22                    v |= rb << (j * wb)23            out.append(v)24    return out252627def grid(w):28    rows, wid = [1], 129    for c in w:30        rows = kron(rows, wid, tile(c), 2)31        wid *= 232    return rows, wid333435def runs(row, n):36    out = []37    j = 038    while j < n:39        if (row >> j) & 1:40            k = j41            while k < n and (row >> k) & 1:42                k += 143            out.append((j, k))44            j = k45        else:46            j += 147    return out484950def label(rows, n, diag):51    parent = []52    segs = []5354    def find(x):55        while parent[x] != x:56            parent[x] = parent[parent[x]]57            x = parent[x]58        return x5960    prev = []61    for i in range(n):62        cur = []63        for (a, b) in runs(rows[i], n):64            idx = len(parent)65            parent.append(idx)66            segs.append((i, a, b))67            for (pa, pb, pi) in prev:68                touch = (pa <= b and a <= pb) if diag else (pa < b and a < pb)69                if touch:70                    ra, rb = find(pi), find(idx)71                    if ra != rb:72                        parent[rb] = ra73            cur.append((a, b, idx))74        prev = cur75    return parent, segs, find767778# OBSERVABLES7980def components(rows, n, diag=False):81    parent, segs, find = label(rows, n, diag)82    return len({find(i) for i in range(len(parent))})838485def enclosed(rows, n, diag):86    full = (1 << n) - 187    parent, segs, find = label([full ^ r for r in rows], n, diag)88    open_roots = set()89    roots = set()90    for i, (r, a, b) in enumerate(segs):91        root = find(i)92        roots.add(root)93        if r == 0 or r == n - 1 or a == 0 or b == n:94            open_roots.add(root)95    return len(roots - open_roots)969798def border_components(rows, n):99    parent, segs, find = label(rows, n, False)100    touching = set()101    for i, (r, a, b) in enumerate(segs):102        if r == 0 or r == n - 1 or a == 0 or b == n:103            touching.add(find(i))104    return len(touching)105106107def fill(rows):108    return sum(bin(r).count("1") for r in rows)109110111def perimeter(rows, n):112    h = sum(bin(r & (r >> 1)).count("1") for r in rows)113    v = sum(bin(rows[i] & rows[i + 1]).count("1") for i in range(n - 1))114    return 4 * fill(rows) - 2 * h - 2 * v115116117def interior(rows, n):118    band = ((1 << (n - 1)) - 2) if n >= 2 else 0119    total = 0120    for i in range(1, n - 1):121        m = rows[i] & (rows[i] >> 1) & (rows[i] << 1) & rows[i - 1] & rows[i + 1] & band122        total += bin(m).count("1")123    return total124125126def boundary(rows, n):127    return fill(rows) - interior(rows, n)128129130def maindiag(rows, n):131    return sum((rows[i] >> i) & 1 for i in range(n))132133134def antidiag(rows, n):135    return sum((rows[i] >> (n - 1 - i)) & 1 for i in range(n))136137138def antiprofile(rows, n):139    p = [0] * (2 * n - 1)140    for i in range(n):141        for j in range(n):142            if (rows[i] >> j) & 1:143                p[i + j] += 1144    return tuple(p)145146147def hcontact(rows, n):148    return sum(1 for i in range(n) if (rows[i] >> (n - 1)) & 1 and rows[i] & 1)149150151def vcontact(rows, n):152    return bin(rows[0] & rows[n - 1]).count("1")153154155NAMES = ("components", "boundary", "euler", "holes", "perimeter", "maindiag", "antidiag", "antiprofile", "fill")156157158def observe(w):159    rows, n = grid(w)160    c = components(rows, n)161    return (c, boundary(rows, n), c - enclosed(rows, n, True), enclosed(rows, n, False),162            perimeter(rows, n), maindiag(rows, n), antidiag(rows, n), antiprofile(rows, n), fill(rows))163164165def landscape(codes, L):166    vals = {w: observe(w) for w in product(codes, repeat=L)}167    groups = {}168    for w in vals:169        groups.setdefault(tuple(sorted(w)), []).append(w)170    multi = [v for v in groups.values() if len(v) >= 2]171    sens = {}172    for i, name in enumerate(NAMES):173        sens[name] = sum(1 for v in multi if len({vals[w][i] for w in v}) > 1)174    return vals, len(groups), len(multi), sens175176177def first_sensitive(vals, index):178    groups = {}179    for w in vals:180        groups.setdefault(tuple(sorted(w)), []).append(w)181    for key in sorted(groups):182        v = groups[key]183        if len(v) >= 2 and len({vals[w][index] for w in v}) > 1:184            return key185    return None186187188# COCYCLE189190CLASSES = {191    (1, 2, 4, 8): ((1, 0, 0, 0), (2, 0, 0, 0), (2, 0, 0, 0), (4, 0, 0, 0)),192    (6, 9): ((2, 0, 0, 0), (4, 0, 0, 0), (4, 0, 0, 0), (8, 0, 0, 0)),193    (3, 12): ((0, 1, 0, 0), (-2, 3, 0, 0), (0, 2, 0, 0), (-4, 6, 0, 0)),194    (5, 10): ((0, 0, 1, 0), (0, 0, 2, 0), (-2, 0, 3, 0), (-4, 0, 6, 0)),195    (7, 11, 13, 14): ((-1, 1, 1, 0), (-4, 3, 2, 0), (-4, 2, 3, 0), (-8, 4, 4, 1)),196    (15,): ((0, 0, 0, 1), (0, 0, -2, 3), (0, -2, 0, 3), (4, -6, -6, 9)),197}198MAT = {}199for klass, m in CLASSES.items():200    for c in klass:201        MAT[c] = m202203204def matmul(A, B):205    return tuple(tuple(sum(A[i][t] * B[t][j] for t in range(4)) for j in range(4)) for i in range(4))206207208def represent(w):209    v = (1, 0, 0, 0)210    for c in w:211        M = MAT[c]212        v = tuple(sum(v[i] * M[i][j] for i in range(4)) for j in range(4))213    return sum(v)214215216PRIME = (1 << 31) - 1217218219def hankel_rank(table, prefixes, index):220    rows = [[table[u + v][index] % PRIME for v in prefixes] for u in prefixes]221    r, col, m, ncol = 0, 0, len(rows), len(prefixes)222    while col < ncol and r < m:223        piv = next((i for i in range(r, m) if rows[i][col]), None)224        if piv is None:225            col += 1226            continue227        rows[r], rows[piv] = rows[piv], rows[r]228        inv = pow(rows[r][col], PRIME - 2, PRIME)229        pr = [x * inv % PRIME for x in rows[r]]230        rows[r] = pr231        for i in range(r + 1, m):232            f = rows[i][col]233            if f:234                rows[i] = [(a - f * b) % PRIME for a, b in zip(rows[i], pr)]235        r += 1236        col += 1237    return r238239240# CHECKS241242def check(name, got, want):243    assert got == want, f"{name}: got {got}, want {want}"244    print(f"{name}: {got}")245246247def main():248    check("comp(3,6)", observe((3, 6))[0], 4)249    check("comp(6,3)", observe((6, 3))[0], 2)250251    sensitive_k2 = sorted(p for p in product(K2, repeat=2) if p[0] < p[1]252                          and observe(p)[0] != observe((p[1], p[0]))[0])253    check("least order-sensitive pair of two-cell designs", sensitive_k2[0], (3, 6))254255    bad = 0256    for a in range(16):257        for b in range(16):258            ra, na = grid((a, b))259            rb, nb = grid((b, a))260            if boundary(ra, na) != boundary(rb, nb):261                bad += 1262    check("boundary-cell asymmetries over all 256 ordered code pairs", bad, 0)263264    def three(c):265        return frozenset({(0, 0), (0, 1), (1, 0), (1, 1)} - {c})266267    def cells(c):268        return {(i, j) for i in range(2) for j in range(2) if (c >> (i + 2 * j)) & 1}269270    bad = 0271    for a in range(16):272        for b in range(16):273            fa, fb = cells(a), cells(b)274            formula = sum(1 for p in [(0, 0), (0, 1), (1, 0), (1, 1)]275                          if three(p) <= fa and three((1 - p[0], 1 - p[1])) <= fb)276            rows, n = grid((a, b))277            if formula != interior(rows, n):278                bad += 1279    check("interior-formula mismatches over all 256 ordered code pairs", bad, 0)280281    bad = 0282    for w in product(CODES, repeat=3):283        rows, n = grid(w)284        prod = 1285        for c in w:286            prod *= sum(1 for d in range(2) if (c >> (d + 2 * d)) & 1)287        if maindiag(rows, n) != prod:288            bad += 1289    check("main-diagonal product-formula mismatches over all 3375 words of length 3", bad, 0)290291    bad = 0292    for w in product(CODES, repeat=3):293        rows, n = grid(w)294        ph, pv = 1, 1295        for c in w:296            t = tile(c)297            ph *= hcontact(t, 2)298            pv *= vcontact(t, 2)299        if (hcontact(rows, n), vcontact(rows, n)) != (ph, pv):300            bad += 1301    check("contact-law mismatches over all 3375 words of length 3", bad, 0)302303    adjacent, diagonal = (3, 5, 10, 12), (6, 9)304    bad = 0305    for a in adjacent:306        for d in diagonal:307            if (observe((a, d))[0], observe((d, a))[0]) != (4, 2):308                bad += 1309    check("adjacent-times-diagonal deviations from (4,2) over all 8 ordered pairs", bad, 0)310311    bad = sum(1 for a in adjacent for b in adjacent if observe((a, b))[0] != observe((b, a))[0])312    check("adjacent-times-adjacent asymmetries over all 16 ordered pairs", bad, 0)313    bad = sum(1 for a in diagonal for b in diagonal if observe((a, b))[0] != 4)314    check("diagonal-times-diagonal deviations from 4 over all 4 ordered pairs", bad, 0)315    bad = sum(1 for a in K3 for b in K3 if observe((a, b))[0] != 1)316    check("k>=3 deviations from 1 over all 25 ordered pairs", bad, 0)317318    bad = 0319    for c1 in CODES:320        for c2 in CODES:321            for L in (2, 3):322                block = grid((c1, c2))323                rows, n = block324                for _ in range(L - 1):325                    rows = kron(rows, n, block[0], block[1])326                    n *= block[1]327                if (rows, n) != grid((c1, c2) * L):328                    bad += 1329    check("block-reduction mismatches over 225 blocks at periods 2 and 3", bad, 0)330331    v2, g2, m2, s2 = landscape(CODES, 2)332    check("length-2 multisets", (len(v2), g2, m2), (225, 120, 105))333    check("length-2 order-sensitive multisets", tuple(s2[k] for k in NAMES),334          (74, 0, 78, 10, 78, 0, 0, 99, 0))335336    v3, g3, m3, s3 = landscape(LIB3, 3)337    check("length-3 multisets over the ten-code library", (len(v3), g3, m3), (1000, 220, 210))338    check("length-3 order-sensitive multisets", tuple(s3[k] for k in NAMES),339          (188, 36, 188, 100, 188, 0, 0, 204, 0))340341    check("length-2 words with euler != components - holes",342          sum(1 for o in v2.values() if o[2] != o[0] - o[3]), 10)343    check("length-3 words with euler != components - holes",344          sum(1 for o in v3.values() if o[2] != o[0] - o[3]), 168)345346    check("anti-diagonal profiles of (1,2) and (2,1)", (observe((1, 2))[7], observe((2, 1))[7]),347          ((0, 1, 0, 0, 0, 0, 0), (0, 0, 1, 0, 0, 0, 0)))348    check("lexicographically first anti-diagonal-profile witness", first_sensitive(v2, 7), (1, 2))349    check("perimeters of (1,3) and (3,1)", (observe((1, 3))[4], observe((3, 1))[4]), (6, 8))350    check("lexicographically first perimeter witness", first_sensitive(v2, 4), (1, 3))351    check("boundary cells of (1,3) and (3,1)", (observe((1, 3))[1], observe((3, 1))[1]), (2, 2))352353    table = {(): (1, 1, 1, 0)}354    peak = {}355    for L in (1, 2, 3, 4):356        best = 0357        for w in product(CODES, repeat=L):358            rows, n = grid(w)359            c = components(rows, n)360            table[w] = (c, boundary(rows, n), c - enclosed(rows, n, True), enclosed(rows, n, False))361            best = max(best, c)362        peak[L] = best363    check("words of length at most 4", len(table), 54241)364    check("representation mismatches over all 54241 words of length at most 4",365          sum(1 for w, o in table.items() if represent(w) != o[0]), 0)366    check("largest component count at lengths 1 to 4", tuple(peak[L] for L in (1, 2, 3, 4)),367          (2, 8, 32, 128))368369    for L in range(1, 7):370        rows, n = grid((15,) * (L - 1) + (6,))371        check(f"components of (15^{L - 1}, 6)", components(rows, n), 2 * 4 ** (L - 1))372    for L in range(1, 8):373        rows, n = grid((15,) * (L - 1) + (3,))374        check(f"border-meeting components of (15^{L - 1}, 3)", border_components(rows, n), 2 ** (L - 1))375376    prefixes = [()] + [(a,) for a in CODES] + [(a, b) for a in CODES for b in CODES]377    check("Hankel ranks of the length-4 truncation",378          tuple(hankel_rank(table, prefixes, i) for i in (0, 2, 1, 3)), (4, 4, 8, 11))379380    klass = list(CLASSES)381    pairs = [(i, j) for i in range(6) for j in range(i + 1, 6)]382    noncommuting = [(klass[i], klass[j]) for (i, j) in pairs383                    if matmul(CLASSES[klass[i]], CLASSES[klass[j]]) != matmul(CLASSES[klass[j]], CLASSES[klass[i]])]384    check("noncommuting class pairs out of 15", len(noncommuting), 14)385    commuting = [(klass[i], klass[j]) for (i, j) in pairs if (klass[i], klass[j]) not in noncommuting]386    check("commuting class pair", commuting, [((1, 2, 4, 8), (6, 9))])387    check("M(6,9) equals twice M(1,2,4,8)", CLASSES[(6, 9)],388          tuple(tuple(2 * x for x in row) for row in CLASSES[(1, 2, 4, 8)]))389390    for (a, b, want) in ((3, 5, 2), (7, 15, 1)):391        check(f"M({a}) and M({b}) commute", matmul(MAT[a], MAT[b]) == matmul(MAT[b], MAT[a]), False)392        check(f"components of ({a},{b}) and of ({b},{a})",393              (observe((a, b))[0], observe((b, a))[0]), (want, want))394395    print("all green")396397398if __name__ == "__main__":399    main()