controls.py

1import itertools2import math3from fractions import Fraction45CODES = range(1, 16)67def tile(code):8    return {(j & 1, (j >> 1) & 1) for j in range(4) if code & (1 << j)}910def kron(a, b, side_b):11    return {(ra * side_b + rb, ca * side_b + cb)12            for ra, ca in a for rb, cb in b}1314def word_array(word):15    cells, side = tile(word[0]), 216    for code in word[1:]:17        cells = kron(cells, tile(code), 2)18        side *= 219    return cells, side2021def profile(cells, side):22    p = [0] * (2 * side - 1)23    for r, c in cells:24        p[r + c] += 125    return p2627def convolve(p, q):28    out = [0] * (len(p) + len(q) - 1)29    for i, x in enumerate(p):30        if x:31            for j, y in enumerate(q):32                out[i + j] += x * y33    return out3435def dilate(p, k):36    out = [0] * ((len(p) - 1) * k + 1)37    for i, x in enumerate(p):38        out[i * k] = x39    return out4041def report_profile_identity():42    bad = 043    for length in (2, 3):44        for word in itertools.product(CODES, repeat=length):45            head, side = word_array(word[:-1])46            tail = tile(word[-1])47            lhs = profile(kron(head, tail, 2), side * 2)48            rhs = convolve(dilate(profile(head, side), 2), profile(tail, 2))49            bad += lhs != rhs50    print(f"profile identity mismatches over words of length 2 and 3: {bad}")5152def digit_polynomial(D):53    p = [0] * (2 * D + 1)54    for v in itertools.product(range(3), repeat=D):55        if sum(1 for x in v if x == 1) <= 1:56            p[sum(v)] += 157    return p5859def central_slice(D, level):60    p = digit_polynomial(D)61    acc = [1]62    for j in range(level):63        acc = convolve(acc, dilate(p, 3 ** j))64    return acc[(len(acc) - 1) // 2]6566def cross_section_vertices(D):67    if D % 2 == 0:68        return math.comb(D, D // 2)69    return math.comb(D, (D - 1) // 2) * ((D + 1) // 2)7071def report_vertex_identity(top):72    row, bad = [], 073    for D in range(2, top + 1):74        lhs, rhs = central_slice(D, 1), cross_section_vertices(D)75        bad += lhs != rhs76        row.append(lhs)77    print(f"level-1 slice against cross-section vertices, D = 2..{top}: "78          f"{bad} mismatches")79    print(f"level-1 slice counts D = 2..{top}: {row}")8081def fit_order_two(seq):82    s0, s1, s2, s3 = (Fraction(x) for x in seq[:4])83    det = s1 * s1 - s0 * s284    return (s1 * s2 - s0 * s3) / det, (s1 * s3 - s2 * s2) / det8586def report_rung(D, levels):87    seq = [central_slice(D, level) for level in range(1, levels + 1)]88    c1, c2 = fit_order_two(seq)89    holds = all(seq[n] == c1 * seq[n - 1] + c2 * seq[n - 2]90                for n in range(2, len(seq)))91    root = (float(c1) + math.sqrt(float(c1 * c1 + 4 * c2))) / 292    print(f"D = {D} census, levels 1..{levels}: {seq}")93    print(f"D = {D} recurrence a(n) = {c1}a(n-1) + {c2}a(n-2), "94          f"fitted on four terms, holds on every term: {holds}")95    print(f"D = {D} dominant root {root:.9f}, "96          f"slice dimension {math.log(root) / math.log(3):.9f}")9798def report_staircase(bases, top):99    for n in range(1, top + 1):100        num = sum((n - j) * math.log(q * q - ((q - 1) // 2) ** 2)101                  for j, q in enumerate(bases[:n]))102        den = sum((n - j) * math.log(q) for j, q in enumerate(bases[:n]))103        print(f"staircase n = {n}, bases {bases[:n]}, "104              f"dimension {num / den:.9f}")105106def main():107    print("fill assumed for base q is q^2 - ((q-1)/2)^2")108    report_profile_identity()109    report_vertex_identity(14)110    report_rung(4, 6)111    report_staircase([3, 5, 7, 9, 11], 5)112113main()