witness.py

1from math import gcd23GASKET = ((0, 0), (1, 0), (0, 1))4FIB = [0, 1, 1]5while len(FIB) < 200:6    FIB.append(FIB[-1] + FIB[-2])789def die(name, got, want):10    raise SystemExit("FAIL %s\n  got  %r\n  want %r" % (name, got, want))111213def check(name, got, want):14    if got != want:15        die(name, got, want)16    print("ok   %s" % name)171819# GASKET POINTS2021def points(n):22    pts = [(0, 0)]23    for _ in range(n):24        nxt = []25        for (x, y) in pts:26            for (dx, dy) in GASKET:27                nxt.append((3 * x + dx, 3 * y + dy))28        pts = nxt29    return pts303132def in_gasket(x, y, n):33    if x >= 3 ** n or y >= 3 ** n:34        return False35    while x or y:36        a, b = x % 3, y % 337        if a > 1 or b > 1 or (a and b):38            return False39        x //= 340        y //= 341    return True424344def is_shift_pair(s, t):45    if min(s, t) != 1:46        return False47    u = max(s, t)48    while u % 3 == 0:49        u //= 350    return u == 1515253# CARRY AUTOMATON5455def carry(s, t):56    index = {(0, 0): 0}57    order = [(0, 0)]58    edges = []59    i = 060    while i < len(order):61        a, b = order[i]62        out = []63        for d in range(3):64            u = (s * d + a) % 365            v = (t * d + b) % 366            if u < 2 and v < 2:67                nxt = ((s * d + a) // 3, (t * d + b) // 3)68                if nxt not in index:69                    index[nxt] = len(order)70                    order.append(nxt)71                out.append((u, v, index[nxt]))72        edges.append(out)73        i += 174    return edges757677def carry_live(edges):78    back = [[] for _ in edges]79    for i, out in enumerate(edges):80        for (u, v, j) in out:81            back[j].append(i)82    seen = {0}83    stack = [0]84    while stack:85        x = stack.pop()86        for y in back[x]:87            if y not in seen:88                seen.add(y)89                stack.append(y)90    keep = sorted(seen)91    place = {v: i for i, v in enumerate(keep)}92    return [[(u, v, place[j]) for (u, v, j) in edges[x] if j in place] for x in keep]939495def ray_mass(z1, z2, n):96    e = carry_live(carry(z1, z2))97    k = len(e)98    adj = [[j for (u, v, j) in e[i] if not (u and v)] for i in range(k)]99    cur = [0] * k100    cur[0] = 1101    out = [1]102    for _ in range(n):103        nxt = [0] * k104        for i in range(k):105            c = cur[i]106            if c:107                for j in adj[i]:108                    nxt[j] += c109        cur = nxt110        out.append(cur[0])111    return [v - 1 for v in out]112113114# FREE AUTOMATON AS A CONSTRAINED TENSOR SQUARE115116def free_automaton(s, t):117    index = {(0, 0, 0, 0): 0}118    order = [(0, 0, 0, 0)]119    edges = []120    i = 0121    while i < len(order):122        a1, a2, b1, b2 = order[i]123        out = []124        for dx in range(3):125            for dy in range(3):126                u = ((s * dx + a1) % 3, (s * dy + a2) % 3)127                v = ((t * dx + b1) % 3, (t * dy + b2) % 3)128                if u in GASKET and v in GASKET:129                    nxt = ((s * dx + a1) // 3, (s * dy + a2) // 3,130                           (t * dx + b1) // 3, (t * dy + b2) // 3)131                    if nxt not in index:132                        index[nxt] = len(order)133                        order.append(nxt)134                    out.append(index[nxt])135        edges.append(out)136        i += 1137    return edges138139140def free_returns(edges, n):141    count = [0] * len(edges)142    count[0] = 1143    out = []144    for _ in range(n + 1):145        out.append(count[0])146        nxt = [0] * len(edges)147        for i, outs in enumerate(edges):148            for e in outs:149                nxt[e] += count[i]150        count = nxt151    return out152153154def tensor_returns(s, t, n):155    e = carry_live(carry(s, t))156    k = len(e)157    S = [[j for (u, v, j) in e[i]] for i in range(k)]158    U = [[j for (u, v, j) in e[i] if u] for i in range(k)]159    V = [[j for (u, v, j) in e[i] if v] for i in range(k)]160    W = [[j for (u, v, j) in e[i] if u and v] for i in range(k)]161    X = [[0] * k for _ in range(k)]162    X[0][0] = 1163    out = [1]164    for _ in range(n):165        half = []166        for A in (S, U, V, W):167            Y = []168            for i in range(k):169                acc = [0] * k170                for p in A[i]:171                    xp = X[p]172                    for q in range(k):173                        if xp[q]:174                            acc[q] += xp[q]175                Y.append(acc)176            half.append(Y)177        Ys, Yu, Yv, Yw = half178        Z = [[0] * k for _ in range(k)]179        for j in range(k):180            for i in range(k):181                a = 0182                for q in S[j]:183                    a += Ys[i][q]184                for q in U[j]:185                    a -= Yu[i][q]186                for q in V[j]:187                    a -= Yv[i][q]188                for q in W[j]:189                    a += Yw[i][q]190                Z[i][j] = a191        X = Z192        out.append(X[0][0])193    return out194195196# BOX CONSTRUCTION197198def box(s, t, n):199    W = (3 ** n - 1) // (2 * max(s, t))200    c = 0201    for z1 in range(1, W):202        for z2 in range(1, W - z1 + 1):203            if in_gasket(s * z1, s * z2, n) and in_gasket(t * z1, t * z2, n):204                c += 1205    return c206207208def fibre(s, t, n):209    binary = set(sum(((u >> i) & 1) * 3 ** i for i in range(n)) for u in range(2 ** n))210    return sum(1 for u in range(1, 3 ** n // t + 1) if s * u in binary and t * u in binary)211212213# CENSUS214215def ray_census(n):216    by = {}217    for (x, y) in points(n):218        if x and y:219            g = gcd(x, y)220            by.setdefault((x // g, y // g), []).append(g)221    return by222223224def residual(n):225    by = ray_census(n)226    R = 0227    weight = {}228    pair = {}229    for r, gs in by.items():230        if len(gs) < 2:231            continue232        rs = r[0] + r[1]233        for a in gs:234            for b in gs:235                if a == b:236                    continue237                d = gcd(a, b)238                if is_shift_pair(a // d, b // d):239                    continue240                R += 1241                weight[d * rs] = weight.get(d * rs, 0) + 1242                k = (a // d, b // d)243                if k[0] > k[1]:244                    k = (k[1], k[0])245                pair[k] = pair.get(k, 0) + 1246    return R, weight, pair, by247248249def no_adjacent(n):250    out = []251    for m in range(1, 3 ** (n - 1)):252        x, prev, ok = m, 0, True253        while x:254            d = x % 3255            if d > 1 or (d and prev):256                ok = False257                break258            prev, x = d, x // 3259        if ok:260            out.append(m)261    return out262263264def is_three_power_ratio(a, b):265    if a > b:266        a, b = b, a267    if b % a:268        return False269    q = b // a270    while q % 3 == 0:271        q //= 3272    return q == 1273274275# LAW 1: THE CONSTRAINED TENSOR SQUARE276277def law_tensor():278    bad = 0279    n = 9280    tested = 0281    for s in range(1, 40):282        for t in range(s + 1, 40):283            if gcd(s, t) != 1:284                continue285            tested += 1286            if free_returns(free_automaton(s, t), n) != tensor_returns(s, t, n):287                bad += 1288    check("tensor square reproduces B(s,t) return counts, n <= 9", (tested, bad), (473, 0))289    sizes = []290    for (s, t) in ((365, 1094), (41, 122), (25, 52), (31, 40)):291        sizes.append((len(carry_live(carry(s, t))), len(free_automaton(s, t))))292    check("carry states against reachable B states",293          sizes, [(729, 26931), (81, 835), (38, 393), (35, 354)])294295296# LAW 2: THE BOX CONSTRUCTION297298def law_box():299    n = 9300    bad = 0301    tested = 0302    for s in range(1, 30):303        for t in range(s + 1, 60):304            if gcd(s, t) != 1:305                continue306            tested += 1307            want = free_returns(free_automaton(s, t), n)[n] - 1 - 2 * fibre(s, t, n)308            if want != box(s, t, n):309                bad += 1310    check("box construction against B(s,t), n = 9", (tested, bad), (812, 0))311    got = []312    for (s, t) in ((365, 1094), (41, 122), (122, 123), (1, 2460), (2431, 2458)):313        for n in (9, 12):314            a = tensor_returns(s, t, n)[n] - 1 - 2 * fibre(s, t, n)315            got.append((a, box(s, t, n)))316    check("box against tensor at large multipliers, n = 9 and 12",317          [x for x in got if x[0] != x[1]], [])318    check("witness counts at the four extreme pairs, n = 12",319          [g[0] for g in got], [2, 180, 50, 172, 50, 172, 2, 12, 2, 8])320321322# LAW 3: WEIGHT FOUR AND THE GOLDEN CEILING323324def law_weights():325    R = {}326    weight = {}327    pair = {}328    for n in range(4, 14):329        R[n], weight[n], pair[n], _ = residual(n)330    check("R(n) for n = 4..13", [R[n] for n in range(4, 14)],331          [20, 88, 432, 1624, 5512, 15896, 46064, 124928, 335704, 863848])332    scaled = 0333    bad = 0334    for n in range(5, 14):335        for w, c in weight[n].items():336            if w % 3 == 0:337                scaled += 1338                if weight[n - 1].get(w // 3, 0) != c:339                    bad += 1340    check("R_{3w}(n) = R_w(n-1)", (scaled, bad), (1869, 0))341    check("no witness of weight below four",342          sorted(set(min(weight[n]) for n in range(4, 14))), [4])343    fours = []344    ok = True345    for n in range(4, 13):346        F = no_adjacent(n)347        if len(F) != FIB[n + 1] - 1:348            ok = False349        c = sum(1 for a in F for b in F350                if a != b and gcd(a, b) == 1 and not is_three_power_ratio(a, b))351        fours.append(2 * c)352        if 2 * c != weight[n].get(4, 0):353            ok = False354    check("R_4(n) counted by the no-adjacent-ones set, n = 4..12",355          (ok, fours), (True, [12, 36, 108, 336, 988, 2596, 6672,356                               17480, 45720]))357    check("|F_n| = Fib(n+1) - 1 for n = 4..12",358          [len(no_adjacent(n)) for n in range(4, 13)],359          [FIB[n + 1] - 1 for n in range(4, 13)])360    tops = []361    for n in range(6, 14):362        thr = (3 ** n - 1) // 10363        big = [v for k, v in pair[n].items() if k[1] > thr]364        tops.append((len(big), sorted(set(big))))365    check("every pair above (3^n-1)/10 contributes exactly 4, n = 6..13",366          tops, [(18, [4]), (57, [4]), (163, [4]), (402, [4]), (1019, [4]),367                 (2702, [4]), (7060, [4]), (18607, [4])])368    check("heaviest ray mass is Fib(n+1) - 1, n = 1..40",369          ray_mass(1, 3, 40), [FIB[n + 1] - 1 for n in range(41)])370    N = 40371    ref = [FIB[n + 1] - 1 for n in range(N + 1)]372    tested = 0373    breaches = 0374    ties = []375    for z1 in range(1, 121):376        for z2 in range(z1, 241):377            if gcd(z1, z2) != 1:378                continue379            tested += 1380            m = ray_mass(z1, z2, N)381            for n in range(1, N + 1):382                if m[n] > ref[n]:383                    breaches += 1384            if m[N] == ref[N]:385                ties.append((z1, z2))386    check("golden ceiling M_n(z) <= Fib(n+1) - 1 over the box, n <= 40",387          (tested, breaches, ties), (13158, 0, [(1, 3)]))388    N = 45389    ref = [FIB[m + 1] - 1 for m in range(N + 1)]390    bin3 = [sum(((u >> i) & 1) * 3 ** i for i in range(8)) for u in range(1, 256)]391    nad = no_adjacent(8)392    families = [393        [(a, b) for a in bin3 for b in bin3 if a < b],394        [(a, b) for a in nad for b in nad if a < b],395        [(1, t) for t in range(2, 3000)],396        [(a, a + 1) for a in range(1, 1500)],397        [(a, 3 * a - 1) for a in range(1, 1200)],398        [(a, 3 * a + 1) for a in range(1, 1200)],399    ]400    sizes = []401    breaches = 0402    for family in families:403        seen = 0404        for (a, b) in family:405            if gcd(a, b) != 1:406                continue407            seen += 1408            m = ray_mass(a, b, N)409            for j in range(1, N + 1):410                if m[j] > ref[j]:411                    breaches += 1412        sizes.append(seen)413    check("golden ceiling on six adversarial families, n <= 45",414          (sizes, sum(sizes), breaches),415          ([16940, 253, 2998, 1499, 1199, 1199], 24088, 0))416417418# LAW 4: THE RESIDUAL TO LEVEL 17419420def law_deep():421    import numpy as np422    E, R = [], []423    for n in range(13, 18):424        X = np.zeros(1, dtype=np.int32)425        Y = np.zeros(1, dtype=np.int32)426        for _ in range(n):427            X = np.concatenate((3 * X, 3 * X + 1, 3 * X))428            Y = np.concatenate((3 * Y, 3 * Y, 3 * Y + 1))429        keep = (X > 0) & (Y > 0)430        X, Y = X[keep], Y[keep]431        del keep432        g = np.gcd(X, Y)433        X //= g434        Y //= g435        del g436        key = X.astype(np.int64)437        key *= 3 ** n438        key += Y439        del X, Y440        key.sort()441        idx = np.flatnonzero(np.concatenate(([True], key[1:] != key[:-1])))442        m = np.diff(np.concatenate((idx, [len(key)]))).astype(np.int64)443        e = int((m * m).sum())444        del key, idx, m445        E.append(e)446        R.append(e - (3 ** n - 2 ** (n + 1) + 1) - (3 ** n - 4 * 2 ** n + 2 * n + 3))447    check("E(n) for n = 13..17", E,448          [4003372, 11679626, 34050692, 99800950, 292848756])449    check("R(n) for n = 13..17", R,450          [863848, 2211960, 5549452, 14100688, 35354824])451    peak3 = max(range(len(R)), key=lambda i: R[i] / 3 ** (13 + i))452    fall3 = all(R[i + 1] / 3 ** (14 + i) < R[i] / 3 ** (13 + i) for i in range(4))453    p = (3 + 5 ** 0.5) / 2454    fallp = all(R[i + 1] / p ** (14 + i) < R[i] / p ** (13 + i) for i in range(4))455    check("R/3^n and R/phi^2n both fall through n = 17",456          (peak3, fall3, fallp), (0, True, True))457    check("R(17)/3^17 and R(17)/phi^34 truncated down",458          (int(R[4] / 3 ** 17 * 10 ** 7), int(R[4] / p ** 17 * 10 ** 7)),459          (2737709, 27724831))460461462# LAW 5: THE CEILING AS A THEOREM463464465def direction_automaton(a, b):466    index = {0: 0}467    order = [0]468    edges = []469    i = 0470    while i < len(order):471        c = order[i]472        out = []473        for eps in (0, b, -a):474            if (c + eps) % 3 == 0:475                nxt = (c + eps) // 3476                if nxt not in index:477                    index[nxt] = len(order)478                    order.append(nxt)479                out.append((eps, index[nxt]))480        edges.append(out)481        i += 1482    return order, edges483484485def direction_live(a, b):486    order, edges = direction_automaton(a, b)487    back = [[] for _ in edges]488    for i, out in enumerate(edges):489        for (eps, j) in out:490            back[j].append(i)491    seen = {0}492    stack = [0]493    while stack:494        x = stack.pop()495        for y in back[x]:496            if y not in seen:497                seen.add(y)498                stack.append(y)499    keep = sorted(seen)500    place = {v: i for i, v in enumerate(keep)}501    return ([order[x] for x in keep],502            [[(eps, place[j]) for (eps, j) in edges[x] if j in place]503             for x in keep])504505506def direction_mass(a, b, n):507    order, edges = direction_live(a, b)508    cur = [0] * len(edges)509    cur[0] = 1510    out = [1]511    for _ in range(n):512        nxt = [0] * len(edges)513        for i in range(len(edges)):514            if cur[i]:515                for (eps, j) in edges[i]:516                    nxt[j] += cur[i]517        cur = nxt518        out.append(cur[0])519    return [v - 1 for v in out]520521522def state_profile(a, b, n):523    order, edges = direction_live(a, b)524    tab = [[1 if order[i] == 0 else 0 for i in range(len(edges))]]525    for _ in range(n):526        prev = tab[-1]527        tab.append([sum(prev[j] for (eps, j) in edges[i])528                    for i in range(len(edges))])529    return [max(r) for r in tab]530531532def first_returns(a, b, n):533    order, edges = direction_live(a, b)534    cur = [0] * len(edges)535    cur[0] = 1536    f = [0] * (n + 1)537    for m in range(1, n + 1):538        nxt = [0] * len(edges)539        for i in range(len(edges)):540            if cur[i]:541                for (eps, j) in edges[i]:542                    nxt[j] += cur[i]543        f[m] = nxt[0]544        nxt[0] = 0545        cur = nxt546    return f547548549def valuation3(x):550    k = 0551    while x % 3 == 0:552        x //= 3553        k += 1554    return k555556557def double_branch(order, edges):558    degs = [len(o) for o in edges]559    return any(degs[i] == 2 and all(degs[j] == 2 for (eps, j) in edges[i])560               for i in range(len(edges)))561562563def is_shift_ray(a, b):564    lo, hi = min(a, b), max(a, b)565    if lo != 1:566        return False567    while hi % 3 == 0:568        hi //= 3569    return hi == 1570571572def certificate_holds(a, b, den, alpha, beta):573    order, edges = direction_live(a, b)574    if len(order) != len(alpha) or len(order) != len(beta):575        return False576    if alpha[0] != den or beta[0] != 0:577        return False578    for i in range(len(edges)):579        if alpha[i] < 0:580            return False581        if sum(alpha[j] for (eps, j) in edges[i]) > alpha[i] + beta[i]:582            return False583        if sum(beta[j] for (eps, j) in edges[i]) > alpha[i]:584            return False585    return True586587588def binary_multiples(w, n):589    cur = [0] * w590    cur[0] = 1591    for i in range(n):592        p = pow(3, i, w)593        nxt = list(cur)594        for r in range(w):595            if cur[r]:596                nxt[(r + p) % w] += cur[r]597        cur = nxt598    return cur[0] - 1599600601def shift_product(j, n):602    if n <= j:603        return 1604    total = 1605    for r in range(j):606        m = len(range(r, n - j, j))607        total *= FIB[m + 2]608    return total609610611CERTIFICATES = {612    (1, 90): (18, [18, 0, 5, 6, 9, 2, 10, 4, 8],613              [0, 18, 4, -4, 5, 6, 8, 1, 2]),614    (4, 117): (40, [40, 0, 5, 14, 6, 13, 11, 27, 17, 1],615               [0, 40, 1, -1, 5, 14, 6, 13, 11, 4]),616    (9, 73): (381,617              [381, 0, 46, 49, 78, 17, 101, 23, 66, 163, 39, 83, 232, 32,618               62, 149],619              [0, 381, 32, -32, 46, 49, 62, 16, 17, 101, 23, 66, 149, 14,620               39, 83]),621    (9, 82): (18, [18, 0, 5, 6, 9, 2, 10, 4, 8],622              [0, 18, 4, -4, 5, 6, 8, 1, 2]),623    (9, 235): (2013,624               [2013, 0, 16, 66, 27, 55, 43, 121, 70, 176, 113, 297, 183,625                473, 296, 770, 479, 1243, 775, 11],626               [0, 2013, 11, -11, 16, 66, 27, 55, 43, 121, 70, 176, 113,627                297, 183, 473, 296, 770, 479, 5]),628    (10, 81): (18, [18, 0, 5, 6, 9, 2, 10, 4, 8],629               [0, 18, 4, -4, 5, 6, 8, 1, 2]),630    (13, 108): (40, [40, 0, 5, 14, 6, 13, 11, 27, 17, 1],631                [0, 40, 1, -1, 5, 14, 6, 13, 11, 4]),632    (27, 217): (2013,633                [2013, 0, 16, 66, 27, 55, 43, 121, 70, 176, 113, 297, 183,634                 473, 296, 770, 479, 1243, 775, 11],635                [0, 2013, 11, -11, 16, 66, 27, 55, 43, 121, 70, 176, 113,636                 297, 183, 473, 296, 770, 479, 5]),637    (27, 226): (34, [34, 0, 0, 1, 1, 0, 1, 1, 2, 1, 3, 5, 8, 13, 20, 1],638                [0, 34, 1, -1, 0, 1, 1, 0, 1, 1, 2, 3, 5, 8, 14, -1]),639}640641642def law_ceiling():643    tested = mismatch = 0644    for a in range(1, 40):645        for b in range(1, 40):646            if gcd(a, b) != 1:647                continue648            tested += 1649            if direction_mass(a, b, 20) != ray_mass(a, b, 20):650                mismatch += 1651    check("direction carry automaton against the gasket-digit automaton",652          (tested, mismatch), (947, 0))653    box = occupied = bounded = single = k1 = nodouble = 0654    twoplus = anomaly = states = renewal = 0655    exceptions = []656    twos = []657    for z1 in range(1, 121):658        for z2 in range(z1, 241):659            if gcd(z1, z2) != 1:660                continue661            box += 1662            d = (z1 % 3 == 0) + (z2 % 3 == 0) + ((z1 + z2) % 3 == 0)663            if d > 1:664                twoplus += 1665            order, edges = direction_live(z1, z2)666            if len(order) == 1:667                continue668            if d == 0:669                anomaly += 1670            occupied += 1671            states = max(states, len(order))672            if all(-z1 <= 2 * c <= z2 for c in order):673                bounded += 1674            degs = [len(o) for o in edges]675            classes = set(order[i] % 3 for i in range(len(order))676                          if degs[i] == 2)677            if max(degs) <= 2 and len(classes) <= 1:678                single += 1679            q = z1 if z1 % 3 == 0 else (z2 if z2 % 3 == 0 else z1 + z2)680            if valuation3(q) == 1:681                k1 += 1682            if not double_branch(order, edges):683                nodouble += 1684            else:685                exceptions.append((z1, z2))686            f = first_returns(z1, z2, 46)687            if f[1] == 1 and all(688                    sum(f[j] * FIB[m + 1 - j] for j in range(2, m + 1))689                    <= FIB[m - 1] for m in range(2, 47)):690                renewal += 1691            if f[2]:692                twos.append((z1, z2))693    check("carry states of a direction lie in [-a/2, b/2]",694          (box, occupied, states, bounded), (13158, 218, 37, 218))695    check("out-degree at most two, branch states in one class mod 3",696          single, 218)697    check("at most one of z1, z2, w is divisible by three, and no occupied "698          "direction has none of them divisible", (twoplus, anomaly), (0, 0))699    check("directions settled by the branch argument over the box",700          (k1, nodouble, exceptions),701          (107, 206, [(1, 9), (1, 27), (1, 81), (1, 90), (4, 117), (9, 73),702                      (9, 82), (9, 235), (10, 81), (13, 108), (27, 217),703                      (27, 226)]))704    check("renewal criterion on every occupied direction of the box, n <= 46",705          (renewal, twos), (218, [(1, 3)]))706    good = sorted(k for k, v in CERTIFICATES.items()707                  if certificate_holds(k[0], k[1], v[0], v[1], v[2]))708    check("Fibonacci certificates for the nine non-shift exceptions",709          good, [(1, 90), (4, 117), (9, 73), (9, 82), (9, 235), (10, 81),710                 (13, 108), (27, 217), (27, 226)])711    identity = all(FIB[p + 2] * FIB[q + 2]712                   == FIB[p + q + 3] - FIB[p + 1] * FIB[q + 1]713                   for p in range(60) for q in range(60))714    cases = over = strict = 0715    for j in range(1, 14):716        for m in range(1, 46):717            cases += 1718            if shift_product(j, m) > FIB[m + 1]:719                over += 1720            if j >= 2 and m >= 2 and shift_product(j, m) >= FIB[m + 1]:721                strict += 1722    check("Fibonacci product identity and the shift-ray ceiling",723          (identity, cases, over, strict), (True, 585, 0, 0))724    tested = breaches = 0725    for z1 in range(1, 31):726        for z2 in range(z1, 61):727            if gcd(z1, z2) != 1:728                continue729            tested += 1730            if direction_mass(z1, z2, 12)[12] > binary_multiples(z1 + z2, 12):731                breaches += 1732    check("ray mass is at most the count of binary multiples of the weight",733          (tested, breaches), (829, 0))734    check("binary multiples of the weight outgrow the ceiling",735          ([binary_multiples(w, 24) for w in (4, 10, 28, 82)], FIB[25] - 1),736          ([4196351, 1683971, 613817, 228519], 75024))737    profile = state_profile(1, 9, 12)738    fails = 0739    for z1 in range(1, 121):740        for z2 in range(z1, 241):741            if gcd(z1, z2) != 1:742                continue743            h = state_profile(z1, z2, 20)744            if any(h[m] > h[m - 1] + h[m - 2] for m in range(2, 21)):745                fails += 1746    check("the state maximum breaks the Fibonacci recursion",747          (profile[:7], profile[4] > profile[3] + profile[2], fails),748          ([1, 1, 1, 2, 4, 6, 9], True, 8))749750751# LAW 6: THE PROVED CASES ON THE ADVERSARIAL FAMILIES752753754def law_families():755    bin3 = [sum(((u >> i) & 1) * 3 ** i for i in range(8)) for u in range(1, 256)]756    nad = no_adjacent(8)757    families = [758        [(a, b) for a in bin3 for b in bin3 if a < b],759        [(a, b) for a in nad for b in nad if a < b],760        [(1, t) for t in range(2, 3000)],761        [(a, a + 1) for a in range(1, 1500)],762        [(a, 3 * a - 1) for a in range(1, 1200)],763        [(a, 3 * a + 1) for a in range(1, 1200)],764    ]765    seen = set()766    multiplicity = 0767    for family in families:768        for (a, b) in family:769            if gcd(a, b) != 1:770                continue771            multiplicity += 1772            seen.add((a, b))773    inside = set((a, b) for (a, b) in seen if a <= 120 and b <= 240)774    zero = nodouble = shift = 0775    left = []776    for (a, b) in sorted(seen - inside):777        order, edges = direction_live(a, b)778        if len(order) == 1:779            zero += 1780        elif not double_branch(order, edges):781            nodouble += 1782        elif is_shift_ray(a, b):783            shift += 1784        else:785            left.append((a, b))786    check("the six families overlap, and their union splits",787          (multiplicity, len(seen), len(inside), len(seen) - len(inside),788           zero, nodouble, shift, len(left)),789          (24088, 23435, 717, 22718, 20945, 1693, 3, 77))790    breaches = 0791    wn, wd = 0, 1792    for (a, b) in left:793        m = direction_mass(a, b, 60)794        for j in range(1, 61):795            top = FIB[j + 1] - 1796            if m[j] > top:797                breaches += 1798            if top and m[j] * wd > wn * top:799                wn, wd = m[j], top800    check("the directions left to the enumeration hold to n = 60",801          (breaches, -(-wn * 10000 // wd)), (0, 1516))802803804# LAW 7: THE GOLDEN POTENTIAL805806807def qnorm(p, q, d):808    if d < 0:809        p, q, d = -p, -q, -d810    g = gcd(gcd(abs(p), abs(q)), d)811    if g > 1:812        p //= g813        q //= g814        d //= g815    return (p, q, d)816817818def qadd(x, y):819    return qnorm(x[0] * y[2] + y[0] * x[2],820                 x[1] * y[2] + y[1] * x[2], x[2] * y[2])821822823def qsub(x, y):824    return qnorm(x[0] * y[2] - y[0] * x[2],825                 x[1] * y[2] - y[1] * x[2], x[2] * y[2])826827828def qmul(x, y):829    return qnorm(x[0] * y[0] + x[1] * y[1],830                 x[0] * y[1] + x[1] * y[0] + x[1] * y[1], x[2] * y[2])831832833def qinv(x):834    p, q, d = x835    n = p * p + p * q - q * q836    return qnorm(d * (p + q), -d * q, n)837838839def qsgn(x):840    p, q, d = x841    hi, lo = 2 * p + q, q842    if hi >= 0 and lo >= 0:843        return 0 if hi == 0 and lo == 0 else 1844    if hi <= 0 and lo <= 0:845        return 0 if hi == 0 and lo == 0 else -1846    s = hi * hi - 5 * lo * lo847    if hi > 0:848        return 1 if s > 0 else (0 if s == 0 else -1)849    return -1 if s > 0 else (0 if s == 0 else 1)850851852QZERO = (0, 0, 1)853QONE = (1, 0, 1)854QPHI = (0, 1, 1)855QINVPHI = (-1, 1, 1)856QINVPHI2 = (2, -1, 1)857858859def golden_potential(a, b):860    order, edges = direction_live(a, b)861    n = len(order)862    if n == 1:863        return None864    m = n - 1865    rows = [[QZERO] * (m + 1) for _ in range(m)]866    for r in range(m):867        rows[r][r] = QPHI868        for (eps, j) in edges[r + 1]:869            if j == 0:870                rows[r][m] = qadd(rows[r][m], QONE)871            else:872                rows[r][j - 1] = qsub(rows[r][j - 1], QONE)873    for col in range(m):874        piv = None875        for r in range(col, m):876            if qsgn(rows[r][col]):877                piv = r878                break879        if piv is None:880            return "singular"881        rows[col], rows[piv] = rows[piv], rows[col]882        scale = qinv(rows[col][col])883        rows[col] = [qmul(v, scale) if qsgn(v) else QZERO for v in rows[col]]884        for r in range(m):885            if r != col and qsgn(rows[r][col]):886                f = rows[r][col]887                rows[r] = [qsub(rows[r][k], qmul(f, rows[col][k]))888                           if qsgn(rows[col][k]) else rows[r][k]889                           for k in range(m + 1)]890    u = [QONE] + [rows[r][m] for r in range(m)]891    return u, order, edges892893894def potential_value(u, edges):895    tot = QZERO896    for (eps, j) in edges[0]:897        if j:898            tot = qadd(tot, u[j])899    return tot900901902def potential_valid(u, edges):903    if u[0] != QONE or any(qsgn(v) <= 0 for v in u):904        return False905    for i in range(1, len(edges)):906        s = QZERO907        for (eps, j) in edges[i]:908            s = qadd(s, u[j])909        if qsgn(qsub(qmul(QPHI, u[i]), s)) < 0:910            return False911    return True912913914def stress_directions():915    out = set()916    for j in range(1, 9):917        out.add((1, 3 ** j))918    for i in range(0, 7):919        for j in range(0, 8):920            a, b = 3 ** i, 3 ** i + 3 ** j921            g = gcd(a, b)922            a, b = a // g, b // g923            if b > a:924                out.add((a, b))925    for k in range(1, 7):926        for m in range(1, 40):927            for (a, b) in ((m, 3 ** k * m + 1), (1, 3 ** k * m),928                           (3 ** k, 3 ** k + m), (3 ** k, 3 ** k + 3 * m + 1)):929                if a > b:930                    a, b = b, a931                if a >= 1 and b > a and b <= 6000 and gcd(a, b) == 1:932                    out.add((a, b))933    return out934935936def law_partition():937    box = []938    for z1 in range(1, 121):939        for z2 in range(z1, 241):940            if gcd(z1, z2) == 1:941                box.append((z1, z2))942    occupied = certified = passing = 0943    over = []944    values = {}945    for (a, b) in box:946        r = golden_potential(a, b)947        if r is None:948            continue949        occupied += 1950        u, order, edges = r951        if potential_valid(u, edges):952            certified += 1953        tot = potential_value(u, edges)954        values.setdefault(tot, []).append((a, b))955        if qsgn(qsub(QINVPHI2, tot)) >= 0:956            passing += 1957        else:958            over.append((a, b))959    check("the golden potential certifies the box away from the shift rays",960          (occupied, certified, passing, over),961          (218, 218, 214, [(1, 3), (1, 9), (1, 27), (1, 81)]))962    ranked = []963    for _ in range(3):964        top = None965        for k in values:966            if k not in ranked and (top is None or qsgn(qsub(k, top)) > 0):967                top = k968        ranked.append(top)969    check("the top of the golden potential over the box",970          ([(k, sorted(values[k])) for k in ranked], len(values)),971          ([((-1, 1, 1), [(1, 3), (1, 9), (1, 27), (1, 81)]),972            ((2, -1, 1), [(1, 12), (3, 10), (4, 9)]),973            ((-14, 9, 2), [(1, 90), (9, 82), (10, 81)])], 57))974    tail = peak = 0975    for (a, b) in box:976        r = golden_potential(a, b)977        if r is None:978            continue979        u, order, edges = r980        tot = potential_value(u, edges)981        f = first_returns(a, b, 46)982        s = QZERO983        p = QONE984        for j in range(1, 47):985            p = qmul(p, QINVPHI)986            if j >= 2 and f[j]:987                s = qadd(s, qmul((f[j], 0, 1), p))988        if f[1] == 1 and qsgn(qsub(qmul(QINVPHI, tot), s)) >= 0:989            tail += 1990        if qsgn(qsub(QINVPHI, tot)) >= 0:991            peak += 1992    check("the potential dominates the first-return series and no direction "993          "beats the shift rays", (tail, peak), (218, 218))994    seen = set(box)995    bin3 = [sum(((u >> i) & 1) * 3 ** i for i in range(8)) for u in range(1, 256)]996    nad = no_adjacent(8)997    families = [998        [(a, b) for a in bin3 for b in bin3 if a < b],999        [(a, b) for a in nad for b in nad if a < b],1000        [(1, t) for t in range(2, 3000)],1001        [(a, a + 1) for a in range(1, 1500)],1002        [(a, 3 * a - 1) for a in range(1, 1200)],1003        [(a, 3 * a + 1) for a in range(1, 1200)],1004    ]1005    for family in families:1006        for (a, b) in family:1007            if gcd(a, b) == 1:1008                seen.add((a, b))1009    seen |= stress_directions()1010    total = len(seen)1011    occupied = certified = passing = states = 01012    over = []1013    gap = []1014    best = QZERO1015    attain = []1016    for (a, b) in sorted(seen):1017        r = golden_potential(a, b)1018        if r is None:1019            continue1020        occupied += 11021        u, order, edges = r1022        states = max(states, len(order))1023        if potential_valid(u, edges):1024            certified += 11025        tot = potential_value(u, edges)1026        if qsgn(qsub(QINVPHI2, tot)) >= 0:1027            passing += 11028            if qsgn(qsub(tot, best)) > 0:1029                best, attain = tot, [(a, b)]1030            elif tot == best:1031                attain.append((a, b))1032        else:1033            over.append((a, b, tot))1034            if tot != QINVPHI:1035                gap.append((a, b, tot))1036    check("the golden potential over the box, the six families and the "1037          "high-valuation stress list",1038          (total, occupied, certified, states, passing),1039          (36037, 1995, 1995, 256, 1987))1040    check("every direction over the criterion is a shift ray at exactly "1041          "one over phi", (over, gap),1042          ([(1, 3 ** j, QINVPHI) for j in range(1, 9)], []))1043    check("the criterion is attained exactly on the supergolden directions",1044          (best, attain), (QINVPHI2, [(1, 12), (3, 10), (4, 9)]))1045    zero = residue = 01046    for (z1, z2) in box:1047        if z1 % 3 and z2 % 3:1048            residue += 11049        if len(direction_live(z1, z2)[0]) == 1:1050            zero += 11051    check("occupancy needs three to divide one coordinate",1052          (len(box), residue, zero, len(box) - zero),1053          (13158, 6566, 12940, 218))105410551056# LAW 8: THE DEGREE POTENTIAL105710581059def split3(a, b):1060    q, p = (a, b) if a % 3 == 0 else (b, a)1061    k = 01062    q1 = q1063    while q1 % 3 == 0:1064        q1 //= 31065        k += 11066    return p, k, q1106710681069def phipow(e):1070    r = QONE1071    for _ in range(abs(e)):1072        r = qmul(r, QPHI if e > 0 else QINVPHI)1073    return r107410751076def degree_potential(edges):1077    return [QONE if len(o) == 2 else QINVPHI for o in edges]107810791080def super_solution(edges, pi):1081    for i in range(1, len(edges)):1082        s = QZERO1083        for (eps, j) in edges[i]:1084            s = qadd(s, pi[j])1085        if qsgn(qsub(qmul(QPHI, pi[i]), s)) < 0:1086            return False1087    return True108810891090def sweep(edges, d):1091    pi = degree_potential(edges)1092    pi[0] = QONE1093    for _ in range(d):1094        nxt = [QONE] + [QZERO] * (len(edges) - 1)1095        for i in range(1, len(edges)):1096            s = QZERO1097            for (eps, j) in edges[i]:1098                s = qadd(s, pi[j])1099            nxt[i] = qmul(QINVPHI, s)1100        pi = nxt1101    tot = QZERO1102    for (eps, j) in edges[0]:1103        if j:1104            tot = qadd(tot, pi[j])1105    return tot110611071108def base3(x):1109    d = []1110    while x:1111        d.append(x % 3)1112        x //= 31113    return d or [0]111411151116def burst_floor(a, b, k, q1):1117    sign = 1 if b % 3 == 0 else -11118    out = []1119    for m in range(1, 3 ** k):1120        if m % 3 != 1:1121            continue1122        x, ok = m, True1123        while x:1124            if x % 3 > 1:1125                ok = False1126                break1127            x //= 31128        if ok:1129            out.append(sign * q1 * m)1130    return out113111321133def law_degree():1134    seen = set()1135    for z1 in range(1, 121):1136        for z2 in range(z1, 241):1137            if gcd(z1, z2) == 1:1138                seen.add((z1, z2))1139    bin3 = [sum(((u >> i) & 1) * 3 ** i for i in range(8)) for u in range(1, 256)]1140    nad = no_adjacent(8)1141    families = [1142        [(a, b) for a in bin3 for b in bin3 if a < b],1143        [(a, b) for a in nad for b in nad if a < b],1144        [(1, t) for t in range(2, 3000)],1145        [(a, a + 1) for a in range(1, 1500)],1146        [(a, 3 * a - 1) for a in range(1, 1200)],1147        [(a, 3 * a + 1) for a in range(1, 1200)],1148    ]1149    for family in families:1150        for (a, b) in family:1151            if gcd(a, b) == 1:1152                seen.add((a, b))1153    seen |= stress_directions()1154    triple = eligible = occupied = 01155    resbad = burst = burstbad = 01156    valid = nodouble = doublevalid = settled = 01157    depths = {}1158    missed = []1159    k1 = k1class = 01160    k1bad = []1161    quantbad = []1162    attain = []1163    for (a, b) in sorted(seen):1164        p, k, q1 = split3(a, b)1165        if a % 3 == 0 or b % 3 == 0:1166            triple += 11167            if (q1 - p) % 3 == 0:1168                eligible += 11169        order, edges = direction_live(a, b)1170        if len(order) == 1:1171            continue1172        occupied += 11173        if (q1 - p) % 3:1174            resbad += 11175        f = first_returns(a, b, max(k, 2))1176        if any(f[j] for j in range(2, k + 1)):1177            burst += 11178        if not double_branch(order, edges):1179            nodouble += 11180        pi = degree_potential(edges)1181        if super_solution(edges, pi):1182            valid += 11183            if double_branch(order, edges):1184                doublevalid += 11185            hit = None1186            for d in range(25):1187                if qsgn(qsub(QINVPHI2, sweep(edges, d))) >= 0:1188                    hit = d1189                    break1190            if hit is None:1191                missed.append((a, b))1192            else:1193                settled += 11194                depths[hit] = depths.get(hit, 0) + 11195        else:1196            missed.append((a, b))1197        u = golden_potential(a, b)[0]1198        tot = potential_value(u, edges)1199        place = {c: i for i, c in enumerate(order)}1200        floor = burst_floor(a, b, k, q1)1201        rung = QZERO1202        for c in floor:1203            if c in place:1204                rung = qadd(rung, u[place[c]])1205        if len(floor) != 2 ** (k - 1) or \1206                qmul(phipow(-(k - 1)), rung) != tot:1207            burstbad += 11208        if 2 * p <= 3 ** k * q1:1209            c = p if b % 3 == 0 else -p1210            if c not in place or u[place[c]] != QINVPHI:1211                burstbad += 11212        if k != 1:1213            continue1214        k1 += 11215        t = valuation3(q1 - p) if q1 != p else None1216        if t is None:1217            continue1218        bound = qmul(QINVPHI, qsub(QONE, phipow(-max(t, 2))))1219        if qsgn(qsub(bound, tot)) < 0:1220            quantbad.append((a, b, t))1221        if tot == bound:1222            attain.append((a, b))1223        if t <= 2:1224            k1class += 11225            if qsgn(qsub(QINVPHI2, tot)) < 0:1226                k1bad.append((a, b, t))1227    check("occupancy needs the residue match q1 = p mod 3",1228          (triple, eligible, occupied, resbad), (20193, 15556, 1995, 0))1229    check("no first return has length between two and v3(q)",1230          (occupied, burst), (1995, 0))1231    check("the burst identity and the value at the near predecessor",1232          (occupied, burstbad), (1995, 0))1233    check("the degree potential is a super-solution beyond the branch case",1234          (occupied, nodouble, valid, doublevalid), (1995, 1902, 1968, 66))1235    check("the swept degree potential settles all but the shift rays and "1236          "twenty-one directions",1237          (settled, sorted(depths.items()),1238           [z for z in missed if is_shift_ray(*z)],1239           [z for z in missed if not is_shift_ray(*z)]),1240          (1966, [(1, 1804), (3, 101), (4, 44), (5, 11), (6, 6)],1241           [(1, 3 ** j) for j in range(1, 9)],1242           [(1, 756), (1, 2196), (1, 2214), (1, 2268), (1, 2430), (9, 2188),1243            (10, 2187), (13, 1080), (13, 3267), (27, 730), (27, 2188),1244            (28, 729), (28, 2187), (40, 1053), (81, 2188), (82, 2187),1245            (91, 2214), (121, 3159), (243, 2188), (244, 2187), (819, 2539)]))1246    check("the golden partition bound at v3(q) = 1",1247          (k1, k1class, k1bad, quantbad, attain),1248          (757, 512, [], [], [(1, 12), (3, 10)]))1249    short = []1250    for a in range(1, 260):1251        for b in range(1, 1100):1252            if gcd(a, b) != 1 or (a % 3 and b % 3):1253                continue1254            if len(direction_live(a, b)[0]) == 1:1255                continue1256            f = first_returns(a, b, 4)1257            if f[2] or f[3]:1258                short.append((min(a, b), max(a, b), f[2], f[3]))1259    check("the first returns of length two and three are classified",1260          sorted(set(short)),1261          [(1, 3, 1, 0), (1, 9, 0, 1), (1, 12, 0, 1), (3, 10, 0, 1),1262           (4, 9, 0, 1)])1263    tested = bad = 01264    quartic = qmul(qmul(QPHI, QPHI), qmul(QPHI, QPHI))1265    for a in range(1, 40):1266        for b in range(1, 120):1267            if gcd(a, b) != 1 or (a % 3 and b % 3):1268                continue1269            order, edges = direction_live(a, b)1270            if len(order) == 1:1271                continue1272            tested += 11273            tot = potential_value(golden_potential(a, b)[0], edges)1274            den = qsub(QINVPHI2, qmul(QINVPHI, tot))1275            closed = qinv(den) if qsgn(den) > 0 else None1276            mass = direction_mass(a, b, 46)1277            partial = QZERO1278            w = QONE1279            for n in range(47):1280                partial = qadd(partial, qmul((mass[n] + 1, 0, 1), w))1281                w = qmul(w, QINVPHI)1282            layer = QZERO1283            for m in range(1, 3 ** 11 // (a + b)):1284                x, y = base3(a * m), base3(b * m)1285                if max(x) < 2 and max(y) < 2 and \1286                        not any(dx and dy for dx, dy in zip(x, y)):1287                    layer = qadd(layer, phipow(-len(base3((a + b) * m))))1288            ok = closed is None or qsgn(qsub(closed, partial)) >= 01289            ok = ok and (closed is not None1290                         and qsgn(qsub(quartic, closed)) >= 0) == \1291                (qsgn(qsub(QINVPHI2, tot)) >= 0)1292            ok = ok and (qsgn(qsub(QPHI, layer)) >= 01293                         or qsgn(qsub(tot, QINVPHI2)) > 0)1294            bad += not ok1295    check("the criterion restated as a golden series and as a multiplier count",1296          (tested, bad), (111, 0))129712981299def main():1300    law_tensor()1301    law_box()1302    law_weights()1303    law_deep()1304    law_ceiling()1305    law_families()1306    law_partition()1307    law_degree()1308    print("gasket witness weights: every law holds")130913101311main()