verify.py

1from fractions import Fraction2from itertools import combinations3from math import gcd45# GASKET67GASKET = ((0, 0), (1, 0), (0, 1))8GASKET_SET = set(GASKET)9CELLS = [(a, b) for a in range(3) for b in range(3)]10PERMS = [(0, 1, 2), (0, 2, 1), (1, 0, 2), (1, 2, 0), (2, 0, 1), (2, 1, 0)]111213def die(name, got, want):14    raise AssertionError("%s: got %r, want %r" % (name, got, want))151617def check(name, got, want):18    if got != want:19        die(name, got, want)202122def points(design, level):23    pts = [(0, 0)]24    for i in range(level):25        p = 3 ** i26        pts = [(x + dx * p, y + dy * p) for (x, y) in pts for (dx, dy) in design]27    return pts282930def rays(pts):31    out = set()32    for (x, y) in pts:33        if x == 0 and y == 0:34            continue35        g = gcd(x, y)36        out.add((x // g, y // g))37    return out383940def diagonal_to(design, level):41    seen = set()42    for L in range(1, level + 1):43        pts = points(design, L)44        n = sum(1 for p in pts if p != (0, 0))45        seen = rays(pts)46        if len(seen) != n:47            return False, 048    return True, sum(1 for (x, y) in seen if x == 0 or y == 0)495051def code(design):52    return sum(2 ** (3 * a + b) for (a, b) in design)535455def graph(phi):56    return tuple((j, phi[j]) for j in range(3))575859def cross_coefficient(phi, d0, u, v):60    return (phi[d0] * (u - v) - d0 * (phi[u] - phi[v])) % 3616263def mass_edges(a, b):64    index = {0: 0}65    order = [0]66    edges = []67    i = 068    while i < len(order):69        c = order[i]70        out = []71        for (dx, dy) in GASKET:72            v = c + b * dx - a * dy73            if v % 3 == 0:74                w = v // 375                if w not in index:76                    index[w] = len(order)77                    order.append(w)78                out.append(index[w])79        edges.append(out)80        i += 181    return edges828384def mass(a, b, n):85    cur = {0: 1}86    for _ in range(n):87        nxt = {}88        for c, k in cur.items():89            for (dx, dy) in GASKET:90                v = c + b * dx - a * dy91                if v % 3 == 0:92                    w = v // 393                    nxt[w] = nxt.get(w, 0) + k94        cur = nxt95    return cur.get(0, 0) - 1969798def mass_brute(a, b, n):99    total = 0100    for (x, y) in points(GASKET, n):101        if (x, y) != (0, 0) and x * b == y * a:102            total += 1103    return total104105106def mass_split(a, b, n):107    half = n // 2108    table = {}109    for (x, y) in points(GASKET, half):110        v = b * x - a * y111        table[v] = table.get(v, 0) + 1112    p = 3 ** half113    total = 0114    for (x, y) in points(GASKET, n - half):115        total += table.get(-p * (b * x - a * y), 0)116    return total - 1117118119FIB = [0, 1, 1]120LUC = [2, 1]121122123def fibonacci(k):124    while len(FIB) <= k:125        FIB.append(FIB[-1] + FIB[-2])126    return FIB[k]127128129def lucas(k):130    while len(LUC) <= k:131        LUC.append(LUC[-1] + LUC[-2])132    return LUC[k]133134135def shift_product(j, n):136    if n <= j:137        return 1138    total = 1139    for r in range(j):140        m = len(range(r, n - j, j))141        total *= fibonacci(m + 2)142    return total143144145def narayana(n):146    a = [1, 1, 1]147    while len(a) <= n:148        a.append(a[-1] + a[-3])149    return a[n]150151152def quartic(n):153    c = [1, 2, 3, 4]154    while len(c) <= n:155        c.append(c[-1] + c[-4])156    return c[n]157158159def automaton(s, t):160    start = (0, 0, 0, 0)161    index = {start: 0}162    order = [start]163    edges = []164    i = 0165    while i < len(order):166        a1, a2, b1, b2 = order[i]167        out = []168        for (dx, dy) in GASKET:169            u = ((s * dx + a1) % 3, (s * dy + a2) % 3)170            v = ((t * dx + b1) % 3, (t * dy + b2) % 3)171            if u in GASKET_SET and v in GASKET_SET:172                nxt = ((s * dx + a1) // 3, (s * dy + a2) // 3,173                       (t * dx + b1) // 3, (t * dy + b2) // 3)174                if nxt not in index:175                    index[nxt] = len(order)176                    order.append(nxt)177                out.append(index[nxt])178        edges.append(out)179        i += 1180    return edges181182183def free_automaton(s, t):184    start = (0, 0, 0, 0)185    index = {start: 0}186    order = [start]187    edges = []188    i = 0189    while i < len(order):190        a1, a2, b1, b2 = order[i]191        out = []192        for dx in range(3):193            for dy in range(3):194                u = ((s * dx + a1) % 3, (s * dy + a2) % 3)195                v = ((t * dx + b1) % 3, (t * dy + b2) % 3)196                if u in GASKET_SET and v in GASKET_SET:197                    nxt = ((s * dx + a1) // 3, (s * dy + a2) // 3,198                           (t * dx + b1) // 3, (t * dy + b2) // 3)199                    if nxt not in index:200                        index[nxt] = len(order)201                        order.append(nxt)202                    out.append(index[nxt])203        edges.append(out)204        i += 1205    return edges206207208def live(edges):209    back = [[] for _ in edges]210    for i, outs in enumerate(edges):211        for e in outs:212            back[e].append(i)213    seen = {0}214    stack = [0]215    while stack:216        v = stack.pop()217        for u in back[v]:218            if u not in seen:219                seen.add(u)220                stack.append(u)221    keep = sorted(seen)222    place = {v: i for i, v in enumerate(keep)}223    return [[place[e] for e in edges[v] if e in place] for v in keep]224225226def pair_brute(s, t, n):227    pts = set(points(GASKET, n))228    return sum(1 for (x, y) in pts229               if (s * x, s * y) in pts and (t * x, t * y) in pts)230231232def return_counts(edges, m):233    count = [0] * len(edges)234    count[0] = 1235    out = []236    for _ in range(m + 1):237        out.append(count[0])238        nxt = [0] * len(edges)239        for i, outs in enumerate(edges):240            for e in outs:241                nxt[e] += count[i]242        count = nxt243    return out244245246def charpoly(edges):247    n = len(edges)248    m = [[1 if i == j else 0 for j in range(n)] for i in range(n)]249    coeffs = [1]250    for k in range(1, n + 1):251        am = [[0] * n for _ in range(n)]252        for i, outs in enumerate(edges):253            row = am[i]254            for e in outs:255                src = m[e]256                for j in range(n):257                    row[j] += src[j]258        trace = sum(am[i][i] for i in range(n))259        if trace % k != 0:260            die("Faddeev-LeVerrier integrality", trace % k, 0)261        ck = -(trace // k)262        coeffs.append(ck)263        for i in range(n):264            am[i][i] += ck265        m = am266    return coeffs267268269def shift_coefficients(coeffs, r):270    p = list(coeffs)271    out = []272    for _ in range(len(coeffs)):273        acc = 0274        nb = []275        for a in p:276            acc = acc * r + a277            nb.append(acc)278        out.append(nb[-1])279        p = nb[:-1]280    return out281282283def evaluate(coeffs, x):284    acc = 0285    for a in coeffs:286        acc = acc * x + a287    return acc288289290def divides(d, f):291    f = list(f)292    while len(f) >= len(d):293        if f[0] % d[0] != 0:294            return False295        q = f[0] // d[0]296        for i in range(len(d)):297            f[i] -= q * d[i]298        if f[0] != 0:299            return False300        f.pop(0)301    return all(c == 0 for c in f)302303304def no_root_above(coeffs, r):305    return all(v >= 0 for v in shift_coefficients(coeffs, r))306307308def bisect(coeffs, steps):309    lo, hi = Fraction(0), Fraction(4)310    for _ in range(steps):311        mid = (lo + hi) / 2312        if no_root_above(coeffs, mid):313            hi = mid314        else:315            lo = mid316    return lo, hi317318319def in_window(coeffs, lo, hi):320    return evaluate(coeffs, lo) < 0 and no_root_above(coeffs, hi)321322323def word_counts(edges, m):324    count = [0] * len(edges)325    count[0] = 1326    out = []327    for _ in range(m + 1):328        out.append(sum(count))329        nxt = [0] * len(edges)330        for i, outs in enumerate(edges):331            for e in outs:332                nxt[e] += count[i]333        count = nxt334    return out335336337# CENSUS OF DESIGNS338339def census_designs():340    survivors = []341    for design in combinations(CELLS, 3):342        ok, fibres = diagonal_to(design, 6)343        if ok:344            survivors.append((design, fibres))345    check("no-two-collinear subsets to level 6", len(survivors), 10)346    perms = [(d, f) for (d, f) in survivors347             if len(set(a for a, b in d)) == 3 and len(set(b for a, b in d)) == 3]348    check("permutation designs among survivors", len(perms), 4)349    check("survivors with exactly two fibre rays",350          sorted(d for (d, f) in survivors if f == 2), sorted(d for (d, f) in perms))351    want_codes = sorted([98, 140, 266, 84])352    check("codes of the four diagonal designs",353          sorted(code(d) for (d, f) in perms), want_codes)354    check("code 148 is not a permutation design",355          sorted(c for c in CELLS if 148 >> (3 * c[0] + c[1]) & 1),356          [(0, 2), (1, 1), (2, 1)])357    good = sorted(graph(p) for p in PERMS if p[0] != 0)358    check("diagonal designs are exactly the phi(0) nonzero graphs",359          sorted(tuple(sorted(d)) for (d, f) in perms),360          sorted(tuple(sorted(g)) for g in good))361    extra = sorted(d for (d, f) in survivors if f != 2)362    check("the six non-permutation survivors", extra, [363        ((0, 1), (1, 1), (2, 1)),364        ((0, 2), (1, 2), (2, 2)),365        ((1, 0), (1, 1), (1, 2)),366        ((1, 1), (1, 2), (2, 1)),367        ((1, 2), (2, 1), (2, 2)),368        ((2, 0), (2, 1), (2, 2)),369    ])370    print("84 subsets to level 6: 10 with no two points collinear, 4 of them permutation designs")371372373def census_cross():374    for phi in PERMS:375        got = set()376        for d0 in range(3):377            for u in range(3):378                for v in range(3):379                    if u != v:380                        got.add(cross_coefficient(phi, d0, u, v) != 0)381                        check("affine collapse of the cross coefficient",382                              cross_coefficient(phi, d0, u, v), (phi[0] * (u - v)) % 3)383        check("18 cross cases for phi=%r" % (phi,), got, {phi[0] != 0})384    print("cross coefficient: nonzero in all 18 cases exactly when phi(0) is nonzero")385386387def census_diagonality():388    for phi in PERMS:389        if phi[0] == 0:390            continue391        design = graph(phi)392        for n in range(1, 9):393            pts = points(design, n)394            check("Z(n) for phi=%r at n=%d" % (phi, n),395                  len(rays(pts)) - 2, 3 ** n - 2)396        check("Z at n=0 for phi=%r" % (phi,), len(rays(points(design, 0))), 0)397    identity = points(graph((0, 1, 2)), 1)398    check("identity witness at level 1", sorted(identity), [(0, 0), (1, 1), (2, 2)])399    check("identity determinant", 1 * 2 - 1 * 2, 0)400    doubling = points(graph((0, 2, 1)), 2)401    check("doubling witnesses at level 2", ((1, 2) in doubling, (3, 6) in doubling), (True, True))402    check("doubling determinant", 1 * 6 - 2 * 3, 0)403    print("Z_F(n) = 3^n - 2 for all four diagonal designs, n = 1..8; both phi(0) = 0 designs collapse")404405406# RAY MASSES407408def census_masses():409    got = [mass(3, 1, n) for n in range(1, 31)]410    check("M(3,1) to n=30", got, [fibonacci(n + 1) - 1 for n in range(1, 31)])411    got = [mass(1, 12, n) for n in range(1, 31)]412    check("M(1,12) to n=30", got, [narayana(n) - 1 for n in range(1, 31)])413    got = [mass(7, 3, n) for n in range(3, 31)]414    check("M(7,3) to n=30", got, [quartic(n - 3) - 1 for n in range(3, 31)])415    for j in range(1, 14):416        got = [mass(3 ** j, 1, n) for n in range(1, 31)]417        check("M(3^%d,1) to n=30" % j, got,418              [shift_product(j, n) - 1 for n in range(1, 31)])419    for (a, b) in [(3, 1), (1, 12), (7, 3), (9, 1), (27, 1), (12, 13)]:420        for n in range(1, 11):421            check("automaton against brute force at (%d,%d), n=%d" % (a, b, n),422                  mass(a, b, n), mass_brute(a, b, n))423        for n in range(1, 15):424            check("automaton against split enumeration at (%d,%d), n=%d" % (a, b, n),425                  mass(a, b, n), mass_split(a, b, n))426    squares = [(n, mass(7, 3, n)) for n in range(1, 31)427               if mass(7, 3, n) > 0 and int(round(mass(7, 3, n) ** 0.5)) ** 2 == mass(7, 3, n)]428    check("perfect-square masses on (7,3) to n=30", squares,429          [(4, 1), (7, 4), (9, 9), (12, 25), (14, 49)])430    check("rational roots of x^4 - x^3 - 1",431          [r for r in (1, -1) if r ** 4 - r ** 3 - 1 == 0], [])432    check("integer quadratic factorisations of x^4 - x^3 - 1",433          [(a, b, c, d) for b, d in ((1, -1), (-1, 1))434           for a in range(-9, 10) for c in range(-9, 10)435           if a + c == -1 and b + d + a * c == 0 and a * d + b * c == 0], [])436    print("mass laws on (3,1), (1,12), (7,3) and the thirteen shift rays: exact to n = 30")437438439def nonneg(x, y):440    if x >= 0 and y >= 0:441        return True442    if x < 0 and y < 0:443        return False444    if y >= 0:445        return 5 * y * y >= x * x446    return x * x >= 5 * y * y447448449def power_root5(a, b, j):450    p, q = 1, 0451    for _ in range(j):452        p, q = a * p + 5 * b * q, a * q + b * p453    return p, q454455456def shift_strings(n, j):457    total = 0458    for z in range(2 ** (n - j)):459        d = [(z >> i) & 1 for i in range(n - j)]460        if all(not (d[i] and d[i - j]) for i in range(j, len(d))):461            total += 1462    return total - 1463464465def census_recurrences():466    for (a, b, want, seeds) in [467            (3, 1, [1, -1, -1], [fibonacci(n + 1) for n in range(0, 3)]),468            (1, 12, [1, -1, 0, -1], [narayana(n) for n in range(0, 4)]),469            (7, 3, [1, -1, 0, 0, -1], [quartic(n) for n in range(0, 5)])]:470        edges = live(mass_edges(a, b))471        check("live carry states at (%d,%d)" % (a, b), len(edges), len(want) - 1)472        check("carry characteristic polynomial at (%d,%d)" % (a, b),473              charpoly(edges), want)474        k = len(want) - 1475        got = return_counts(edges, 2 * k + 12)476        off = 3 if (a, b) == (7, 3) else 0477        check("automaton seeds at (%d,%d)" % (a, b), got[off:off + k], seeds[:k])478        for m in range(k, len(got)):479            check("annihilator residual at (%d,%d), m=%d" % (a, b, m),480                  sum(want[i] * got[m - i] for i in range(k + 1)), 0)481    print("mass recurrences are Cayley-Hamilton on 2, 3 and 4 live carry states")482483484def census_shift():485    for k in range(2, 400):486        check("F(%d) <= 2 F(%d)" % (k, k - 1),487              fibonacci(k) <= 2 * fibonacci(k - 1), True)488    for n in range(1, 12):489        for j in range(1, n + 1):490            check("shift strings at n=%d, j=%d" % (n, j),491                  (mass(3 ** j, 1, n), shift_strings(n, j)),492                  (shift_product(j, n) - 1, shift_product(j, n) - 1))493    for n in range(1, 61):494        for j in range(1, n + 1):495            q, rem = divmod(n - j, j)496            check("block form at n=%d, j=%d" % (n, j), shift_product(j, n),497                  fibonacci(q + 3) ** rem * fibonacci(q + 2) ** (j - rem))498            p, r = power_root5(3, -1, j)499            m = shift_product(j, n) - 1500            check("per-j bound at n=%d, j=%d" % (n, j),501                  nonneg(p * lucas(n) + 5 * r * fibonacci(n) - 2 * m,502                         p * fibonacci(n) + r * lucas(n)), True)503    for n in range(1, 161):504        sh = 2 * sum((shift_product(j, n) - 1) ** 2 for j in range(1, n + 1))505        el, ef = lucas(2 * n), fibonacci(2 * n)506        check("family bound at n=%d" % n,507              nonneg(4 * el + 60 * ef - 22 * sh, 12 * el + 4 * ef), True)508        if n >= 100:509            check("family limit window below at n=%d" % n,510                  nonneg(2000000 * sh - 2198212 * el, -2198212 * ef), True)511            check("family limit window above at n=%d" % n,512                  nonneg(2198214 * el - 2000000 * sh, 2198214 * ef), True)513    print("shift family: M_n(3^j,1) < (3 - sqrt5)^j phi^n and sum of squares "514          "in (2.198212, 2.198214) phi^2n from n = 100")515516517def census_gasket(n):518    pts = points(GASKET, n)519    check("gasket size at n=%d" % n, len(pts), 3 ** n)520    counts = {}521    for (x, y) in pts:522        if x == 0 or y == 0:523            continue524        g = gcd(x, y)525        r = (x // g, y // g)526        counts[r] = counts.get(r, 0) + 1527    total = sum(counts.values())528    check("non-fibre mass at n=%d" % n, total, 3 ** n - 2 ** (n + 1) + 1)529    return counts530531532def census_rays():533    second, residual = [], []534    for n in range(1, 13):535        counts = census_gasket(n)536        z = sum(m * m for m in counts.values())537        second.append(z)538        residual.append(z - (3 ** n - 2 ** (n + 1) + 1)539                        - (3 ** n - 4 * 2 ** n + 2 * n + 3))540        check("shift-ray second moment at n=%d" % n,541              sum(m * m for r, m in counts.items() if is_shift_pair(*r)),542              2 * sum((shift_product(j, n) - 1) ** 2 for j in range(1, n + 1)))543    for n in range(1, 21):544        check("shift-multiplier closed form at n=%d" % n,545              2 * sum(3 ** (n - j) - 2 ** (n - j + 1) + 1 for j in range(1, n)),546              3 ** n - 4 * 2 ** n + 2 * n + 3)547    check("E(n) for n = 1..12", second,548          [0, 2, 16, 98, 396, 1522, 5248, 17118, 52212, 158042, 466960, 1374038])549    check("R(n) for n = 1..12", residual,550          [0, 0, 0, 20, 88, 432, 1624, 5512, 15896, 46064, 124928, 335704])551    check("occupied non-fibre rays at n=12", len(counts), 345318)552    check("non-fibre mass at n=12", sum(counts.values()), 523250)553    top = sorted(counts.items(), key=lambda kv: (-kv[1], kv[0]))[:10]554    check("ten heaviest rays at n=12", top, [555        ((1, 3), 232), ((3, 1), 232), ((1, 9), 168), ((9, 1), 168),556        ((1, 27), 124), ((27, 1), 124), ((1, 81), 80), ((81, 1), 80),557        ((1, 243), 71), ((243, 1), 71)])558    check("mirror symmetry of ray masses at n=12",559          all(counts[(b, a)] == m for (a, b), m in counts.items()), True)560    print("gasket census to n = 12: 345318 occupied non-fibre rays carrying 523250 points, "561          "E(12) = 1374038, R(12) = 335704")562563564# PAIR CENSUS565566def ray_masses(n):567    out = {}568    for (x, y) in points(GASKET, n):569        if x == 0 or y == 0:570            continue571        g = gcd(x, y)572        r = (x // g, y // g)573        if r in out:574            out[r].append(g)575        else:576            out[r] = [g]577    return out578579580def is_shift_pair(s, t):581    if min(s, t) != 1:582        return False583    u = max(s, t)584    while u % 3 == 0:585        u //= 3586    return u == 1587588589def census_pairs():590    n = 9591    pts = set(points(GASKET, n))592    by_ray = ray_masses(n)593    check("occupied non-fibre rays at n=9", len(by_ray), 12170)594    Z = sum(len(v) ** 2 for v in by_ray.values())595    T = 3 ** n - 2 ** (n + 1) + 1596    S = 3 ** n - 4 * 2 ** n + 2 * n + 3597    check("E(9), diagonal, 3-power family and residual at n=9",598          (Z, T, S, Z - T - S), (52212, 18660, 17656, 15896))599    check("shift-ray share of the second moment at n=9",600          sum(len(v) ** 2 for r, v in by_ray.items() if is_shift_pair(*r)), 11852)601    census, inside = {}, {}602    for ray, gs in by_ray.items():603        for g in gs:604            for h in gs:605                if g == h:606                    continue607                d = gcd(g, h)608                key = (g // d, h // d)609                census[key] = census.get(key, 0) + 1610                hit = 1 if (d * ray[0], d * ray[1]) in pts else 0611                inside[key] = inside.get(key, 0) + hit612    check("ordered off-diagonal collinear pairs at n=9",613          sum(census.values()), Z - T)614    check("shift-multiplier family at n=9",615          sum(v for k, v in census.items() if is_shift_pair(*k)), S)616    check("active ordered multiplier pairs at n=9", len(census), 2656)617    check("active unordered multiplier pairs at n=9", len(census) // 2, 1328)618    check("shift pairs among them",619          sorted(k for k in census if is_shift_pair(*k)),620          [(1, 3 ** j) for j in range(1, 8)] + [(3 ** j, 1) for j in range(1, 8)])621    check("largest multiplier at n=9", max(max(k) for k in census), 2460)622    short = [k for k in census if census[k] != inside[k]]623    check("ordered pairs whose witness leaves the gasket", len(short), 482)624    check("witnesses outside the gasket",625          sum(census[k] - inside[k] for k in short), 2540)626    check("every such pair has both multipliers above one",627          min(min(k) for k in short), 2)628    worst = max(census[k] - inside[k] for k in short)629    check("worst such pairs",630          (worst, sorted(k for k in short if census[k] - inside[k] == worst),631           inside[(41, 122)]),632          (50, [(41, 122), (122, 41), (122, 123), (123, 122)], 0))633    binary = set()634    for u in range(2 ** n):635        binary.add(sum(((u >> i) & 1) * 3 ** i for i in range(n)))636    for (s, t) in sorted(k for k in census if k[0] < k[1]):637        got = return_counts(live(automaton(s, t)), n)[n]638        fib = sum(1 for u in binary639                  if u > 0 and s * u in binary and t * u in binary)640        check("A(%d,%d) against brute force at n=9" % (s, t),641              got - 1 - 2 * fib, inside[(s, t)])642        if t > 52:643            continue644        got = return_counts(live(free_automaton(s, t)), n)[n]645        fib = sum(1 for u in range(1, 3 ** n // t + 1)646                  if s * u in binary and t * u in binary)647        check("B(%d,%d) against brute force at n=9" % (s, t),648              got - 1 - 2 * fib, census[(s, t)])649    small = [k for k in census if k[0] < k[1] and k[1] <= 52]650    check("unordered pairs of height at most 52 at n=9", len(small), 103)651    check("collinear pairs they carry", sum(census[k] for k in small), 11330)652    print("n = 9 pair census: 2656 active ordered multiplier pairs, 482 of them "653          "with a witness off the gasket")654655656# SPECTRAL GAP657658def census_spectrum():659    pairs = [(s, t) for s in range(1, 53) for t in range(s + 1, 53) if gcd(s, t) == 1]660    check("coprime unordered pairs with max at most 52", len(pairs), 829)661    reachable = {p: automaton(*p) for p in pairs}662    trimmed = {p: live(reachable[p]) for p in pairs}663    check("largest reachable state set", max(len(reachable[p]) for p in pairs), 45)664    check("largest live state set", max(len(trimmed[p]) for p in pairs), 33)665    check("(1,16) reachable states", len(reachable[(1, 16)]), 11)666    check("(1,16) live states", len(trimmed[(1, 16)]), 1)667    plastic = bisect(charpoly(reachable[(1, 16)]), 52)668    check("(1,16) reachable radius is the plastic number",669          (Fraction(13247179572, 10 ** 10) < plastic[0],670           plastic[1] < Fraction(13247179573, 10 ** 10)), (True, True))671    polys = {p: charpoly(trimmed[p]) for p in pairs}672    at_three, at_two, between, below = [], [], [], []673    for p in pairs:674        c = polys[p]675        if not no_root_above(c, 2):676            if evaluate(c, 3) == 0 and no_root_above(c, 3):677                at_three.append(p)678            else:679                between.append(p)680        elif evaluate(c, 2) == 0:681            at_two.append(p)682        else:683            below.append(p)684    check("pairs of spectral radius 3", at_three, [(1, 3), (1, 9), (1, 27)])685    check("pairs in the open interval (2,3)", between, [])686    check("non-shift pairs attaining exactly 2", len(at_two), 20)687    check("the twenty attainers", at_two, [688        (1, 4), (1, 7), (1, 10), (1, 12), (1, 21), (1, 28), (1, 30), (1, 36),689        (3, 4), (3, 7), (3, 10), (3, 28), (4, 9), (4, 27), (7, 9), (7, 27),690        (9, 10), (9, 28), (10, 27), (27, 28)])691    top_lo, top_hi = Fraction(16956207695598, 10 ** 13), Fraction(16956207695599, 10 ** 13)692    next_lo, next_hi = Fraction(16769719158912, 10 ** 13), Fraction(16769719158913, 10 ** 13)693    ties = [p for p in below if not no_root_above(polys[p], Fraction(42, 25))]694    check("pairs with radius above 42/25", len(ties), 25)695    check("the simplest attainer", ties[0], (1, 13))696    check("(12,13) is among them", (12, 13) in ties, True)697    for p in ties:698        check("largest radius below 2, window at %r" % (p,),699              in_window(polys[p], top_lo, top_hi), True)700    second = [p for p in below if p not in ties701              and not no_root_above(polys[p], Fraction(167, 100))]702    check("pairs with radius above 167/100 and below the largest", second,703          [(1, 31), (3, 31), (9, 31), (27, 31)])704    for p in second:705        check("next radius down, window at %r" % (p,),706              in_window(polys[p], next_lo, next_hi), True)707    lo, hi = bisect(polys[(1, 13)], 52)708    check("exact bisection lands inside the quoted window",709          (top_lo <= lo, hi <= top_hi), (True, True))710    for p in [(2, 9), (1, 18)]:711        check("live automaton at %r" % (p,), (trimmed[p], polys[p]), ([[0]], [1, -1]))712    check("(1,4) characteristic polynomial", polys[(1, 4)], [1, -1, -2, 0])713    check("(1,4) live states", len(trimmed[(1, 4)]), 3)714    counts = word_counts(trimmed[(1, 4)], 15)715    check("(1,4) admissible word counts to m=15", counts,716          [(2 ** (m + 2) - (-1) ** m) // 3 for m in range(16)])717    for m in range(2, 16):718        check("(1,4) recurrence residual at m=%d" % m,719              counts[m] - counts[m - 1] - 2 * counts[m - 2], 0)720    returns = return_counts(trimmed[(1, 4)], 15)721    check("(1,4) return counts to m=15", returns,722          [(2 ** (m + 1) + (-1) ** m) // 3 for m in range(16)])723    for (s, t) in [(1, 4), (1, 16), (3, 7), (12, 13)]:724        got = return_counts(trimmed[(s, t)], 8)725        check("return counts against brute force at (%d,%d)" % (s, t), got,726              [pair_brute(s, t, n) for n in range(9)])727    print("829 coprime pairs with max at most 52: radius 3 on three pairs, exactly 2 on twenty, nothing in between")728    free = {p: live(free_automaton(*p)) for p in pairs}729    check("largest free live state set", max(len(free[p]) for p in pairs), 167)730    check("(1,16) free live states", len(free[(1, 16)]), 1)731    for p in pairs:732        if p[0] == 1:733            check("free and gasket-digit automata agree at %r" % (p,),734                  free[p], trimmed[p])735    fpolys = {p: charpoly(free[p]) for p in pairs}736    f3, f2, fmid = [], [], []737    for p in pairs:738        c = fpolys[p]739        if not no_root_above(c, 2):740            (f3 if evaluate(c, 3) == 0 and no_root_above(c, 3) else fmid).append(p)741        elif evaluate(c, 2) == 0:742            f2.append(p)743    check("free pairs of spectral radius 3", f3, at_three)744    check("free pairs in the open interval (2,3)", fmid, [])745    check("free pairs attaining exactly 2", f2, at_two)746    fstrict = [p for p in pairs if p not in f3 and p not in f2]747    flo, fhi = Fraction(18488475886, 10 ** 10), Fraction(18488475887, 10 ** 10)748    ftop = [p for p in fstrict if not no_root_above(fpolys[p], flo)]749    check("free pairs above 1.8488475886", ftop,750          [(4, 13), (4, 39), (12, 13), (13, 36)])751    for p in ftop:752        check("largest free radius below 2, window at %r" % (p,),753              in_window(fpolys[p], flo, fhi), True)754    check("free bisection at (4,13) lands inside the quoted window",755          (flo <= bisect(fpolys[(4, 13)], 52)[0],756           bisect(fpolys[(4, 13)], 52)[1] <= fhi), (True, True))757    theta = [1, -1, 0, -2]758    check("theta is the gasket-digit ceiling below 2", in_window(theta, top_lo, top_hi), True)759    fabove = [p for p in fstrict if not no_root_above(fpolys[p], top_hi)]760    fat = [p for p in fstrict if p not in fabove and not no_root_above(fpolys[p], top_lo)]761    check("free pairs strictly above theta and below 2", len(fabove), 19)762    check("free pairs whose radius is exactly theta", len(fat), 25)763    check("free pairs reaching or beating theta", len(fabove) + len(fat), 44)764    check("the 25 carry the theta factor and the 19 do not",765          (all(divides(theta, fpolys[p]) for p in fat),766           [p for p in fabove if divides(theta, fpolys[p])]), (True, []))767    for p in fat:768        check("radius exactly theta, window at %r" % (p,),769              in_window(fpolys[p], top_lo, top_hi), True)770    check("(4,13) free closed paths at n=60",771          return_counts(free[(4, 13)], 60)[60], 4583352807133551)772    check("(4,13) gasket-digit closed paths at n=60 are smaller",773          return_counts(trimmed[(4, 13)], 60)[60] < 4583352807133551, True)774    print("the free-digit automaton keeps the gap at three, two and the empty "775          "interval, but its largest radius below 2 is 1.8488475886 on four pairs")776777778# WITNESS WEIGHTS779780def carry_automaton(s, t):781    index = {(0, 0): 0}782    order = [(0, 0)]783    edges = []784    i = 0785    while i < len(order):786        a, b = order[i]787        out = []788        for d in range(3):789            u = (s * d + a) % 3790            v = (t * d + b) % 3791            if u < 2 and v < 2:792                nxt = ((s * d + a) // 3, (t * d + b) // 3)793                if nxt not in index:794                    index[nxt] = len(order)795                    order.append(nxt)796                out.append((u, v, index[nxt]))797        edges.append(out)798        i += 1799    return edges800801802def carry_live(edges):803    back = [[] for _ in edges]804    for i, out in enumerate(edges):805        for (u, v, j) in out:806            back[j].append(i)807    seen = {0}808    stack = [0]809    while stack:810        x = stack.pop()811        for y in back[x]:812            if y not in seen:813                seen.add(y)814                stack.append(y)815    keep = sorted(seen)816    place = {v: i for i, v in enumerate(keep)}817    return [[(u, v, place[j]) for (u, v, j) in edges[x] if j in place] for x in keep]818819820def tensor_returns(s, t, n):821    e = carry_live(carry_automaton(s, t))822    k = len(e)823    S = [[j for (u, v, j) in e[i]] for i in range(k)]824    U = [[j for (u, v, j) in e[i] if u] for i in range(k)]825    V = [[j for (u, v, j) in e[i] if v] for i in range(k)]826    W = [[j for (u, v, j) in e[i] if u and v] for i in range(k)]827    X = [[0] * k for _ in range(k)]828    X[0][0] = 1829    out = [1]830    for _ in range(n):831        half = []832        for A in (S, U, V, W):833            Y = []834            for i in range(k):835                acc = [0] * k836                for p in A[i]:837                    xp = X[p]838                    for q in range(k):839                        if xp[q]:840                            acc[q] += xp[q]841                Y.append(acc)842            half.append(Y)843        Ys, Yu, Yv, Yw = half844        Z = [[0] * k for _ in range(k)]845        for j in range(k):846            for i in range(k):847                a = 0848                for q in S[j]:849                    a += Ys[i][q]850                for q in U[j]:851                    a -= Yu[i][q]852                for q in V[j]:853                    a -= Yv[i][q]854                for q in W[j]:855                    a += Yw[i][q]856                Z[i][j] = a857        X = Z858        out.append(X[0][0])859    return out860861862def witness_mass(z1, z2, n):863    e = carry_live(carry_automaton(z1, z2))864    k = len(e)865    adj = [[j for (u, v, j) in e[i] if not (u and v)] for i in range(k)]866    cur = [0] * k867    cur[0] = 1868    out = [1]869    for _ in range(n):870        nxt = [0] * k871        for i in range(k):872            c = cur[i]873            if c:874                for j in adj[i]:875                    nxt[j] += c876        cur = nxt877        out.append(cur[0])878    return [v - 1 for v in out]879880881def in_gasket(x, y, n):882    if x >= 3 ** n or y >= 3 ** n:883        return False884    while x or y:885        a, b = x % 3, y % 3886        if a > 1 or b > 1 or (a and b):887            return False888        x //= 3889        y //= 3890    return True891892893def box_count(s, t, n):894    w = (3 ** n - 1) // (2 * max(s, t))895    c = 0896    for z1 in range(1, w):897        for z2 in range(1, w - z1 + 1):898            if in_gasket(s * z1, s * z2, n) and in_gasket(t * z1, t * z2, n):899                c += 1900    return c901902903def fibre_count(s, t, n):904    binary = set(sum(((u >> i) & 1) * 3 ** i for i in range(n)) for u in range(2 ** n))905    return sum(1 for u in range(1, 3 ** n // t + 1) if s * u in binary and t * u in binary)906907908def no_adjacent(n):909    out = []910    for m in range(1, 3 ** (n - 1)):911        x, prev, ok = m, 0, True912        while x:913            d = x % 3914            if d > 1 or (d and prev):915                ok = False916                break917            prev, x = d, x // 3918        if ok:919            out.append(m)920    return out921922923def three_power_ratio(a, b):924    if a > b:925        a, b = b, a926    if b % a:927        return False928    q = b // a929    while q % 3 == 0:930        q //= 3931    return q == 1932933934def residual_layers(n):935    by = {}936    for (x, y) in points(GASKET, n):937        if x and y:938            g = gcd(x, y)939            by.setdefault((x // g, y // g), []).append(g)940    total = 0941    weight = {}942    pair = {}943    for r, gs in by.items():944        if len(gs) < 2:945            continue946        rs = r[0] + r[1]947        for a in gs:948            for b in gs:949                if a == b:950                    continue951                d = gcd(a, b)952                if is_shift_pair(a // d, b // d):953                    continue954                total += 1955                weight[d * rs] = weight.get(d * rs, 0) + 1956                k = (a // d, b // d)957                if k[0] > k[1]:958                    k = (k[1], k[0])959                pair[k] = pair.get(k, 0) + 1960    return total, weight, pair961962963def ceiling_families():964    bin3 = [sum(((u >> i) & 1) * 3 ** i for i in range(7)) for u in range(1, 128)]965    nad = []966    for m in range(1, 3 ** 7):967        x, prev, ok = m, 0, True968        while x:969            d = x % 3970            if d > 1 or (d and prev):971                ok = False972                break973            prev, x = d, x // 3974        if ok:975            nad.append(m)976    return [977        [(a, b) for a in bin3 for b in bin3 if a < b],978        [(a, b) for a in nad for b in nad if a < b],979        [(1, t) for t in range(2, 3000)],980        [(a, a + 1) for a in range(1, 1500)],981        [(a, 3 * a - 1) for a in range(1, 1200)],982        [(a, 3 * a + 1) for a in range(1, 1200)],983    ]984985986def second_moment(n):987    keys = []988    for (x, y) in points(GASKET, n):989        if x and y:990            g = gcd(x, y)991            keys.append((x // g) * 3 ** n + y // g)992    keys.sort()993    total = 0994    run = 1995    for i in range(1, len(keys)):996        if keys[i] == keys[i - 1]:997            run += 1998        else:999            total += run * run1000            run = 11001    return total + run * run100210031004# THE PROVED CEILING100510061007def direction_automaton(a, b):1008    index = {0: 0}1009    order = [0]1010    edges = []1011    i = 01012    while i < len(order):1013        c = order[i]1014        out = []1015        for eps in (0, b, -a):1016            if (c + eps) % 3 == 0:1017                nxt = (c + eps) // 31018                if nxt not in index:1019                    index[nxt] = len(order)1020                    order.append(nxt)1021                out.append((eps, index[nxt]))1022        edges.append(out)1023        i += 11024    return order, edges102510261027def direction_live(a, b):1028    order, edges = direction_automaton(a, b)1029    back = [[] for _ in edges]1030    for i, out in enumerate(edges):1031        for (eps, j) in out:1032            back[j].append(i)1033    seen = {0}1034    stack = [0]1035    while stack:1036        x = stack.pop()1037        for y in back[x]:1038            if y not in seen:1039                seen.add(y)1040                stack.append(y)1041    keep = sorted(seen)1042    place = {v: i for i, v in enumerate(keep)}1043    return ([order[x] for x in keep],1044            [[(eps, place[j]) for (eps, j) in edges[x] if j in place]1045             for x in keep])104610471048def direction_mass(a, b, n):1049    order, edges = direction_live(a, b)1050    cur = [0] * len(edges)1051    cur[0] = 11052    out = [1]1053    for _ in range(n):1054        nxt = [0] * len(edges)1055        for i in range(len(edges)):1056            if cur[i]:1057                for (eps, j) in edges[i]:1058                    nxt[j] += cur[i]1059        cur = nxt1060        out.append(cur[0])1061    return [v - 1 for v in out]106210631064def state_profile(order, edges, n):1065    tab = [[1 if order[i] == 0 else 0 for i in range(len(edges))]]1066    for _ in range(n):1067        prev = tab[-1]1068        tab.append([sum(prev[j] for (eps, j) in edges[i])1069                    for i in range(len(edges))])1070    return [max(r) for r in tab]107110721073def valuation3(x):1074    k = 01075    while x % 3 == 0:1076        x //= 31077        k += 11078    return k107910801081def shift_ray(a, b):1082    lo, hi = min(a, b), max(a, b)1083    if lo != 1:1084        return False1085    while hi % 3 == 0:1086        hi //= 31087    return hi == 1108810891090def certificate_holds(a, b, den, alpha, beta):1091    order, edges = direction_live(a, b)1092    if len(order) != len(alpha) or len(order) != len(beta):1093        return False1094    if alpha[0] != den or beta[0] != 0:1095        return False1096    for i in range(len(edges)):1097        if alpha[i] < 0:1098            return False1099        if sum(alpha[j] for (eps, j) in edges[i]) > alpha[i] + beta[i]:1100            return False1101        if sum(beta[j] for (eps, j) in edges[i]) > alpha[i]:1102            return False1103    return True110411051106def binary_multiples(w, n):1107    cur = [0] * w1108    cur[0] = 11109    for i in range(n):1110        p = pow(3, i, w)1111        nxt = list(cur)1112        for r in range(w):1113            if cur[r]:1114                nxt[(r + p) % w] += cur[r]1115        cur = nxt1116    return cur[0] - 1111711181119def qnorm(p, q, d):1120    if d < 0:1121        p, q, d = -p, -q, -d1122    g = gcd(gcd(abs(p), abs(q)), d)1123    if g > 1:1124        p //= g1125        q //= g1126        d //= g1127    return (p, q, d)112811291130def qadd(x, y):1131    return qnorm(x[0] * y[2] + y[0] * x[2],1132                 x[1] * y[2] + y[1] * x[2], x[2] * y[2])113311341135def qsub(x, y):1136    return qnorm(x[0] * y[2] - y[0] * x[2],1137                 x[1] * y[2] - y[1] * x[2], x[2] * y[2])113811391140def qmul(x, y):1141    return qnorm(x[0] * y[0] + x[1] * y[1],1142                 x[0] * y[1] + x[1] * y[0] + x[1] * y[1], x[2] * y[2])114311441145def qinv(x):1146    p, q, d = x1147    return qnorm(d * (p + q), -d * q, p * p + p * q - q * q)114811491150def qsgn(x):1151    p, q, d = x1152    hi, lo = 2 * p + q, q1153    if hi >= 0 and lo >= 0:1154        return 0 if hi == 0 and lo == 0 else 11155    if hi <= 0 and lo <= 0:1156        return 0 if hi == 0 and lo == 0 else -11157    s = hi * hi - 5 * lo * lo1158    if hi > 0:1159        return 1 if s > 0 else (0 if s == 0 else -1)1160    return -1 if s > 0 else (0 if s == 0 else 1)116111621163QZERO = (0, 0, 1)1164QONE = (1, 0, 1)1165QPHI = (0, 1, 1)1166QINVPHI = (-1, 1, 1)1167QINVPHI2 = (2, -1, 1)116811691170def golden_potential(order, edges):1171    n = len(order)1172    if n == 1:1173        return None1174    m = n - 11175    rows = [[QZERO] * (m + 1) for _ in range(m)]1176    for r in range(m):1177        rows[r][r] = QPHI1178        for (eps, j) in edges[r + 1]:1179            if j == 0:1180                rows[r][m] = qadd(rows[r][m], QONE)1181            else:1182                rows[r][j - 1] = qsub(rows[r][j - 1], QONE)1183    for col in range(m):1184        piv = None1185        for r in range(col, m):1186            if qsgn(rows[r][col]):1187                piv = r1188                break1189        if piv is None:1190            return "singular"1191        rows[col], rows[piv] = rows[piv], rows[col]1192        scale = qinv(rows[col][col])1193        rows[col] = [qmul(v, scale) if qsgn(v) else QZERO for v in rows[col]]1194        for r in range(m):1195            if r != col and qsgn(rows[r][col]):1196                f = rows[r][col]1197                rows[r] = [qsub(rows[r][k], qmul(f, rows[col][k]))1198                           if qsgn(rows[col][k]) else rows[r][k]1199                           for k in range(m + 1)]1200    return [QONE] + [rows[r][m] for r in range(m)]120112021203def potential_value(u, edges):1204    tot = QZERO1205    for (eps, j) in edges[0]:1206        if j:1207            tot = qadd(tot, u[j])1208    return tot120912101211def potential_valid(u, edges):1212    if u[0] != QONE or any(qsgn(v) <= 0 for v in u):1213        return False1214    for i in range(1, len(edges)):1215        s = QZERO1216        for (eps, j) in edges[i]:1217            s = qadd(s, u[j])1218        if qsgn(qsub(qmul(QPHI, u[i]), s)) < 0:1219            return False1220    return True122112221223CERTIFICATES = {1224    (1, 90): (18, [18, 0, 5, 6, 9, 2, 10, 4, 8],1225              [0, 18, 4, -4, 5, 6, 8, 1, 2]),1226    (4, 117): (40, [40, 0, 5, 14, 6, 13, 11, 27, 17, 1],1227               [0, 40, 1, -1, 5, 14, 6, 13, 11, 4]),1228    (9, 73): (381,1229              [381, 0, 46, 49, 78, 17, 101, 23, 66, 163, 39, 83, 232, 32,1230               62, 149],1231              [0, 381, 32, -32, 46, 49, 62, 16, 17, 101, 23, 66, 149, 14,1232               39, 83]),1233    (9, 82): (18, [18, 0, 5, 6, 9, 2, 10, 4, 8],1234              [0, 18, 4, -4, 5, 6, 8, 1, 2]),1235    (9, 235): (2013,1236               [2013, 0, 16, 66, 27, 55, 43, 121, 70, 176, 113, 297, 183,1237                473, 296, 770, 479, 1243, 775, 11],1238               [0, 2013, 11, -11, 16, 66, 27, 55, 43, 121, 70, 176, 113,1239                297, 183, 473, 296, 770, 479, 5]),1240    (10, 81): (18, [18, 0, 5, 6, 9, 2, 10, 4, 8],1241               [0, 18, 4, -4, 5, 6, 8, 1, 2]),1242    (13, 108): (40, [40, 0, 5, 14, 6, 13, 11, 27, 17, 1],1243                [0, 40, 1, -1, 5, 14, 6, 13, 11, 4]),1244    (27, 217): (2013,1245                [2013, 0, 16, 66, 27, 55, 43, 121, 70, 176, 113, 297, 183,1246                 473, 296, 770, 479, 1243, 775, 11],1247                [0, 2013, 11, -11, 16, 66, 27, 55, 43, 121, 70, 176, 113,1248                 297, 183, 473, 296, 770, 479, 5]),1249    (27, 226): (34, [34, 0, 0, 1, 1, 0, 1, 1, 2, 1, 3, 5, 8, 13, 20, 1],1250                [0, 34, 1, -1, 0, 1, 1, 0, 1, 1, 2, 3, 5, 8, 14, -1]),1251}125212531254def census_ceiling():1255    tested = mismatch = 01256    for a in range(1, 40):1257        for b in range(1, 40):1258            if gcd(a, b) != 1:1259                continue1260            tested += 11261            if direction_mass(a, b, 20) != witness_mass(a, b, 20):1262                mismatch += 11263    check("direction carry automaton against the gasket-digit automaton",1264          (tested, mismatch), (947, 0))1265    box = occupied = bounded = single = k1 = nodouble = 01266    twoplus = anomaly = states = fails = 01267    residue = certified = passing = 01268    best = QZERO1269    attain = []1270    over = []1271    exceptions = []1272    for z1 in range(1, 121):1273        for z2 in range(z1, 241):1274            if gcd(z1, z2) != 1:1275                continue1276            box += 11277            d = (z1 % 3 == 0) + (z2 % 3 == 0) + ((z1 + z2) % 3 == 0)1278            if d > 1:1279                twoplus += 11280            if z1 % 3 and z2 % 3:1281                residue += 11282            order, edges = direction_live(z1, z2)1283            if len(order) == 1:1284                continue1285            if d == 0:1286                anomaly += 11287            occupied += 11288            states = max(states, len(order))1289            if all(-z1 <= 2 * c <= z2 for c in order):1290                bounded += 11291            degs = [len(o) for o in edges]1292            classes = set(order[i] % 3 for i in range(len(order))1293                          if degs[i] == 2)1294            if max(degs) <= 2 and len(classes) <= 1:1295                single += 11296            q = z1 if z1 % 3 == 0 else (z2 if z2 % 3 == 0 else z1 + z2)1297            if valuation3(q) == 1:1298                k1 += 11299            if not any(degs[i] == 21300                       and all(degs[j] == 2 for (eps, j) in edges[i])1301                       for i in range(len(edges))):1302                nodouble += 11303            else:1304                exceptions.append((z1, z2))1305            h = state_profile(order, edges, 20)1306            if any(h[m] > h[m - 1] + h[m - 2] for m in range(2, 21)):1307                fails += 11308            u = golden_potential(order, edges)1309            if potential_valid(u, edges):1310                certified += 11311            tot = potential_value(u, edges)1312            if qsgn(qsub(QINVPHI2, tot)) >= 0:1313                passing += 11314                if qsgn(qsub(tot, best)) > 0:1315                    best, attain = tot, [(z1, z2)]1316                elif tot == best:1317                    attain.append((z1, z2))1318            else:1319                over.append((z1, z2, tot))1320    check("carry states of a direction lie in [-a/2, b/2]",1321          (box, occupied, states, bounded), (13158, 218, 37, 218))1322    check("out-degree at most two, branch states in one class mod 3",1323          single, 218)1324    check("at most one of z1, z2, z1+z2 is divisible by three, and no "1325          "occupied direction has none of them divisible",1326          (twoplus, anomaly), (0, 0))1327    check("three divides a coordinate of every occupied direction",1328          (residue, box - residue), (6566, 6592))1329    check("the golden potential certifies the box away from the shift rays",1330          (certified, passing, [(z1, z2) for (z1, z2, tot) in over],1331           sorted(set(tot for (z1, z2, tot) in over))),1332          (218, 214, [(1, 3), (1, 9), (1, 27), (1, 81)], [QINVPHI]))1333    check("the golden potential is largest on the supergolden directions",1334          (best, attain), (QINVPHI2, [(1, 12), (3, 10), (4, 9)]))1335    check("directions settled by the branch argument over the box",1336          (k1, nodouble, exceptions),1337          (107, 206, [(1, 9), (1, 27), (1, 81), (1, 90), (4, 117), (9, 73),1338                      (9, 82), (9, 235), (10, 81), (13, 108), (27, 217),1339                      (27, 226)]))1340    good = sorted(k for k, v in CERTIFICATES.items()1341                  if certificate_holds(k[0], k[1], v[0], v[1], v[2]))1342    check("Fibonacci certificates for the nine non-shift exceptions",1343          good, [(1, 90), (4, 117), (9, 73), (9, 82), (9, 235), (10, 81),1344                 (13, 108), (27, 217), (27, 226)])1345    identity = all(fibonacci(p + 2) * fibonacci(q + 2)1346                   == fibonacci(p + q + 3) - fibonacci(p + 1) * fibonacci(q + 1)1347                   for p in range(60) for q in range(60))1348    cases = 01349    over = 01350    strict = 01351    for j in range(1, 14):1352        for m in range(1, 46):1353            cases += 11354            if shift_product(j, m) > fibonacci(m + 1):1355                over += 11356            if j >= 2 and m >= 2 and shift_product(j, m) >= fibonacci(m + 1):1357                strict += 11358    check("Fibonacci product identity and the shift-ray ceiling",1359          (identity, cases, over, strict), (True, 585, 0, 0))1360    tested = breaches = 01361    for z1 in range(1, 31):1362        for z2 in range(z1, 61):1363            if gcd(z1, z2) != 1:1364                continue1365            tested += 11366            if direction_mass(z1, z2, 12)[12] > binary_multiples(z1 + z2, 12):1367                breaches += 11368    check("ray mass is at most the count of binary multiples of the weight",1369          (tested, breaches), (829, 0))1370    check("binary multiples of the weight outgrow the ceiling",1371          ([binary_multiples(w, 24) for w in (4, 10, 28, 82)],1372           fibonacci(25) - 1),1373          ([4196351, 1683971, 613817, 228519], 75024))1374    seen = set()1375    multiplicity = 01376    for family in ceiling_families():1377        for (a, b) in family:1378            if gcd(a, b) != 1:1379                continue1380            multiplicity += 11381            seen.add((a, b))1382    inside = set((a, b) for (a, b) in seen if a <= 120 and b <= 240)1383    zero = case2 = case3 = enumonly = 01384    outside = outpass = 01385    outover = []1386    for (a, b) in sorted(seen - inside):1387        order, edges = direction_live(a, b)1388        if len(order) == 1:1389            zero += 11390            continue1391        outside += 11392        u = golden_potential(order, edges)1393        if potential_valid(u, edges) and \1394                qsgn(qsub(QINVPHI2, potential_value(u, edges))) >= 0:1395            outpass += 11396        else:1397            outover.append((a, b, potential_value(u, edges)))1398        degs = [len(o) for o in edges]1399        if not any(degs[i] == 21400                   and all(degs[j] == 2 for (eps, j) in edges[i])1401                   for i in range(len(edges))):1402            case2 += 11403        elif shift_ray(a, b):1404            case3 += 11405        else:1406            enumonly += 11407    check("the six adversarial families overlap, and their union splits",1408          (multiplicity, len(seen), len(inside), len(seen) - len(inside),1409           zero, case2, case3, enumonly),1410          (11369, 10862, 717, 10145, 9498, 608, 3, 36))1411    check("the golden potential outside the box, shift rays excepted",1412          (outside, outpass, outover),1413          (647, 644, [(1, 243, QINVPHI), (1, 729, QINVPHI),1414                      (1, 2187, QINVPHI)]))1415    profile = state_profile(*direction_live(1, 9), 12)1416    check("the state maximum breaks the Fibonacci recursion",1417          (profile[:7], profile[4] > profile[3] + profile[2], fails),1418          ([1, 1, 1, 2, 4, 6, 9], True, 8))1419    print("golden ceiling: proved for every direction of the box")142014211422def split_three(a, b):1423    q, p = (a, b) if a % 3 == 0 else (b, a)1424    k = 01425    q1 = q1426    while q1 % 3 == 0:1427        q1 //= 31428        k += 11429    return p, k, q1143014311432def phi_power(e):1433    r = QONE1434    for _ in range(abs(e)):1435        r = qmul(r, QPHI if e > 0 else QINVPHI)1436    return r143714381439def first_returns(edges, n):1440    cur = [0] * len(edges)1441    cur[0] = 11442    f = [0] * (n + 1)1443    for m in range(1, n + 1):1444        nxt = [0] * len(edges)1445        for i in range(len(edges)):1446            if cur[i]:1447                for (eps, j) in edges[i]:1448                    nxt[j] += cur[i]1449        f[m] = nxt[0]1450        nxt[0] = 01451        cur = nxt1452    return f145314541455def degree_potential(edges):1456    return [QONE if len(o) == 2 else QINVPHI for o in edges]145714581459def super_solution(edges, pi):1460    for i in range(1, len(edges)):1461        s = QZERO1462        for (eps, j) in edges[i]:1463            s = qadd(s, pi[j])1464        if qsgn(qsub(qmul(QPHI, pi[i]), s)) < 0:1465            return False1466    return True146714681469def swept_bound(edges, d):1470    pi = degree_potential(edges)1471    pi[0] = QONE1472    for _ in range(d):1473        nxt = [QONE] + [QZERO] * (len(edges) - 1)1474        for i in range(1, len(edges)):1475            s = QZERO1476            for (eps, j) in edges[i]:1477                s = qadd(s, pi[j])1478            nxt[i] = qmul(QINVPHI, s)1479        pi = nxt1480    tot = QZERO1481    for (eps, j) in edges[0]:1482        if j:1483            tot = qadd(tot, pi[j])1484    return tot148514861487def burst_floor(a, b, k, q1):1488    sign = 1 if b % 3 == 0 else -11489    out = []1490    for m in range(1, 3 ** k):1491        if m % 3 != 1:1492            continue1493        x, ok = m, True1494        while x:1495            if x % 3 > 1:1496                ok = False1497                break1498            x //= 31499        if ok:1500            out.append(sign * q1 * m)1501    return out150215031504def census_degree():1505    seen = set()1506    for z1 in range(1, 121):1507        for z2 in range(z1, 241):1508            if gcd(z1, z2) == 1:1509                seen.add((z1, z2))1510    for family in ceiling_families():1511        for (a, b) in family:1512            if gcd(a, b) == 1:1513                seen.add((a, b))1514    triple = eligible = occupied = resbad = burst = 01515    burstbad = 01516    nodouble = valid = doublevalid = settled = 01517    depths = {}1518    missed = []1519    k1 = k1class = 01520    k1bad = []1521    quantbad = []1522    attain = []1523    for (a, b) in sorted(seen):1524        p, k, q1 = split_three(a, b)1525        if a % 3 == 0 or b % 3 == 0:1526            triple += 11527            if (q1 - p) % 3 == 0:1528                eligible += 11529        order, edges = direction_live(a, b)1530        if len(order) == 1:1531            continue1532        occupied += 11533        if (q1 - p) % 3:1534            resbad += 11535        f = first_returns(edges, max(k, 2))1536        if any(f[j] for j in range(2, k + 1)):1537            burst += 11538        degs = [len(o) for o in edges]1539        double = any(degs[i] == 21540                     and all(degs[j] == 2 for (eps, j) in edges[i])1541                     for i in range(len(edges)))1542        if not double:1543            nodouble += 11544        if super_solution(edges, degree_potential(edges)):1545            valid += 11546            if double:1547                doublevalid += 11548            hit = None1549            for d in range(25):1550                if qsgn(qsub(QINVPHI2, swept_bound(edges, d))) >= 0:1551                    hit = d1552                    break1553            if hit is None:1554                missed.append((a, b))1555            else:1556                settled += 11557                depths[hit] = depths.get(hit, 0) + 11558        else:1559            missed.append((a, b))1560        u = golden_potential(order, edges)1561        tot = potential_value(u, edges)1562        place = {c: i for i, c in enumerate(order)}1563        floor = burst_floor(a, b, k, q1)1564        rung = QZERO1565        for c in floor:1566            if c in place:1567                rung = qadd(rung, u[place[c]])1568        if len(floor) != 2 ** (k - 1) or \1569                qmul(phi_power(-(k - 1)), rung) != tot:1570            burstbad += 11571        if 2 * p <= 3 ** k * q1:1572            c = p if b % 3 == 0 else -p1573            if c not in place or u[place[c]] != QINVPHI:1574                burstbad += 11575        if k != 1:1576            continue1577        k1 += 11578        if q1 == p:1579            continue1580        t = valuation3(q1 - p)1581        bound = qmul(QINVPHI, qsub(QONE, phi_power(-max(t, 2))))1582        if qsgn(qsub(bound, tot)) < 0:1583            quantbad.append((a, b, t))1584        if tot == bound:1585            attain.append((a, b))1586        if t <= 2:1587            k1class += 11588            if qsgn(qsub(QINVPHI2, tot)) < 0:1589                k1bad.append((a, b, t))1590    check("occupancy needs the residue match q1 = p mod 3",1591          (len(seen), triple, eligible, occupied, resbad),1592          (23303, 11691, 7103, 865, 0))1593    check("no first return has length between two and v3(q)",1594          (occupied, burst), (865, 0))1595    check("the burst identity and the value at the near predecessor",1596          (occupied, burstbad), (865, 0))1597    check("the degree potential is a super-solution beyond the branch case",1598          (nodouble, valid, doublevalid), (814, 851, 37))1599    check("the swept degree potential settles all but sixteen directions",1600          (settled, sorted(depths.items()),1601           [z for z in missed if shift_ray(*z)],1602           [z for z in missed if not shift_ray(*z)]),1603          (849, [(1, 760), (3, 48), (4, 31), (5, 7), (6, 3)],1604           [(1, 3 ** j) for j in range(1, 8)],1605           [(1, 756), (1, 2196), (1, 2214), (1, 2268), (1, 2430),1606            (13, 1080), (27, 730), (28, 729), (40, 1053)]))1607    check("the golden partition bound at v3(q) = 1",1608          (k1, k1class, k1bad, quantbad, attain),1609          (360, 261, [], [], [(1, 12), (3, 10)]))1610    short = []1611    for a in range(1, 130):1612        for b in range(1, 800):1613            if gcd(a, b) != 1 or (a % 3 and b % 3):1614                continue1615            order, edges = direction_live(a, b)1616            if len(order) == 1:1617                continue1618            f = first_returns(edges, 4)1619            if f[2] or f[3]:1620                short.append((min(a, b), max(a, b), f[2], f[3]))1621    check("the first returns of length two and three are classified",1622          sorted(set(short)),1623          [(1, 3, 1, 0), (1, 9, 0, 1), (1, 12, 0, 1), (3, 10, 0, 1),1624           (4, 9, 0, 1)])1625    print("degree potential: the golden partition bound on an infinite class")162616271628def census_witness():1629    n = 91630    bad = 01631    tested = 01632    for s in range(1, 40):1633        for t in range(s + 1, 40):1634            if gcd(s, t) != 1:1635                continue1636            tested += 11637            if free_returns(free_automaton(s, t), n) != tensor_returns(s, t, n):1638                bad += 11639    check("tensor square against B(s,t) return counts to n=9", (tested, bad), (473, 0))1640    check("carry states against reachable B states",1641          [(len(carry_live(carry_automaton(s, t))), len(free_automaton(s, t)))1642           for (s, t) in ((365, 1094), (41, 122), (25, 52), (31, 40))],1643          [(729, 26931), (81, 835), (38, 393), (35, 354)])1644    bad = 01645    tested = 01646    for s in range(1, 30):1647        for t in range(s + 1, 60):1648            if gcd(s, t) != 1:1649                continue1650            tested += 11651            want = free_returns(free_automaton(s, t), n)[n] - 1 - 2 * fibre_count(s, t, n)1652            if want != box_count(s, t, n):1653                bad += 11654    check("box construction against B(s,t) at n=9", (tested, bad), (812, 0))1655    got = []1656    for (s, t) in ((365, 1094), (41, 122), (122, 123), (1, 2460), (2431, 2458)):1657        for m in (9, 12):1658            got.append((tensor_returns(s, t, m)[m] - 1 - 2 * fibre_count(s, t, m),1659                        box_count(s, t, m)))1660    check("box against tensor at large multipliers",1661          ([x for x in got if x[0] != x[1]], [g[0] for g in got]),1662          ([], [2, 180, 50, 172, 50, 172, 2, 12, 2, 8]))1663    R, weight, pair = {}, {}, {}1664    for m in range(4, 14):1665        R[m], weight[m], pair[m] = residual_layers(m)1666    check("R(n) for n=4..13", [R[m] for m in range(4, 14)],1667          [20, 88, 432, 1624, 5512, 15896, 46064, 124928, 335704, 863848])1668    scaled = bad = 01669    for m in range(5, 14):1670        for w, c in weight[m].items():1671            if w % 3 == 0:1672                scaled += 11673                if weight[m - 1].get(w // 3, 0) != c:1674                    bad += 11675    check("weight layers scale as R_3w(n) = R_w(n-1)", (scaled, bad), (1869, 0))1676    check("no witness of weight below four",1677          sorted(set(min(weight[m]) for m in range(4, 14))), [4])1678    check("largest multiplier is floor(3^n/8) for n = 4..13",1679          [max(k[1] for k in pair[m]) for m in range(4, 14)],1680          [3 ** m // 8 for m in range(4, 14)])1681    fours = []1682    sets = []1683    for m in range(4, 13):1684        F = no_adjacent(m)1685        sets.append(len(F))1686        fours.append(2 * sum(1 for a in F for b in F if a != b1687                             and gcd(a, b) == 1 and not three_power_ratio(a, b)))1688    check("R_4(n) from the no-adjacent-ones set, n=4..12",1689          (fours, [weight[m].get(4, 0) for m in range(4, 13)]),1690          ([12, 36, 108, 336, 988, 2596, 6672, 17480, 45720],1691           [12, 36, 108, 336, 988, 2596, 6672, 17480, 45720]))1692    check("size of the no-adjacent-ones set is Fibonacci",1693          sets, [fibonacci(m + 1) - 1 for m in range(4, 13)])1694    check("weight-four orbit at n=13", sum(fours) + weight[13].get(4, 0), 194096)1695    tops = []1696    for m in range(6, 14):1697        thr = (3 ** m - 1) // 101698        big = [v for k, v in pair[m].items() if k[1] > thr]1699        tops.append((len(big), sorted(set(big))))1700    check("pairs above (3^n-1)/10 each contribute exactly four",1701          tops, [(18, [4]), (57, [4]), (163, [4]), (402, [4]), (1019, [4]),1702                 (2702, [4]), (7060, [4]), (18607, [4])])1703    N = 401704    ref = [fibonacci(m + 1) - 1 for m in range(N + 1)]1705    check("heaviest ray mass is Fibonacci", witness_mass(1, 3, N), ref)1706    tested = breaches = 01707    ties = []1708    for z1 in range(1, 121):1709        for z2 in range(z1, 241):1710            if gcd(z1, z2) != 1:1711                continue1712            tested += 11713            m = witness_mass(z1, z2, N)1714            for j in range(1, N + 1):1715                if m[j] > ref[j]:1716                    breaches += 11717            if m[N] == ref[N]:1718                ties.append((z1, z2))1719    check("golden ceiling on ray mass over the box, n <= 40",1720          (tested, breaches, ties), (13158, 0, [(1, 3)]))1721    N = 451722    ref = [fibonacci(m + 1) - 1 for m in range(N + 1)]1723    sizes = []1724    breaches = 01725    for family in ceiling_families():1726        seen = 01727        for (a, b) in family:1728            if gcd(a, b) != 1:1729                continue1730            seen += 11731            m = witness_mass(a, b, N)1732            for j in range(1, N + 1):1733                if m[j] > ref[j]:1734                    breaches += 11735        sizes.append(seen)1736    check("golden ceiling on six adversarial families, n <= 45",1737          (sizes, sum(sizes), breaches),1738          ([4221, 253, 2998, 1499, 1199, 1199], 11369, 0))1739    moments = []1740    for m in (13, 14):1741        E = second_moment(m)1742        D = 3 ** m - 2 ** (m + 1) + 11743        S = 3 ** m - 4 * 2 ** m + 2 * m + 31744        moments.append((E, E - D - S))1745    check("second moment and residual at n = 13 and 14",1746          moments, [(4003372, 863848), (11679626, 2211960)])1747    print("witness weights: weight four is Fibonacci and the ceiling holds")174817491750def free_returns(edges, n):1751    count = [0] * len(edges)1752    count[0] = 11753    out = []1754    for _ in range(n + 1):1755        out.append(count[0])1756        nxt = [0] * len(edges)1757        for i, outs in enumerate(edges):1758            for e in outs:1759                nxt[e] += count[i]1760        count = nxt1761    return out176217631764def main():1765    census_designs()1766    census_cross()1767    census_diagonality()1768    census_masses()1769    census_recurrences()1770    census_shift()1771    census_rays()1772    census_pairs()1773    census_spectrum()1774    census_witness()1775    census_ceiling()1776    census_degree()1777    print("all green")177817791780if __name__ == "__main__":1781    main()