verify.py

1from fractions import Fraction2from math import comb, gcd, log, sqrt34# PINCER56F = ((0, 0), (1, 0), (0, 1))7FS = set(F)8PHI = (1 + sqrt(5)) / 29KAPPA = 3 - log(5, 3)1011def near(got, want, tol, what):12    assert abs(got - want) <= tol, "%s: got %.12f want %.12f" % (what, got, want)1314def same(got, want, what):15    assert got == want, "%s: got %r want %r" % (what, got, want)1617def klass(a, b):18    if a % 3 and b % 3:19        if (a - b) % 3 == 0:20            return "eq", 021        k = 022        s = a + b23        while s % 3 == 0:24            s //= 325            k += 126        return "opp", k27    x = a if a % 3 == 0 else b28    k = 029    while x % 3 == 0:30        x //= 331        k += 132    return "div", k3334def oriented(a, b):35    if a % 3 == 0:36        return a, b37    if b % 3 == 0:38        return b, a39    return a, b4041def edges(a, b):42    out = [[] for _ in range(a * b)]43    for c1 in range(a):44        for c2 in range(b):45            s = c1 * b + c246            for d in range(3):47                x = a * d + c148                y = b * d + c249                if (x % 3, y % 3) in FS:50                    out[s].append((d, (x // 3) * b + (y // 3)))51    return out5253def rowsums(out, w):54    v = [1] * len(out)55    for _ in range(w):56        v = [sum(v[t] for _, t in row) for row in out]57    return v5859def fib(n):60    x, y = 0, 161    for _ in range(n):62        x, y = y, x + y63    return x6465def dfree(k, w):66    p = 167    for r in range(k):68        m = len([i for i in range(1, w + 1) if i % k == r])69        p *= fib(m + 2)70    return p7172def dfree_brute(k, w):73    c = 074    for mask in range(1 << w):75        ok = True76        for i in range(w - k):77            if (mask >> i) & 1 and (mask >> (i + k)) & 1:78                ok = False79                break80        if ok:81            c += 182    return c8384def rays(h):85    return [(a, b) for a in range(1, h + 1) for b in range(a, h + 1) if gcd(a, b) == 1]8687def typ(cls, a, b, s):88    c1, c2 = divmod(s, b)89    return c1 % 3 if cls == "div" else (c1 + c2) % 39091# LADDER9293def rung(lam, order):94    k = order - log(lam, 3)95    return k / (2 * k + 2 - KAPPA)9697def delta_counts(K):98    types = [(b, c, comb(K, b) * comb(K - b, c)) for b in range(K + 1) for c in range(K + 1 - b)]99    out = {}100    for b, c, m in types:101        for u, v, w in types:102            out[(b - u, c - v)] = out.get((b - u, c - v), 0) + m * w103    return out104105def carry(K):106    r = (K - 1) // 2107    states = [(i, j) for i in range(-r, r + 1) for j in range(-r, r + 1)]108    at = {s: i for i, s in enumerate(states)}109    M = [[0] * len(states) for _ in states]110    for s in states:111        for (d1, d2), w in delta_counts(K).items():112            z1, z2 = s[0] + d1, s[1] + d2113            if z1 % 3 == 0 and z2 % 3 == 0:114                M[at[s]][at[(z1 // 3, z2 // 3)]] += w115    return M, at[(0, 0)], r116117def energies(K, levels):118    M, zero, _ = carry(K)119    v = [0] * len(M)120    v[zero] = 1121    out = [1]122    for _ in range(levels):123        v = [sum(v[i] * M[i][j] for i in range(len(v))) for j in range(len(v))]124        out.append(v[zero])125    return out126127def energy_direct(K, a):128    pts = [(0, 0)]129    for l in range(a):130        pts = [(x + dx * 3 ** l, y + dy * 3 ** l) for x, y in pts for dx, dy in F]131    dist = {(0, 0): 1}132    for _ in range(K):133        fresh = {}134        for (x, y), c in dist.items():135            for dx, dy in pts:136                key = (x + dx, y + dy)137                fresh[key] = fresh.get(key, 0) + c138        dist = fresh139    return sum(c * c for c in dist.values())140141def charpoly(M):142    n = len(M)143    A = [[Fraction(x) for x in row] for row in M]144    coeffs = [Fraction(1)]145    N = [[Fraction(int(i == j)) for j in range(n)] for i in range(n)]146    for k in range(1, n + 1):147        AN = [[sum(A[i][t] * N[t][j] for t in range(n)) for j in range(n)] for i in range(n)]148        c = -sum(AN[i][i] for i in range(n)) / k149        coeffs.append(c)150        N = [[AN[i][j] + (c if i == j else 0) for j in range(n)] for i in range(n)]151    return [int(c) for c in coeffs]152153def polymul(a, b):154    out = [0] * (len(a) + len(b) - 1)155    for i, x in enumerate(a):156        for j, y in enumerate(b):157            out[i + j] += x * y158    return out159160def trim(p):161    i = 0162    while i < len(p) - 1 and p[i] == 0:163        i += 1164    return p[i:]165166def deriv(p):167    n = len(p) - 1168    return trim([p[i] * (n - i) for i in range(n)]) if n else [Fraction(0)]169170def polyrem(a, b):171    a = [Fraction(x) for x in a]172    b = [Fraction(x) for x in b]173    while len(a) >= len(b) and any(a):174        f = a[0] / b[0]175        for i in range(len(b)):176            a[i] -= f * b[i]177        a = trim(a)178    return a179180def sturm(p):181    chain = [trim([Fraction(x) for x in p])]182    chain.append(deriv(chain[0]))183    while len(chain[-1]) > 1:184        r = polyrem(chain[-2], chain[-1])185        if not any(r):186            break187        chain.append([-c for c in r])188    return chain189190def variations(vals):191    s = [v for v in vals if v != 0]192    return sum(1 for i in range(len(s) - 1) if (s[i] > 0) != (s[i + 1] > 0))193194def value(p, t):195    v = Fraction(0)196    for c in p:197        v = v * t + c198    return v199200def sign_changes(chain, t):201    return variations([value(q, t) for q in chain])202203def sign_changes_infinity(chain):204    return variations([q[0] for q in chain])205206def squarefree(p):207    chain = sturm(p)208    return len(chain[-1]) == 1 and chain[-1][0] != 0209210def real_roots_above(p, t):211    chain = sturm(p)212    return sign_changes(chain, t) - sign_changes_infinity(chain)213214def real_roots_between(p, lo, hi):215    chain = sturm(p)216    return sign_changes(chain, lo) - sign_changes(chain, hi)217218QUARTIC = [1, -7833, 7916949, -850684437, 13054946580]219220OTHERS = [[1, 0], [1, -120], [1, -450, 12231], [1, -2190, 282096, -5186835], [1, -990, 116154, -2569725]]221222MULTIPLICITY = [6, 1, 1, 2, 2]223224ENERGIES = [1, 4653, 28967859, 190911254427, 1270015973323281, 8461182216374750493]225226LO = Fraction(66641136625, 10 ** 7)227228HI = Fraction(66641136626, 10 ** 7)229230def ladder():231    for K in range(1, 9):232        r = (K - 1) // 2233        assert (r + K) // 3 <= r, "order %d carry box: got (r + K) // 3 = %d want <= %d" % (2 * K, (r + K) // 3, r)234        assert max(max(abs(d1), abs(d2)) for d1, d2 in delta_counts(K)) == K, "order %d digit spread: want %d" % (2 * K, K)235    same(energies(2, 6), [15 ** a for a in range(7)], "E_4(G_a) = 15^a")236    M, zero, r = carry(5)237    same((len(M), r), (25, 2), "order-10 carry matrix shape")238    same(energies(5, 5), ENERGIES, "E_10(G_a) from the carry matrix")239    for a in range(1, 4):240        same(energy_direct(5, a), ENERGIES[a], "E_10(G_%d) by direct convolution" % a)241    poly = charpoly(M)242    same(len(poly) - 1, 25, "charpoly degree")243    product = QUARTIC244    for f, m in zip(OTHERS, MULTIPLICITY):245        for _ in range(m):246            product = polymul(product, f)247    same(product, poly, "charpoly factorisation")248    assert squarefree(QUARTIC), "quartic factor: got a repeated root want squarefree"249    same(real_roots_above(QUARTIC, HI), 0, "quartic real roots above 6664.1136626")250    same(real_roots_between(QUARTIC, LO, HI), 1, "quartic real roots in the bracket")251    for f in OTHERS[1:]:252        assert squarefree(f), "factor %r: got a repeated root want squarefree" % f253        same(real_roots_above(f, HI), 0, "factor %r real roots above 6664.1136626" % f)254    k10 = 10 - log(float(HI), 3)255    b10 = rung(float(HI), 10)256    assert k10 > 1.985805792698, "certified kappa_10: got %.12f want > 1.985805792698" % k10257    assert b10 > 0.447597813453, "certified rung 10: got %.12f want > 0.447597813453" % b10258    assert b10 > 0.4475978, "short edge: got %.7f want > 0.4475978" % b10259    assert b10 < rung(6664.113662506, 10), "certified rung 10 must sit below the floating value"260    print("ladder: order-10 carry matrix 25 states, energies to a = 5, exact charpoly factorisation, Sturm gives lambda_10 < 6664.1136626 and beta_0^(10) > 0.447597813453")261262# PINCER263264def constants():265    L = log(PHI, 3)266    near(1 / (2 - L), 0.6402121938, 5e-10, "top edge")267    near((1 - L) / (2 - L), 0.3597878, 5e-7, "c star")268    x = 1.5269    for _ in range(80):270        x -= (x ** 3 - x * x - 1) / (3 * x * x - 2 * x)271    near(x, 1.4655712319, 5e-10, "supergolden root")272    near(1 / (2 - log(x, 3)), 0.6053028664, 5e-10, "supergolden edge")273    near(2 / (3 + log(5, 3)), 0.447930988, 5e-9, "ladder cap")274    near(rung(456 + 3 * sqrt(11017), 8), 0.446717310462, 5e-12, "rung 8")275    near(rung(6664.113662506, 10), 0.447597813454, 5e-12, "rung 10")276    near(1 / (2 - log(1.0639086, 3)), 0.5145062, 5e-8, "quarantined 0.5145062")277    near(1 / (2 - log(1.0997454, 3)), 0.5226147, 5e-8, "quarantined 0.5226147")278    print("constants: edge 0.6402121938, rungs 0.446717310462 / 0.447597813454, cap 0.447930988")279280def subsets():281    for k in range(1, 6):282        for w in range(1, 15):283            same(dfree(k, w), dfree_brute(k, w), "D_%d(%d)" % (k, w))284    for k in range(1, 7):285        for w in range(1, 61):286            assert dfree(k, w) <= PHI ** (w + k), "D_%d(%d) exceeds phi^(w+k)" % (k, w)287    print("subsets: D_k(w) exact for w <= 14, k <= 5; D_k(w) <= phi^(w+k) for k <= 6, w <= 60")288289def branching(h):290    n2 = 0291    for (a, b) in rays(h):292        cls, k = klass(a, b)293        a, b = oriented(a, b)294        out = edges(a, b)295        for s, row in enumerate(out):296            if cls == "eq":297                assert len(row) <= 1, "eq ray (%d,%d) state %d: got %d digits want <= 1" % (a, b, s, len(row))298                continue299            t = typ(cls, a, b, s)300            want = {0: 2, 1: 1, 2: 0}[t] if cls == "div" else {0: 1, 1: 2, 2: 0}[t]301            same(len(row), want, "ray (%d,%d) state %d type %d" % (a, b, s, t))302            if want == 2:303                n2 += 1304                same(len({d % 3 for d, _ in row}), 2, "ray (%d,%d) state %d digits" % (a, b, s))305    print("branching: exact 2/1/0 counts at every state of every ray of height <= %d (%d branching states)" % (h, n2))306307def delay(h):308    tested = 0309    for (a, b) in rays(h):310        cls, k = klass(a, b)311        if cls == "eq":312            continue313        a, b = oriented(a, b)314        out = edges(a, b)315        for s in range(len(out)):316            paths = [(s, [])]317            for _ in range(k):318                nxt = []319                for st, seq in paths:320                    for d, t in out[st]:321                        nxt.append((t, seq + [(d, typ(cls, a, b, t))]))322                paths = nxt323            if not paths:324                continue325            for j in range(k - 1):326                vals = {seq[j][1] for _, seq in paths}327                same(len(vals), 1, "ray (%d,%d) state %d step %d predetermination" % (a, b, s, j + 1))328            heads = {}329            for _, seq in paths:330                heads.setdefault(seq[0][0], set()).add(seq[k - 1][1])331            if len(heads) == 2:332                (d1, t1), (d2, t2) = list(heads.items())333                same(len(t1) * len(t2), 1, "ray (%d,%d) state %d split determinism" % (a, b, s))334                assert t1 != t2, "ray (%d,%d) state %d: got equal types %r want distinct" % (a, b, s, t1)335                tested += 1336    print("delay: predetermination and the k-step split at every state of every div/opp ray of height <= %d (%d splits)" % (h, tested))337338def burst(h, wmax):339    for (a, b) in rays(h):340        cls, k = klass(a, b)341        a, b = oriented(a, b)342        out = edges(a, b)343        v = [1] * len(out)344        for w in range(1, wmax + 1):345            v = [sum(v[t] for _, t in row) for row in out]346            cap = 1 if cls == "eq" else dfree(k, w)347            got = max(v)348            assert got <= cap, "ray (%d,%d) w=%d: got P_w %d want <= %d" % (a, b, w, got, cap)349    print("burst: P_w(s) <= D_k(w) at every start state, every ray of height <= %d, every w <= %d" % (h, wmax))350351def saturation():352    same(rowsums(edges(9, 1), 36)[0], 45765225, "P_36 at (9,1)")353    same(dfree(2, 36), 45765225, "D_2(36)")354    same(rowsums(edges(27, 1), 36)[0], 53582633, "P_36 at (27,1)")355    same(dfree(3, 36), 53582633, "D_3(36)")356    print("saturation: P_36 = D_2(36) = 45765225 at (9,1) and P_36 = D_3(36) = 53582633 at (27,1)")357358def catalogue():359    cat = rays(40)360    same(len(cat), 490, "primitive rays of height <= 40")361    shifts = [(1, 3), (1, 9), (1, 27)]362    body = [r for r in cat if r != (1, 1)]363    same(len(body), 489, "catalogue after removing (1,1)")364    nonshift = [r for r in body if r not in shifts]365    same(len(nonshift), 486, "non-shift rays")366    best = []367    wsum = 0.0368    wnum = 0.0369    for (a, b) in body:370        w = 1.0 / (a * a + b * b)371        if (a, b) in shifts:372            r = PHI373        else:374            r = max(rowsums(edges(*oriented(a, b)), 55)) ** (1.0 / 55)375            assert r < PHI, "ray (%d,%d): got certificate %.10f want < phi" % (a, b, r)376            best.append((r, (a, b)))377        wsum += w378        wnum += w * r379    best.sort(reverse=True)380    top = [(p, round(r, 10)) for r, p in best[:3]]381    same(top, [((4, 9), 1.481203426), ((3, 10), 1.4765525267), ((1, 12), 1.4728541511)], "top three certificates")382    near(wnum / wsum, 1.0997454, 5e-8, "inverse-square-weighted upper mean")383    print("catalogue: 490 rays of height <= 40, 489 = 486 + 3 after removing (1,1), top certificate 1.4812034260 at (4,9), weighted upper mean 1.0997454")384385def census(jmax):386    for j in range(1, jmax + 1):387        lo, hi = 3 ** j, 3 ** (j + 1)388        counts = {}389        for (a, b) in rays(hi - 1):390            if max(a, b) < lo:391                continue392            cls, k = klass(a, b)393            counts[(cls, k)] = counts.get((cls, k), 0) + 1394        for (cls, k), c in counts.items():395            if cls == "eq":396                cap = 9 * 3 ** (2 * j)397            elif cls == "div":398                cap = 9 * 3 ** (2 * j - k)399                assert k <= j, "octave %d div: got k %d want <= %d" % (j, k, j)400            else:401                cap = 18 * 3 ** (2 * j - k)402                assert k <= j + 1, "octave %d opp: got k %d want <= %d" % (j, k, j + 1)403            assert c <= cap, "octave %d %s k=%d: got %d want <= %d" % (j, cls, k, c, cap)404        d1 = counts.get(("div", 1), 0)405        want = 3 ** (2 * j - 1)406        assert d1 >= want, "octave %d depth-one census: got %d want >= %d" % (j, d1, want)407    print("census: div, opp and eq upper bounds and the depth-one count against 3^(2j-1) for every octave 1 <= j <= %d" % jmax)408409def main():410    constants()411    ladder()412    subsets()413    branching(24)414    delay(24)415    burst(24, 20)416    saturation()417    catalogue()418    census(4)419    print("all green")420421if __name__ == "__main__":422    main()