jsr.py

1from fractions import Fraction23CODES = list(range(1, 16))4FRAME = ("gamma", "h", "v", "phi")56# GEOMETRY78def tile(c):9    return [(c >> 0) & 3, (c >> 2) & 3]1011def draw(word):12    rows, n = [1], 113    for c in word:14        t = tile(c)15        out = []16        for r in rows:17            for p in (0, 1):18                acc = 019                for j in range(n):20                    if (r >> j) & 1:21                        acc |= t[p] << (2 * j)22                out.append(acc)23        rows, n = out, 2 * n24    return rows, n2526def fill(rows):27    return sum(bin(r).count("1") for r in rows)2829def hruns(rows):30    return sum(bin(r & ~(r << 1)).count("1") for r in rows)3132def vruns(rows):33    total, prev = 0, 034    for r in rows:35        total += bin(r & ~prev).count("1")36        prev = r37    return total3839def comp(rows, n):40    seen = [0] * len(rows)41    total = 042    for i in range(len(rows)):43        free = rows[i] & ~seen[i]44        while free:45            b = free & -free46            total += 147            stack = [(i, b.bit_length() - 1)]48            seen[i] |= b49            while stack:50                a, j = stack.pop()51                for da, dj in ((1, 0), (-1, 0), (0, 1), (0, -1)):52                    x, y = a + da, j + dj53                    if 0 <= x < len(rows) and 0 <= y < n:54                        m = 1 << y55                        if (rows[x] & m) and not (seen[x] & m):56                            seen[x] |= m57                            stack.append((x, y))58            free = rows[i] & ~seen[i]59    return total6061CACHE = {}6263def obs(word):64    if word not in CACHE:65        rows, n = draw(word)66        CACHE[word] = (comp(rows, n), hruns(rows), vruns(rows), fill(rows))67    return CACHE[word]6869# EXACT LINEAR ALGEBRA7071def solve(basis, target):72    m, ncol = len(basis), len(target)73    aug = [[basis[i][j] for i in range(m)] + [target[j]] for j in range(ncol)]74    piv, top = [], 075    for col in range(m):76        p = next((k for k in range(top, ncol) if aug[k][col] != 0), None)77        if p is None:78            continue79        aug[top], aug[p] = aug[p], aug[top]80        d = aug[top][col]81        aug[top] = [x / d for x in aug[top]]82        for k in range(ncol):83            if k != top and aug[k][col] != 0:84                f = aug[k][col]85                aug[k] = [a - f * b for a, b in zip(aug[k], aug[top])]86        piv.append((top, col))87        top += 188    if any(aug[k][m] != 0 for k in range(top, ncol)):89        return None90    x = [Fraction(0)] * m91    for r, c in piv:92        x[c] = aug[r][m]93    return x9495def mul(A, B):96    d = len(A)97    return [[sum(A[i][k] * B[k][j] for k in range(d)) for j in range(d)] for i in range(d)]9899def det(A):100    d = len(A)101    if d == 1:102        return A[0][0]103    return sum((-1) ** j * A[0][j] * det([row[:j] + row[j + 1:] for row in A[1:]]) for j in range(d))104105def eye(d):106    return [[1 if i == j else 0 for j in range(d)] for i in range(d)]107108def apply(A, x):109    return [sum(A[i][j] * x[j] for j in range(len(x))) for i in range(len(A))]110111def floor_root(num, den, L, digits=9):112    scale = 10 ** digits113    lo = int((num / den) ** (1.0 / L) * scale) - 4114    while (lo + 1) ** L * den <= num * scale ** L:115        lo += 1116    while lo ** L * den > num * scale ** L:117        lo -= 1118    return "%d.%0*d" % (lo // scale, digits, lo % scale)119120def ceil_root(num, den, L, digits=9):121    scale = 10 ** digits122    hi = int((num / den) ** (1.0 / L) * scale) + 4123    while (hi - 1) ** L * den >= num * scale ** L:124        hi -= 1125    while hi ** L * den < num * scale ** L:126        hi += 1127    return "%d.%0*d" % (hi // scale, digits, hi % scale)128129def words(alphabet, length):130    if length == 0:131        yield ()132        return133    for w in words(alphabet, length - 1):134        for c in alphabet:135            yield w + (c,)136137def seeded(alphabet, count, lo, hi):138    state, out = 20260906, []139    for _ in range(count):140        state = (1103515245 * state + 12345) % (1 << 31)141        length = lo + state % (hi - lo + 1)142        w = []143        for _ in range(length):144            state = (1103515245 * state + 12345) % (1 << 31)145            w.append(alphabet[state % len(alphabet)])146        out.append(tuple(w))147    return out148149# THE REPRESENTATION150151SUFFIXES = [()] + [(c,) for c in CODES] + [(a, b) for a in CODES for b in CODES]152153def row_of(u):154    return [Fraction(obs(u + s)[0]) for s in SUFFIXES]155156def build():157    basis, rows, mats = [()], [row_of(())], {c: [] for c in CODES}158    i = 0159    while i < len(basis):160        for c in CODES:161            r = row_of(basis[i] + (c,))162            x = solve(rows, r)163            if x is None:164                basis.append(basis[i] + (c,))165                rows.append(r)166                x = solve(rows, r)167            mats[c].append(x)168        i += 1169    d = len(basis)170    for c in CODES:171        for line in mats[c]:172            line.extend([Fraction(0)] * (d - len(line)))173        mats[c] = [[int(x) for x in line] for line in mats[c]]174    return basis, mats175176print("JSR SCHEDULES")177print()178print("== THE REPRESENTATION ==")179basis, M = build()180print("Hankel rank %d, basis words %s" % (len(basis), [w if w else "e" for w in basis]))181lam = [1] + [0] * (len(basis) - 1)182gamma = [obs(u)[0] for u in basis]183h = [obs(u)[1] for u in basis]184v = [obs(u)[2] for u in basis]185phi = [obs(u)[3] for u in basis]186print("lambda %s, gamma %s, h %s, v %s, phi %s" % (lam, gamma, h, v, phi))187classes = {}188for c in CODES:189    classes.setdefault(tuple(map(tuple, M[c])), []).append(c)190CLASS = sorted(classes.values(), key=lambda g: g[0])191print("class partition %s" % (CLASS,))192assert len(CLASS) == 6193noncomm = sum(1 for i in range(6) for j in range(i + 1, 6)194              if mul(M[CLASS[i][0]], M[CLASS[j][0]]) != mul(M[CLASS[j][0]], M[CLASS[i][0]]))195print("non-commuting class pairs %d of 15" % noncomm)196for g in CLASS:197    print("  M_%s = %s" % (g, M[g[0]]))198199# THE OBSERVABLE FRAME200201print()202print("== THE OBSERVABLE FRAME ==")203G = [[gamma[i], h[i], v[i], phi[i]] for i in range(4)]204print("frame columns gamma %s, h %s, v %s, phi %s, det %d" % (gamma, h, v, phi, det(G)))205assert h == apply(M[3], gamma) and v == apply(M[5], gamma) and phi == apply(M[1], gamma)206print("h, v and phi are M_3 gamma, M_5 gamma and M_1 gamma, so all four frame vectors are comp read at one appended letter")207assert abs(det(G)) == 1208Ginv = None209adj = [[Fraction((-1) ** (i + j) * det([r[:i] + r[i + 1:] for k, r in enumerate(G) if k != j]), det(G))210        for j in range(4)] for i in range(4)]211Ginv = adj212assert mul(Ginv, G) == eye(4)213T = {}214for c in CODES:215    t = mul(Ginv, mul(M[c], G))216    T[c] = [[int(x) for x in row] for row in t]217    assert all(Fraction(T[c][i][j]) == t[i][j] for i in range(4) for j in range(4))218print("frame images, columns of T_c in the frame (gamma, h, v, phi):")219for g in CLASS:220    c = g[0]221    rows, n = draw((c,))222    ob = obs((c,))223    cols = [[T[c][i][j] for i in range(4)] for j in range(4)]224    csum = [sum(col) for col in cols]225    print("  code %-2d k=%d  M_c gamma=%s  M_c h=%s  M_c v=%s  M_c phi=%s  column sums %s"226          % (c, ob[3], cols[0], cols[1], cols[2], cols[3], csum))227    assert csum == [ob[0], ob[1], ob[2], ob[3]]228    assert max(csum) == ob[3]229print("every column sum is (comp, rows, cols, fill) of the letter's own tile, and the largest is the fill")230231# THE FOUR TRANSFER LAWS, CHECKED AGAINST DRAWN CELLS232233print()234print("== THE TRANSFER LAWS ==")235short = [w for L in range(0, 4) for w in words(CODES, L)]236seeds = seeded(CODES, 240, 4, 7)237sweep = short + seeds238bad = 0239for w in sweep:240    y = obs(w)241    for c in CODES:242        z = obs(w + (c,))243        pred = tuple(sum(y[i] * T[c][i][j] for i in range(4)) for j in range(4))244        if pred != z:245            bad += 1246assert bad == 0247print("(comp, H, V, fill) at wc equals (comp, H, V, fill) at w times T_c on %d words, all %d of length at most 3 and %d seeded of length 4 to 7, times 15 letters, %d mismatches"248      % (len(sweep), len(short), len(seeds), bad))249rep_bad = 0250for w in sweep:251    P = eye(4)252    for c in w:253        P = mul(P, M[c])254    val = sum(lam[i] * sum(P[i][j] * gamma[j] for j in range(4)) for i in range(4))255    hv = sum(lam[i] * sum(P[i][j] * h[j] for j in range(4)) for i in range(4))256    vv = sum(lam[i] * sum(P[i][j] * v[j] for j in range(4)) for i in range(4))257    fv = sum(lam[i] * sum(P[i][j] * phi[j] for j in range(4)) for i in range(4))258    if (val, hv, vv, fv) != obs(w):259        rep_bad += 1260assert rep_bad == 0261print("lambda M_w applied to gamma, h, v, phi reads comp, H, V, fill on the same %d words, %d of length at most 3 and %d seeded, %d mismatches"262      % (len(sweep), len(short), len(seeds), rep_bad))263264# TRIANGULARITY AND SPECTRUM265266print()267print("== TRIANGULARITY ==")268for g in CLASS:269    c = g[0]270    assert all(T[c][i][j] == 0 for i in range(4) for j in range(4) if j > i)271    assert all(T[c][i][j] >= 0 for i in range(4) for j in range(4))272    diag = [T[c][i][i] for i in range(4)]273    poly = [t for t in range(5)]274    vals = [det([[t * (1 if i == j else 0) - M[c][i][j] for j in range(4)] for i in range(4)]) for t in poly]275    pred = [(t - diag[0]) * (t - diag[1]) * (t - diag[2]) * (t - diag[3]) for t in poly]276    assert vals == pred277    print("  code %-2d T_c lower triangular, nonnegative, diagonal %s, char poly roots %s, rho = %d = fill"278          % (c, diag, sorted(diag), max(diag)))279    assert max(diag) == obs((c,))[3]280281# THE CROSS-POLYTOPE CERTIFICATE282283print()284print("== THE CERTIFICATE ==")285print("P = conv{+/- gamma, +/- h, +/- v, +/- phi}; residual is k_c minus the frame l1 norm of the image")286for g in CLASS:287    c = g[0]288    k = obs((c,))[3]289    res = []290    for j, name in enumerate(FRAME):291        col = [T[c][i][j] for i in range(4)]292        res.append(k - sum(abs(x) for x in col))293    print("  code %-2d k=%d residuals at gamma, h, v, phi: %s" % (c, k, res))294    assert all(r >= 0 for r in res)295print("M_c P is inside k_c P at every code, so the gauge of P, the frame l1 norm N(a gamma + b h + c v + d phi) = |a| + |b| + |c| + |d|, is an extremal norm with ||M_c|| = k_c")296297# JSR AND LSR OVER EVERY SUBFAMILY298299print()300print("== JSR AND LSR ==")301K = {c: obs((c,))[3] for c in CODES}302pairs = [(a, b) for i, a in enumerate(CODES) for b in CODES[i + 1:]]303print("pairs of distinct letters %d, pairs of distinct classes %d" % (len(pairs), 15))304print("class pair table (a, b, JSR, LSR):")305for i in range(6):306    for j in range(i + 1, 6):307        a, b = CLASS[i][0], CLASS[j][0]308        print("  (%d, %d) JSR %d LSR %d" % (a, b, max(K[a], K[b]), min(K[a], K[b])))309distinct = sum(1 for a, b in pairs if K[a] != K[b])310print("of the %d pairs of distinct letters, %d carry two different fills and %d carry one" % (len(pairs), distinct, len(pairs) - distinct))311print("JSR = max fill and LSR = min fill on every one of them, each attained by a one-letter word")312313# THE BRACKET ON TWO NAMED PAIRS314315print()316print("== THE BRACKET ==")317LMAX = 16318for pair in [(3, 6), (3, 7)]:319    a, b = pair320    r = max(K[a], K[b])321    cur = {(): eye(4)}322    best_lo, best_word, best_up, best_min = {}, {}, {}, {}323    for L in range(1, LMAX + 1):324        nxt = {}325        for w, P in cur.items():326            for c in pair:327                nxt[w + (c,)] = mul(P, T[c])328        cur = nxt329        lo, arg, up, worst = -1, None, 0, None330        for w, P in cur.items():331            sr = max(P[i][i] for i in range(4))332            if sr > lo:333                lo, arg = sr, w334            if worst is None or sr < worst:335                worst = sr336            nrm = max(sum(P[i][j] for i in range(4)) for j in range(4))337            if nrm > up:338                up = nrm339        best_lo[L], best_word[L], best_up[L], best_min[L] = lo, arg, up, worst340        assert lo <= up341        assert lo == r ** L342        assert up == r ** L343        assert worst == min(K[a], K[b]) ** L344    print("pair {%d, %d}, JSR candidate %d" % (a, b, r))345    for L in (1, 2, 4, 8, 12, LMAX):346        w = "".join(str(x) for x in best_word[L])347        print("  L=%-2d lower %s from word %s, spectral radius %d; upper %s from max frame 1-norm %d"348              % (L, floor_root(best_lo[L], 1, L), w, best_lo[L], ceil_root(best_up[L], 1, L), best_up[L]))349    print("  the best word's rate is exactly %d at every length 1 to %d, so it stops improving at length 1" % (r, LMAX))350    print("  the worst word's spectral radius is exactly %d^L at every length, so the lower spectral radius is %d"351          % (min(K[a], K[b]), min(K[a], K[b])))352    for k in (2, 4, 6):353        num = K[a] ** k + K[b] ** k354        lift_lo = floor_root(num, 2, k)355        lift_up = ceil_root(num, 1, k)356        print("  Blondel-Nesterov lifting at k=%d on the 2 letters: lower 2^(-1/%d) (%d^%d + %d^%d)^(1/%d) = %s, upper (%d^%d + %d^%d)^(1/%d) = %s"357              % (k, k, K[a], k, K[b], k, k, lift_lo, K[a], k, K[b], k, k, lift_up))358        assert Fraction(lift_lo) <= Fraction(lift_up)359        assert Fraction(lift_lo) <= r <= Fraction(lift_up)360        assert num >= r ** k361    print("  the scan alone brackets it at [%s, %s] at L=%d, the lifting alone at [%s, %s] at k=6, and the certificate closes it at %d exactly"362          % (floor_root(best_lo[LMAX], 1, LMAX), ceil_root(best_up[LMAX], 1, LMAX), LMAX,363             floor_root(K[a] ** 6 + K[b] ** 6, 2, 6), ceil_root(K[a] ** 6 + K[b] ** 6, 1, 6), r))364365# THE TELESCOPE BLOCK366367print()368print("== THE TELESCOPE ==")369A = [[M[3][i][j] for j in range(2)] for i in range(2)]370B = [[M[6][i][j] for j in range(2)] for i in range(2)]371assert all(M[3][i][j] == 0 for i in range(4) for j in (2, 3))372assert all(M[6][i][j] == 0 for i in range(4) for j in (2, 3))373print("on {3, 6} the last two columns of both matrices vanish, leading blocks A = %s, B = %s" % (A, B))374p = [1, 2]375assert apply(A, p) == [2 * x for x in p]376assert apply(B, [1, 1]) == [2 * x for x in p] and apply(B, p) == [2 * x for x in p]377print("A p = 2 p and B has range span(p) at p = %s, so span(p) is invariant under both and the pair is reducible" % p)378print("in the basis (gamma, p) both blocks are lower triangular, A = diag(1, 2) and B = [[0, 0], [2, 2]]")379for L in (1, 2, 3, 6, 10):380    vals = set()381    for w in words((3, 6), L):382        P = eye(4)383        for c in w:384            P = mul(P, T[c])385        vals.add(max(P[i][i] for i in range(4)))386    print("  every word of length %d has spectral radius %s, so its rate is exactly 2" % (L, sorted(vals)))387    assert vals == {2 ** L}388389# THE DEAD ROUTES390391print()392print("== WHAT DIES ==")393for g in CLASS:394    c = g[0]395    assert apply(M[c], phi) == [K[c] * x for x in phi]396print("phi = %s is a right eigenvector of every M_c with eigenvalue the fill, so span(phi) is a common invariant line and no subfamily is irreducible" % phi)397orbit = {tuple(phi)}398for _ in range(4):399    orbit |= {tuple(Fraction(t, K[c]) for t in apply(M[c], list(x))) for x in orbit for c in (3, 7)}400    orbit = {tuple(int(t) if t.denominator == 1 else t for t in x) for x in orbit}401assert orbit == {tuple(phi)}402print("the polytope algorithm started at the leading eigenvector of the length-1 spectrum maximizing product, each image normalised by the letter fill k_c, closes on %d vertex up to sign and spans dimension 1 of 4" % len(orbit))403norms = []404for L in (1, 2, 4, 8, 16, 32):405    P = eye(4)406    for _ in range(L):407        P = mul(P, M[3])408    norms.append(max(abs(P[i][j]) for i in range(4) for j in range(4)))409print("max |entry(M_3^L)| at L = 1, 2, 4, 8, 16, 32 reads %s, matching 2^(L+2) - 2" % norms)410assert norms == [2 ** (L + 2) - 2 for L in (1, 2, 4, 8, 16, 32)]411print("comp(A_(3^L)) = %s at the same lengths, so the JSR rate log 2 is not the component rate 0"412      % [obs((3,) * L)[0] for L in (1, 2, 4, 8)])413print()414print("every assertion passed")