compose.py

1import sys2import time3from fractions import Fraction4from itertools import product56import sympy as sp78BASES = list(range(3, 26, 2))91011def digit_polynomial(b):12    p = [0] * (3 * (b - 1) + 1)13    for v in product(range(b), repeat=3):14        if sum(d % 2 for d in v) <= 1:15            p[sum(v)] += 116    return p171819def generating_polynomial(b):20    n = (b + 1) // 221    e = [1 if k % 2 == 0 else 0 for k in range(b)]22    o = [1 if k % 2 == 1 else 0 for k in range(b)]23    ee = convolve(e, e)24    return [x + 3 * y for x, y in zip(convolve(ee, e), convolve(ee, o))]252627def convolve(p, q):28    out = [0] * (len(p) + len(q) - 1)29    for i, x in enumerate(p):30        if x:31            for j, y in enumerate(q):32                out[i + j] += x * y33    return out343536def coefficient(p, s):37    return p[s] if 0 <= s < len(p) else 0383940def automaton(b, p):41    g = 3 * (b - 1) // 242    return {(c, d): coefficient(p, c + g - b * d) for c in (-1, 0, 1) for d in (-1, 0, 1)}434445def leak(b, p):46    g = 3 * (b - 1) // 247    return sum(coefficient(p, c + g - b * d) for c in (-1, 0, 1) for d in (-2, 2))484950def even_block(m):51    return [[m[0, 0], m[0, 1]], [m[1, 0] + m[-1, 0], m[1, 1] + m[-1, 1]]]525354def power_column(m, levels):55    vec = {-1: 0, 0: 1, 1: 0}56    out = [(vec[0], vec[1] + vec[-1])]57    for _ in range(levels):58        vec = {c: sum(m[c, d] * vec[d] for d in (-1, 0, 1)) for c in (-1, 0, 1)}59        out.append((vec[0], vec[1] + vec[-1]))60    return out616263def layer_census(b, p, level):64    dist = {0: 1}65    for k in range(level):66        shift = b**k67        nxt = {}68        for s0, v0 in dist.items():69            for s1, v1 in enumerate(p):70                if v1:71                    nxt[s0 + s1 * shift] = nxt.get(s0 + s1 * shift, 0) + v0 * v172        dist = nxt73    h = 3 * (b**level - 1) // 274    return dist.get(h, 0), dist.get(h - 1, 0) + dist.get(h + 1, 0)757677def spectra_closed_form(b):78    if b % 4 == 3:79        return [80            [Fraction(3 * (b + 1) * (3 * b - 1), 16), Fraction((b + 1) * (b + 5), 32)],81            [Fraction(3 * (b + 1) ** 2, 8), Fraction(3 * (b + 1) ** 2, 16)],82        ]83    return [84        [Fraction(3 * b * b + 6 * b + 7, 16), Fraction(3 * (b - 1) * (b + 3), 32)],85        [Fraction(3 * (b - 1) * (3 * b + 5), 8), Fraction((b + 3) ** 2, 16)],86    ]878889def verb_block():90    print("BLOCK: even carry block from the automaton against spectra's closed forms and a layer census")91    bad = []92    for b in BASES:93        p = digit_polynomial(b)94        if p != generating_polynomial(b):95            bad.append((b, "polynomial"))96        m = automaton(b, p)97        block = even_block(m)98        if leak(b, p):99            bad.append((b, "leak"))100        if [[Fraction(x) for x in row] for row in block] != spectra_closed_form(b):101            bad.append((b, "form"))102        counts = power_column(m, 8)103        for k in range(8):104            h, t = counts[k]105            if counts[k + 1] != (block[0][0] * h + block[0][1] * t, block[1][0] * h + block[1][1] * t):106                bad.append((b, "step", k))107        top = {3: 6, 5: 5, 7: 4, 9: 4}.get(b, 3)108        for k in range(top + 1):109            if layer_census(b, p, k) != counts[k]:110                bad.append((b, "census", k))111        print("b=%d g=%d block=%s hexagons=%s" % (b, 3 * (b - 1) // 2, block, [h for h, _ in counts[:5]]))112    print("odd bases 3..25: %d bases, mismatches=%s" % (len(BASES), bad))113114115def binomial2(x):116    return x * (x - 1) / 2117118119def extraction(kind, n, idx):120    if kind == "A":121        terms = [(1, 0), (-3, n), (3, 2 * n), (-1, 3 * n)]122    else:123        terms = [(1, 0), (-2, n), (1, 2 * n), (-1, n - 1), (2, 2 * n - 1), (-1, 3 * n - 1)]124    return [(sign, idx - shift) for sign, shift in terms]125126127def settle(expr, n):128    poly = sp.Poly(sp.expand(expr), n)129    alpha, beta = poly.coeff_monomial(n), poly.coeff_monomial(1)130    if alpha > 0:131        return True, sp.ceiling(-beta / alpha)132    if alpha < 0:133        return False, sp.floor(-beta / alpha) + 1134    return beta >= 0, 0135136137def closed_coefficient(parity, n, m, s_of_n):138    nn = 2 * m + parity139    s = s_of_n.subs(n, nn)140    even = sp.Poly(sp.expand(s), m).coeff_monomial(m) % 2 == 0 and int(s.subs(m, 0)) % 2 == 0141    idx = (s / 2) if even else ((s - 1) / 2)142    kind = "A" if even else "B"143    total, threshold = 0, 2144    for sign, arg in extraction(kind, nn, idx):145        active, since = settle(arg, m)146        threshold = max(threshold, 2 * since + parity)147        if active:148            total += sign * binomial2(arg + 2)149    if not even:150        total *= 3151    return sp.expand(total), threshold152153154ENTRIES = {155    "P[g]": lambda n: 3 * (n - 1),156    "P[g+1]": lambda n: 3 * (n - 1) + 1,157    "P[g+b]": lambda n: 3 * (n - 1) + 2 * n - 1,158    "P[g+b-1]": lambda n: 3 * (n - 1) + 2 * n - 2,159    "P[g+b+1]": lambda n: 3 * (n - 1) + 2 * n,160}161162163def derive_forms():164    n, m, b = sp.symbols("n m b")165    forms = {}166    for parity in (0, 1):167        for name, s_of_n in ENTRIES.items():168            poly_m, threshold = closed_coefficient(parity, n, m, s_of_n(n))169            poly_b = sp.factor(sp.expand(poly_m.subs(m, ((b + 1) / 2 - parity) / 2)))170            forms[(parity, name)] = (poly_b, threshold)171    return forms172173174def evaluate(expr, b_value):175    b = sp.Symbol("b")176    return sp.Rational(expr.subs(b, b_value))177178179def verb_forms():180    print("FORMS: the five coefficients as polynomials in b, by coefficient extraction on E^2 (E + 3 O)")181    forms = derive_forms()182    b = sp.Symbol("b")183    for parity, klass in ((0, "b = 3 mod 4"), (1, "b = 1 mod 4")):184        for name in ENTRIES:185            poly_b, threshold = forms[(parity, name)]186            print("%s: %s = %s, extraction stable from b >= %d" % (klass, name, poly_b, 2 * threshold - 1))187    bad = []188    for b_value in range(3, 102, 2):189        p = generating_polynomial(b_value)190        parity = 0 if b_value % 4 == 3 else 1191        for name, s_of_n in ENTRIES.items():192            s = s_of_n((b_value + 1) // 2)193            if evaluate(forms[(parity, name)][0], b_value) != coefficient(p, s):194                bad.append((b_value, name))195    print("closed forms against the digit polynomial at odd b = 3..101: mismatches=%s, every extraction threshold is at most b = %d" % (bad, max(2 * t - 1 for _, t in forms.values())))196    block = {}197    for parity in (0, 1):198        a = forms[(parity, "P[g]")][0]199        c = forms[(parity, "P[g+b]")][0]200        d = 2 * forms[(parity, "P[g+1]")][0]201        e = forms[(parity, "P[g+b-1]")][0] + forms[(parity, "P[g+b+1]")][0]202        block[parity] = [[sp.factor(a), sp.factor(c)], [sp.factor(d), sp.factor(e)]]203    spectra = {204        0: [[3 * (b + 1) * (3 * b - 1) / 16, (b + 1) * (b + 5) / 32], [3 * (b + 1) ** 2 / 8, 3 * (b + 1) ** 2 / 16]],205        1: [[(3 * b * b + 6 * b + 7) / 16, 3 * (b - 1) * (b + 3) / 32], [3 * (b - 1) * (3 * b + 5) / 8, (b + 3) ** 2 / 16]],206    }207    for parity, klass in ((0, "b = 3 mod 4"), (1, "b = 1 mod 4")):208        same = all(sp.expand(block[parity][i][j] - spectra[parity][i][j]) == 0 for i in range(2) for j in range(2))209        tr = sp.factor(block[parity][0][0] + block[parity][1][1])210        det = sp.factor(block[parity][0][0] * block[parity][1][1] - block[parity][0][1] * block[parity][1][0])211        print("%s: block %s equals spectra's closed form: %s, trace %s, det %s" % (klass, block[parity], same, tr, det))212    return block213214215def positive_from(poly, n, start):216    p = sp.Poly(sp.expand(poly), n)217    roots = p.intervals()218    top = max((hi for (lo, hi), _ in roots), default=sp.Integer(-10**9))219    value = p.eval(start)220    return top < start and value > 0, top, value221222223def verb_split():224    print("SPLIT: the sign of rho_b - fill/b as a polynomial inequality in n = (b + 1)/2")225    b, n = sp.symbols("b n")226    block = verb_forms()227    fill = n**2 * (4 * n - 3)228    q = fill / (2 * n - 1)229    for parity, klass, start, want in ((0, "b = 3 mod 4", 2, "above"), (1, "b = 1 mod 4", 3, "below")):230        blk = [[e.subs(b, 2 * n - 1) for e in row] for row in block[parity]]231        tr = sp.expand(blk[0][0] + blk[1][1])232        det = sp.expand(blk[0][0] * blk[1][1] - blk[0][1] * blk[1][0])233        disc = sp.expand(tr**2 - 4 * det)234        u = sp.expand(sp.cancel((tr - 2 * q) * (2 * n - 1)))235        chi = sp.expand(sp.cancel((q**2 - tr * q + det) * (2 * n - 1) ** 2))236        ok_disc, top_disc, _ = positive_from(disc, n, start)237        ok_u, top_u, _ = positive_from(-u, n, start)238        sign = -1 if want == "above" else 1239        ok_chi, top_chi, val_chi = positive_from(sign * chi, n, start)240        print("%s: trace %s, det %s, n >= %d" % (klass, sp.factor(tr), sp.factor(det), start))241        print("  disc > 0: %s (largest root bound %s)" % (ok_disc, top_disc))242        print("  (2n-1)(tr - 2q) = %s < 0: %s (largest root bound %s)" % (sp.factor(u), ok_u, top_u))243        print("  (2n-1)^2 chi(q) = %s, needs sign %+d: %s (largest root bound %s, value at n=%d is %s)" % (sp.factor(chi), sign, ok_chi, top_chi, start, sign * val_chi))244        print("  hence rho_b %s fill/b for every b = %s: %s" % (">" if want == "above" else "<", klass[4:], ok_disc and ok_u and ok_chi))245246247def verb_classes():248    print("CLASSES: residue-class sums of P_b are the row sums of the automaton, and P_b(w) = -2/(1 + w)^3 at w^b = 1, w != 1")249    w = sp.Symbol("w")250    bad = []251    for b in BASES:252        p = generating_polynomial(b)253        poly = sum(c * w**s for s, c in enumerate(p))254        rem = sp.rem(sp.expand(poly * (1 + w) ** 3 + 2), sum(w**k for k in range(b)), w)255        if rem != 0:256            bad.append((b, "root"))257        g = 3 * (b - 1) // 2258        m = automaton(b, p)259        for c in (-1, 0, 1):260            row = sum(m[c, d] for d in (-1, 0, 1))261            klass = sum(v for s, v in enumerate(p) if (s - g - c) % b == 0)262            if row != klass:263                bad.append((b, "row", c))264    print("odd b = 3..25: mismatches=%s" % bad)265    b = sp.Symbol("b")266    forms = derive_forms()267    for parity, klass in ((0, "b = 3 mod 4"), (1, "b = 1 mod 4")):268        f = {name: forms[(parity, name)][0] for name in ENTRIES}269        mid = sp.factor(f["P[g]"] + 2 * f["P[g+b]"])270        side = sp.factor(f["P[g+1]"] + f["P[g+b-1]"] + f["P[g+b+1]"])271        wrong = []272        for b_value in range(3 + parity * 2, 102, 4):273            p = generating_polynomial(b_value)274            g = 3 * (b_value - 1) // 2275            sums = [sum(v for s, v in enumerate(p) if (s - g - c) % b_value == 0) for c in (0, 1, -1)]276            if sums != [evaluate(mid, b_value), evaluate(side, b_value), evaluate(side, b_value)]:277                wrong.append(b_value)278        print("%s: class of g sums to %s, classes of g +- 1 sum to %s each, mismatches to b = 101: %s" % (klass, mid, side, wrong))279280281VERBS = {"block": verb_block, "forms": verb_forms, "split": verb_split, "classes": verb_classes}282283284def main():285    names = sys.argv[1:] or ["block", "split", "classes"]286    for name in names:287        start = time.time()288        VERBS[name]()289        print("verb %s took %.1fs" % (name, time.time() - start))290        print("")291292293main()