slice_grammar.py

1from fractions import Fraction23import numpy as np4from mpmath import mp, mpf5from mpmath import log as mlog6from mpmath import sqrt as msqrt78mp.dps = 4091011def digit_sums(b, rule):12    if rule == "odd":13        marked = {d for d in range(b) if d % 2 == 1}14    else:15        marked = {(b - 1) // 2}16    plain = np.zeros(b, dtype=np.int64)17    special = np.zeros(b, dtype=np.int64)18    for d in range(b):19        if d in marked:20            special[d] = 121        else:22            plain[d] = 123    square = np.convolve(plain, plain)24    total = np.convolve(square, plain) + 3 * np.convolve(special, square)25    return [int(v) for v in total]262728def at(c, i):29    if i < 0 or i >= len(c):30        return 031    return c[i]323334def transfer(b, c):35    row_mid = (at(c, (3 * b - 3) // 2), at(c, (3 * b - 1) // 2), at(c, (3 * b - 5) // 2))36    row_low = (at(c, (b - 3) // 2), at(c, (b - 1) // 2), at(c, (b - 5) // 2))37    row_high = (at(c, (5 * b - 3) // 2), at(c, (5 * b - 1) // 2), at(c, (5 * b - 5) // 2))38    return row_mid, row_low, row_high394041def apply_rows(rows, v):42    return tuple(r[0] * v[0] + r[1] * v[1] + r[2] * v[2] for r in rows)434445def tile_counts(b, c, levels):46    rows = transfer(b, c)47    states = [(0, 1, 0), (1, 0, 0), (0, 0, 1)]48    hexes = [states[1][0]]49    tris = [states[0][0] + states[2][0]]50    for _ in range(levels):51        states = [apply_rows(rows, v) for v in states]52        hexes.append(states[1][0])53        tris.append(states[0][0] + states[2][0])54    return hexes, tris555657def layer_counts(b, c, level):58    dist = {0: 1}59    nz = [(s, v) for s, v in enumerate(c) if v]60    for k in range(level):61        shift = b ** k62        nxt = {}63        for s0, v0 in dist.items():64            for s1, v1 in nz:65                key = s0 + s1 * shift66                nxt[key] = nxt.get(key, 0) + v0 * v167        dist = nxt68    m = (3 * b ** level - 1) // 269    return dist.get(m - 1, 0), dist.get(m, 0) + dist.get(m - 2, 0)707172def substitution(b, c, levels=9):73    hexes, tris = tile_counts(b, c, levels)74    m00, m10 = hexes[1], tris[1]75    denom = tris[1]76    m01 = Fraction(hexes[2] - m00 * hexes[1], denom)77    m11 = Fraction(tris[2] - m10 * hexes[1], denom)78    mat = [[Fraction(m00), m01], [Fraction(m10), m11]]79    ok = all(x.denominator == 1 for row in mat for x in row)80    for n in range(1, levels + 1):81        ph = mat[0][0] * hexes[n - 1] + mat[0][1] * tris[n - 1]82        pt = mat[1][0] * hexes[n - 1] + mat[1][1] * tris[n - 1]83        if ph != hexes[n] or pt != tris[n]:84            ok = False85    ints = [[int(x) for x in row] for row in mat] if ok else None86    return ints, hexes, tris, ok878889def spectral(tr, det, b):90    lam = (mpf(tr) + msqrt(mpf(tr) ** 2 - 4 * det)) / 291    return lam, mlog(lam) / mlog(b)929394def side(tr, det, fill, b):95    q = Fraction(fill, b)96    disc = Fraction(tr * tr - 4 * det)97    u = Fraction(tr) - 2 * q98    if u > 0:99        return 1100    if u == 0:101        return 1 if disc > 0 else 0102    val = disc - u * u103    return 1 if val > 0 else (0 if val == 0 else -1)104105106def closed_form(b):107    if b % 4 == 3:108        return [109            [Fraction(3 * (b + 1) * (3 * b - 1), 16), Fraction((b + 1) * (b + 5), 32)],110            [Fraction(3 * (b + 1) ** 2, 8), Fraction(3 * (b + 1) ** 2, 16)],111        ]112    return [113        [Fraction(3 * b * b + 6 * b + 7, 16), Fraction(3 * (b - 1) * (b + 3), 32)],114        [Fraction(3 * (b - 1) * (3 * b + 5), 8), Fraction((b + 3) ** 2, 16)],115    ]116117118def rule_text(tr, det):119    return "x%d %s%d" % (tr, "+" if -det >= 0 else "-", abs(det))120121122def report(b, rule, levels=9):123    c = digit_sums(b, rule)124    fill = sum(c)125    mat, hexes, tris, ok = substitution(b, c, levels)126    tr = mat[0][0] + mat[1][1]127    det = mat[0][0] * mat[1][1] - mat[0][1] * mat[1][0]128    lam, dim = spectral(tr, det, b)129    solid = mlog(mpf(fill)) / mlog(b)130    return {131        "b": b,132        "rule": rule,133        "fill": fill,134        "cells": b ** 3,135        "matrix": mat,136        "closes": ok,137        "trace": tr,138        "det": det,139        "recurrence": rule_text(tr, det),140        "lam": lam,141        "dim": dim,142        "solid": solid,143        "minus_one": solid - 1,144        "hexes": hexes,145        "tris": tris,146        "side": side(tr, det, fill, b),147    }148149150def main():151    print("SLICE GRAMMAR, ODD BASES, AT MOST ONE ODD COORDINATE")152    bases = [3, 5, 7, 9]153    rows = {b: report(b, "odd") for b in bases}154    for b in bases:155        r = rows[b]156        print(157            "b=%d fill=%d/%d matrix=%s closes=%s rule=%s dim_slice=%.4f d-1=%.4f"158            % (159                b,160                r["fill"],161                r["cells"],162                r["matrix"],163                r["closes"],164                r["recurrence"],165                float(r["dim"]),166                float(r["minus_one"]),167            )168        )169    print(170        "four rules: %s"171        % " / ".join(rows[b]["recurrence"] for b in bases)172    )173    print(174        "dimensions: %s"175        % " / ".join("%.4f" % float(rows[b]["dim"]) for b in bases)176    )177    print(178        "d-1: %s" % " / ".join("%.4f" % float(rows[b]["minus_one"]) for b in bases)179    )180181    print("")182    print("BASE THREE EXTERNAL TARGET")183    print("matrix %s" % rows[3]["matrix"])184    print("hexagons %s" % ", ".join(str(v) for v in rows[3]["hexes"][:7]))185    print("triangles %s" % ", ".join(str(v) for v in rows[3]["tris"][:7]))186187    print("")188    print("INDEPENDENT LAYER CENSUS, NO TILE REDUCTION")189    for b, top in ((3, 6), (5, 5), (7, 4), (9, 4)):190        c = digit_sums(b, "odd")191        hexes, tris = tile_counts(b, c, top)192        agree = True193        for n in range(top + 1):194            dh, dt = layer_counts(b, c, n)195            if (dh, dt) != (hexes[n], tris[n]):196                agree = False197        print("b=%d levels 0..%d agree=%s" % (b, top, agree))198199    print("")200    print("CLOSED FORM AGAINST CENSUS")201    bad = []202    for b in range(3, 22, 2):203        r = report(b, "odd")204        cf = closed_form(b)205        got = [[Fraction(x) for x in row] for row in r["matrix"]]206        if cf != got or not r["closes"]:207            bad.append(b)208    print("odd bases 3..21 matched=%d mismatched=%s" % (10 - len(bad), bad))209210    print("")211    print("MOD FOUR SPLIT")212    wrong = []213    tested = 0214    for b in range(3, 402, 2):215        r = report(b, "odd", levels=5)216        tested += 1217        want = 1 if b % 4 == 3 else -1218        if r["side"] != want or not r["closes"]:219            wrong.append(b)220    print("odd bases 3..401 tested=%d exceptions=%s" % (tested, wrong))221222    print("")223    print("BASE FIVE, THE TWO DIGIT RULES")224    odd5 = rows[5]225    mid5 = report(5, "middle")226    print(227        "middle-digit rule fills %d of %d dim_slice=%.6f d-1=%.6f excess=%+.3e"228        % (229            mid5["fill"],230            mid5["cells"],231            float(mid5["dim"]),232            float(mid5["minus_one"]),233            float(mid5["dim"] - mid5["minus_one"]),234        )235    )236    print(237        "odd-coordinate rule fills %d of %d dim_slice=%.4f d-1=%.4f excess=%+.3e"238        % (239            odd5["fill"],240            odd5["cells"],241            float(odd5["dim"]),242            float(odd5["minus_one"]),243            float(odd5["dim"] - odd5["minus_one"]),244        )245    )246    print("middle-digit matrix %s rule %s" % (mid5["matrix"], mid5["recurrence"]))247248249main()