memory.py

1import csv2import math3import tempfile4import subprocess5import sys6import time7from itertools import permutations8from pathlib import Path910HERE = Path(__file__).resolve().parent11ROOT = HERE.parents[3]1213CASES = [(1, 1), (1, 2), (1, 3), (1, 4), (2, 1), (2, 2)]14A000616 = {1: 3, 2: 6, 3: 22, 4: 402}15NAMED = {16    (1, -1, -1): "golden",17    (1, -1, 0, -1): "supergolden",18    (1, 0, -1, -1): "plastic",19    (1, -1, -1, -1): "tribonacci",20}2122# WINDOWS2324def digits_of(w, d, k):25    a = 1 << d26    return [(w >> (d * (k - 1 - j))) & (a - 1) for j in range(k)]272829def word_of(ds, d):30    w = 031    for c in ds:32        w = (w << d) | c33    return w3435# GROUP3637def bd_apply(c, perm, flip, d):38    out = 039    for i in range(d):40        b = ((c >> i) & 1) ^ ((flip >> i) & 1)41        out |= b << perm[i]42    return out434445def window_tables(d, k, reversal):46    nw = 1 << (d * k)47    revs = (False, True) if reversal else (False,)48    tables = set()49    for perm in permutations(range(d)):50        for flip in range(1 << d):51            for rev in revs:52                t = []53                for w in range(nw):54                    ds = [bd_apply(c, perm, flip, d) for c in digits_of(w, d, k)]55                    if rev:56                        ds = ds[::-1]57                    t.append(word_of(ds, d))58                tables.add(tuple(t))59    return sorted(tables)606162def code_array(table, nbits):63    n = 1 << nbits64    arr = [0] * n65    for code in range(1, n):66        low = code & -code67        arr[code] = arr[code ^ low] | (1 << table[low.bit_length() - 1])68    return arr697071def orbit_walk(maps, ncodes):72    seen = bytearray(ncodes)73    reps, sizes = [], []74    for c in range(ncodes):75        if seen[c]:76            continue77        orb = {m[c] for m in maps}78        for x in orb:79            seen[x] = 180        reps.append(c)81        sizes.append(len(orb))82    return reps, sizes8384# TRANSFER MATRIX8586def transfer(code, d, k):87    a = 1 << d88    if k == 1:89        return [[bin(code).count("1")]]90    s_n = a ** (k - 1)91    m = [[0] * s_n for _ in range(s_n)]92    for s in range(s_n):93        for c in range(a):94            w = s * a + c95            if (code >> w) & 1:96                m[s][w % s_n] = 197    return m9899100def matmul(x, y):101    n = len(x)102    return [[sum(x[i][t] * y[t][j] for t in range(n)) for j in range(n)] for i in range(n)]103104105def charpoly(m):106    n = len(m)107    cs = [1]108    mj = [row[:] for row in m]109    for j in range(1, n + 1):110        tr = sum(mj[i][i] for i in range(n))111        cj = -tr // j112        cs.append(cj)113        if j < n:114            t = [row[:] for row in mj]115            for i in range(n):116                t[i][i] += cj117            mj = matmul(m, t)118    return tuple(cs)119120121def poly_str(cs):122    n = len(cs) - 1123    parts = []124    for i, c in enumerate(cs):125        if c == 0:126            continue127        e = n - i128        mono = "1" if e == 0 else ("x" if e == 1 else f"x^{e}")129        if abs(c) != 1 or e == 0:130            mono = f"{abs(c)}*{mono}" if e else f"{abs(c)}"131        parts.append(("- " if c < 0 else "+ ") + mono)132    s = " ".join(parts)133    return s[2:] if s.startswith("+ ") else "-" + s[2:]134135136def divides(p, q):137    p = list(p)138    q = list(q)139    while len(q) >= len(p):140        if q[0] % p[0]:141            return False142        f = q[0] // p[0]143        for i, c in enumerate(p):144            q[i] -= f * c145        if q[0] != 0:146            return False147        q.pop(0)148    return all(c == 0 for c in q)149150# PARI151152def pari_roots(polys, scratch):153    body = (154        "for(i = 1, #V, p = V[i]; f = factor(p); best = -1; bg = 0; "155        "for(j = 1, #f~, g = f[j, 1]; if(poldegree(g) > 0, r = polrootsreal(g); "156        "if(#r > 0, m = vecmax(r); if(m > best, best = m; bg = g)))); "157        "s = 0; r = polroots(bg); for(t = 1, #r, if(abs(r[t] - best) > 1e-25, s = max(s, abs(r[t])))); "158        'print(i, "|", Vec(bg), "|", best, "|", s))'159    )160    lines = [f"V = [{', '.join(poly_str(p) for p in polys)}];", body, "quit()"]161    src = scratch / "perron.gp"162    src.write_text("\n".join(lines) + "\n")163    out = subprocess.run(164        ["gp", "-q", str(src)],165        capture_output=True,166        text=True,167        stdin=subprocess.DEVNULL,168    )169    if out.returncode != 0:170        sys.exit(out.stderr)171    got = {}172    for line in out.stdout.strip().splitlines():173        idx, vec, val, sec = line.split("|")174        cs = tuple(int(t) for t in vec.strip().strip("[]").split(","))175        got[int(idx) - 1] = (cs, val.strip(), float(sec))176    return [got[i] for i in range(len(polys))]177178# CENSUS179180def burnside(d, k, reversal):181    total = 0182    tables = window_tables(d, k, reversal)183    for t in tables:184        seen = [False] * len(t)185        cycles = 0186        for w in range(len(t)):187            if seen[w]:188                continue189            cycles += 1190            x = w191            while not seen[x]:192                seen[x] = True193                x = t[x]194        total += 1 << cycles195    return total // len(tables)196197198def product_codes(d, k):199    a = 1 << d200    out = set()201    for f in range(1 << a):202        full = 0203        for w in range(1 << (d * k)):204            if all((f >> c) & 1 for c in digits_of(w, d, k)):205                full |= 1 << w206        out.add(full)207    return out208209210def main():211    scratch = Path(sys.argv[1]) if len(sys.argv) > 1 else Path(tempfile.gettempdir())212    t_all = time.time()213    census, table = [], []214    for d, k in CASES:215        t0 = time.time()216        nbits = 1 << (d * k)217        ncodes = 1 << nbits218        maps_g = [code_array(t, nbits) for t in window_tables(d, k, True)]219        maps_b = [code_array(t, nbits) for t in window_tables(d, k, False)]220        reps, sizes = orbit_walk(maps_g, ncodes)221        reps_b, _ = orbit_walk(maps_b, ncodes)222        orbit_of = dict(zip(reps, sizes))223        products = product_codes(d, k)224        t_orbit = time.time() - t0225226        t0 = time.time()227        polys = {}228        for c in reps:229            polys.setdefault(charpoly(transfer(c, d, k)), []).append(c)230        t_poly = time.time() - t0231232        t0 = time.time()233        keys = sorted(polys)234        roots = pari_roots(keys, scratch)235        t_pari = time.time() - t0236237        rows = []238        strict = {}239        for key, (minp, val, sec) in zip(keys, roots):240            strict[minp] = sec241            for c in polys[key]:242                w = bin(c).count("1")243                rho = float(val)244                zero = w > 0 and divides(minp, [1] + [0] * (k - 1) + [-w])245                rows.append((c, w, minp, val, rho, zero))246        n_zero = sum(1 for r in rows if r[5])247        bad = [r for r in rows if r[5] and r[0] not in products]248        live = [r for r in rows if r[4] > 0]249        kappas = [math.log2(r[1]) / k - math.log2(r[4]) for r in live]250        k_max = max(kappas)251        tied = [r for kap, r in zip(kappas, live) if kap > k_max - 1e-12]252        rho1 = [r for r in live if abs(r[4] - 1.0) < 1e-12]253        by_poly = {}254        for c, w, minp, val, rho, zero in rows:255            e = by_poly.setdefault(minp, {"rho": val, "classes": 0, "rules": 0, "code": c, "w": w, "ks": []})256            e["classes"] += 1257            e["rules"] += orbit_of[c]258            if rho > 0 and w > 0:259                e["ks"].append(math.log2(w) / k - math.log2(rho))260            if c < e["code"]:261                e["code"], e["w"] = c, w262        weak = {263            m: e for m, e in by_poly.items() if float(e["rho"]) > 0 and strict[m] >= float(e["rho"]) - 1e-20264        }265        n_weak = len(weak)266        roots_seen = {m for m in by_poly if float(by_poly[m]["rho"]) > 0 and strict[m] < float(by_poly[m]["rho"]) - 1e-20}267        n_radical = 0268        for m in weak:269            step = 0270            for i, c in enumerate(m):271                if c:272                    step = math.gcd(step, len(m) - 1 - i)273            if step > 1 and tuple(m[i] for i in range(0, len(m), step)) in roots_seen:274                n_radical += 1275        entry = {276                "D": d,277                "k": k,278                "rules": ncodes,279                "classes_G": len(reps),280                "classes_B": len(reps_b),281                "group_G": len(maps_g),282                "group_B": len(maps_b),283                "charpolys": len(keys),284                "perron_polys": len(by_poly),285                "classes_kappa0": n_zero,286                "kappa0_nonproduct": len(bad),287                "kappa_min": round(min(kappas), 6),288                "kappa_max": round(k_max, 6),289                "kappa_max_code": min(r[0] for r in tied),290                "kappa_max_ties": len(tied),291                "kappa_max_windows": max(r[1] for r in tied),292                "rho1_windows": max((r[1] for r in rho1), default=0),293                "weak_perron_polys": n_weak,294                "weak_are_radicals": n_radical,295                "burnside_G": burnside(d, k, True),296                "burnside_B": burnside(d, k, False),297                "dead_classes": len(rows) - len(live),298                "a000616": A000616.get(d * k, ""),299                "t_orbit": round(t_orbit, 2),300                "t_poly": round(t_poly, 2),301                "t_pari": round(t_pari, 2),302        }303        census.append(entry)304        for minp, e in sorted(by_poly.items(), key=lambda kv: -float(kv[1]["rho"])):305            table.append(306                {307                    "D": d,308                    "k": k,309                    "minpoly": poly_str(minp),310                    "degree": len(minp) - 1,311                    "rho": e["rho"][:14],312                    "classes": e["classes"],313                    "rules": e["rules"],314                    "code": e["code"],315                    "windows": e["w"],316                    "kappa_min": round(min(e["ks"]), 6) if e["ks"] else "",317                    "kappa_max": round(max(e["ks"]), 6) if e["ks"] else "",318                    "strict": int(strict[minp] < float(e["rho"]) - 1e-20) if float(e["rho"]) > 0 else "",319                    "name": NAMED.get(minp, ""),320                }321            )322        print(323            f"D={d} k={k} rules={ncodes} classes(G)={len(reps)} classes(B)={len(reps_b)} "324            f"charpoly={len(keys)} perron={len(by_poly)} weak={n_weak} kappa0={n_zero} nonproduct={len(bad)} "325            f"kappa_max={round(k_max, 6)} ties={len(tied)} "326            f"orbit={t_orbit:.2f}s poly={t_poly:.2f}s pari={t_pari:.2f}s"327        )328        for r in bad:329            print(f"  NONPRODUCT kappa=0 code={r[0]} windows={r[1]} minpoly={poly_str(r[2])}")330331    with (HERE / "census.csv").open("w", newline="") as fh:332        w = csv.DictWriter(fh, fieldnames=list(census[0]))333        w.writeheader()334        w.writerows(census)335    with (HERE / "classes.csv").open("w", newline="") as fh:336        w = csv.DictWriter(fh, fieldnames=list(table[0]))337        w.writeheader()338        w.writerows(table)339    print(f"live classes {sum(r['classes_G'] - r['dead_classes'] for r in census)}")340    seen_polys = {}341    for r in table:342        seen_polys.setdefault((r["D"], r["k"]), set()).add(r["minpoly"])343    for a, b in (((1, 3), (1, 4)), ((1, 2), (1, 3)), ((1, 1), (1, 2)), ((2, 1), (2, 2))):344        miss = seen_polys[a] - seen_polys[b]345        print(f"nest {a} -> {b}: {len(seen_polys[a])} of {len(seen_polys[b])}, missing {len(miss)}")346    t0 = time.time()347    print("burnside D=1 G  " + ", ".join(str(burnside(1, k, True)) for k in range(1, 9)))348    print("burnside D=1 B  " + ", ".join(str(burnside(1, k, False)) for k in range(1, 9)))349    print("burnside D=2 G  " + ", ".join(str(burnside(2, k, True)) for k in range(1, 5)))350    print("burnside D=2 B  " + ", ".join(str(burnside(2, k, False)) for k in range(1, 5)))351    print(f"burnside extension {time.time() - t0:.2f}s")352    print(f"total {time.time() - t_all:.2f}s")353354# CLOSED FORM355356def classes_closed(k):357    if k % 2 == 0:358        m = k // 2359        return (360            (1 << (2 ** (2 * m) - 2))361            + (1 << (2 ** (2 * m - 1) - 2))362            + (1 << (2 ** (2 * m - 1) + 2 ** (m - 1) - 1))363        )364    m = (k - 1) // 2365    return (366        (1 << (2 ** (2 * m + 1) - 2))367        + (1 << (2 ** (2 * m) - 1))368        + (1 << (2 ** (2 * m) + 2 ** m - 2))369    )370371372def word_orbits(d, k):373    tables = window_tables(d, k, True)374    seen = [False] * (1 << (d * k))375    n = 0376    for w in range(len(seen)):377        if seen[w]:378            continue379        n += 1380        for t in tables:381            seen[t[w]] = True382    return n383384385def verb_burnside(top=11):386    t0 = time.time()387    cycle = [burnside(1, k, True) for k in range(1, top + 1)]388    closed = [classes_closed(k) for k in range(1, top + 1)]389    words = [word_orbits(1, k) for k in range(1, top + 1)]390    print("classes G_(1,k) k=1..8  " + ", ".join(str(v) for v in cycle[:8]))391    for k in range(9, top + 1):392        print(f"classes G_(1,{k})  {cycle[k - 1]}")393    agree = sum(1 for a, b in zip(cycle, closed) if a == b)394    print(f"closed form against the cycle index k=1..{top}: {agree} of {top} agree")395    for k in range(1, top + 1):396        if cycle[k - 1] != closed[k - 1]:397            print(f"  MISMATCH k={k}")398    print("digits k=9,10,11  " + ", ".join(str(len(str(cycle[k - 1]))) for k in (9, 10, 11)))399    print("words under the same group k=1..8  " + ", ".join(str(v) for v in words[:8]))400    print(f"burnside {time.time() - t0:.2f}s")401402# DETERMINANT403404def det_bareiss(m):405    n = len(m)406    a = [row[:] for row in m]407    sign = 1408    prev = 1409    for i in range(n - 1):410        if a[i][i] == 0:411            p = -1412            for r in range(i + 1, n):413                if a[r][i]:414                    p = r415                    break416            if p < 0:417                return 0418            a[i], a[p] = a[p], a[i]419            sign = -sign420        for r in range(i + 1, n):421            for c in range(i + 1, n):422                a[r][c] = (a[r][c] * a[i][i] - a[r][i] * a[i][c]) // prev423            a[r][i] = 0424        prev = a[i][i]425    return sign * a[n - 1][n - 1]426427428def perm_sign(p):429    seen = [False] * len(p)430    s = 1431    for i in range(len(p)):432        if seen[i]:433            continue434        j, c = i, 0435        while not seen[j]:436            seen[j] = True437            j = p[j]438            c += 1439        if c % 2 == 0:440            s = -s441    return s442443444def block_frame(k):445    s_n = 1 << (k - 1)446    h = 1 << (k - 2)447    rows, cols = [], []448    for s in range(h):449        rows += [s, s + h]450        cols += [(2 * s) % s_n, (2 * s + 1) % s_n]451    return rows, cols, perm_sign(rows) * perm_sign(cols)452453454def block_det(m, rows, cols, sign):455    d = 1456    for i in range(0, len(rows), 2):457        d *= m[rows[i]][cols[i]] * m[rows[i + 1]][cols[i + 1]] - m[rows[i]][cols[i + 1]] * m[rows[i + 1]][cols[i]]458        if d == 0:459            return 0460    return sign * d461462463def block_support(k):464    s_n = 1 << (k - 1)465    h = 1 << (k - 2)466    m = transfer((1 << (1 << k)) - 1, 1, k)467    cover = []468    for s in range(h):469        pair = {(2 * s) % s_n, (2 * s + 1) % s_n}470        for r in (s, s + h):471            if {c for c in range(s_n) if m[r][c]} - pair:472                return False473        cover += sorted(pair)474    return sorted(cover) == list(range(s_n))475476# PERRON477478def pari_minpolys(polys, scratch):479    body = (480        "for(i = 1, #V, p = V[i]; rho = vecmax(abs(polroots(p))); "481        "f = factor(p); g = 0; n = 0; w = 0; "482        "for(j = 1, #f~, h = f[j, 1]; if(poldegree(h) > 0, s = polrootsreal(h); "483        "if(#s > 0 && abs(vecmax(s) - rho) < 1e-20, g = h; n = n + 1; "484        "w = vecmax(abs(polroots(h)))))); "485        'print(i, "|", Vec(g), "|", n, "|", rho, "|", w))'486    )487    lines = [f"V = [{', '.join(poly_str(p) for p in polys)}];", body, "quit()"]488    src = scratch / "lemmas.gp"489    src.write_text("\n".join(lines) + "\n")490    out = subprocess.run(491        ["gp", "-q", str(src)],492        capture_output=True,493        text=True,494        stdin=subprocess.DEVNULL,495    )496    if out.returncode != 0:497        sys.exit(out.stderr)498    got = {}499    for line in out.stdout.strip().splitlines():500        idx, vec, hits, rho, wide = line.split("|")501        cs = tuple(int(t) for t in vec.strip().strip("[]").split(","))502        got[int(idx) - 1] = (cs, int(hits), rho.strip(), float(wide) - float(rho))503    return [got[i] for i in range(len(polys))]504505506def verb_lemmas(scratch):507    t0 = time.time()508    keys, dets = {}, {}509    for k in range(1, 5):510        ncodes = 1 << (1 << k)511        seen = set()512        if k < 2:513            for code in range(ncodes):514                seen.add(charpoly(transfer(code, 1, k)))515        else:516            rows, cols, sign = block_frame(k)517            small, agree = 0, 0518            for code in range(ncodes):519                m = transfer(code, 1, k)520                seen.add(charpoly(m))521                d = det_bareiss(m)522                small += d in (-1, 0, 1)523                agree += block_det(m, rows, cols, sign) == d524            dets[k] = (ncodes, small, agree, 1 << (k - 2), block_support(k))525        keys[k] = sorted(seen)526    for k in sorted(dets):527        ncodes, small, agree, blocks, support = dets[k]528        print(529            f"det k={k} rules={ncodes} blocks={blocks} det in [-1,0,1] {small} of {ncodes} "530            f"block product agrees {agree} of {ncodes} support={int(support)}"531        )532    minpolys, radii = [], []533    for k in range(1, 5):534        got = pari_minpolys(keys[k], scratch)535        ones = sum(1 for _, n, _, _ in got if n == 1)536        tight = sum(1 for _, _, _, gap in got if abs(gap) < 1e-20)537        mp = {cs for cs, _, _, _ in got}538        rr = {rho for _, _, rho, _ in got}539        minpolys.append(len(mp))540        radii.append(len(rr))541        print(542            f"perron k={k} charpolys={len(keys[k])} minpolys={len(mp)} distinct_rho={len(rr)} "543            f"unique factor {ones} of {len(got)} no conjugate above rho {tight} of {len(got)}"544        )545    print("minpolys " + ", ".join(str(v) for v in minpolys) + " and distinct rho " + ", ".join(str(v) for v in radii))546    print(f"lemmas {time.time() - t0:.2f}s")547548549if __name__ == "__main__":550    verb = sys.argv[1] if len(sys.argv) > 1 else ""551    if verb == "burnside":552        verb_burnside()553    elif verb == "lemmas":554        verb_lemmas(Path(sys.argv[2]) if len(sys.argv) > 2 else Path(tempfile.gettempdir()))555    else:556        main()