large.py

1import math2import sys3import time4from fractions import Fraction5from functools import lru_cache6from itertools import combinations78import mpmath as mp9import sympy as sp1011CODE = (1, 3, 0, 0)12R = sp.Rational13N = sp.symbols("n")14B = sp.symbols("b", positive=True)15FIVE = {"g": (0, 0), "g+b": (0, -1), "g+1": (1, 0), "g+b-1": (-1, -1), "g+b+1": (1, -1)}1617# DIGIT POLYNOMIALS1819def weights(code):20    w = [0, 0, 0, 0]21    for p in range(8):22        if code >> p & 1:23            w[bin(p).count("1")] += 124    return tuple(w)2526def convolve(p, q):27    out = [0] * (len(p) + len(q) - 1)28    for i, x in enumerate(p):29        if x:30            for j, y in enumerate(q):31                if y:32                    out[i + j] += x * y33    return out3435BASIS = {}3637def digit_polynomial(b, w):38    if b not in BASIS:39        e = [1 - k % 2 for k in range(b)]40        o = [k % 2 for k in range(b)]41        ee, eo, oo = convolve(e, e), convolve(e, o), convolve(o, o)42        BASIS[b] = [convolve(ee, e), convolve(ee, o), convolve(eo, o), convolve(oo, o)]43    return [sum(w[k] * BASIS[b][k][s] for k in range(4)) for s in range(3 * b - 2)]4445def at(p, s):46    return p[s] if 0 <= s < len(p) else 04748def automaton(b, w=CODE):49    p, g = digit_polynomial(b, w), 3 * (b - 1) // 250    return [[at(p, c + g - b * d) for d in (-1, 0, 1)] for c in (-1, 0, 1)]5152def block(b, w=CODE):53    p, g = digit_polynomial(b, w), 3 * (b - 1) // 254    return [[at(p, g), at(p, g + b)], [2 * at(p, g + 1), at(p, g + b - 1) + at(p, g + b + 1)]]5556def fill(b, w):57    n = (b + 1) // 258    return sum(w[k] * n ** (3 - k) * (n - 1) ** k for k in range(4))5960def spectra_block(b):61    if b % 4 == 3:62        return [[3 * (b + 1) * (3 * b - 1) // 16, (b + 1) * (b + 5) // 32], [3 * (b + 1) ** 2 // 8, 3 * (b + 1) ** 2 // 16]]63    return [[(3 * b * b + 6 * b + 7) // 16, 3 * (b - 1) * (b + 3) // 32], [3 * (b - 1) * (3 * b + 5) // 8, (b + 3) ** 2 // 16]]6465def matmul(a, c):66    return [[sum(a[i][k] * c[k][j] for k in range(len(c))) for j in range(len(c[0]))] for i in range(len(a))]6768def word_count(word, letter, centre):69    m = None70    for b in reversed(word):71        x = letter(b)72        m = x if m is None else matmul(m, x)73    return m[centre][centre]7475def brute(word, code):76    side = math.prod(word)77    target = 3 * (side - 1) // 278    keep = [code >> p & 1 for p in range(8)]7980    def ok(x, y, z):81        for b in reversed(word):82            if not keep[(x % b & 1) | (y % b & 1) << 1 | (z % b & 1) << 2]:83                return False84            x, y, z = x // b, y // b, z // b85        return True8687    return sum(ok(x, y, target - x - y) for x in range(side) for y in range(max(0, target - x - side + 1), min(side, target - x + 1)))8889# THE WORDS9091def verb_words():92    print("WORDS: brute-force slice counts against the carry product, finest letter on the left")93    expect = {(3, 5): 60, (5, 3): 72, (3, 5, 7): 2412, (7, 5, 3): 2688, (5, 7): 300}94    bad = []95    for word in [(3, 5), (5, 3), (3, 5, 7), (7, 5, 3), (5, 7), (7, 5), (3, 3, 5), (5, 3, 3), (3, 7, 5), (5, 7, 9), (9, 7, 5), (3, 3, 3)]:96        x = brute(word, 23)97        y = word_count(word, automaton, 1)98        z = word_count(word, block, 0)99        if not x == y == z or expect.get(word, x) != x:100            bad.append(word)101        print("word=%s brute=%d automaton=%d block=%d" % (word, x, y, z))102    for code in (105, 150, 126, 127, 7, 129, 255):103        for word in [(3, 5), (5, 3), (3, 5, 7), (7, 5, 3)]:104            x, y = brute(word, code), word_count(word, lambda b: automaton(b, weights(code)), 1)105            if x != y:106                bad.append((code, word))107        print("code=%d words (3,5) (5,3) (3,5,7) (7,5,3): %s" % (code, [brute(wd, code) for wd in [(3, 5), (5, 3), (3, 5, 7), (7, 5, 3)]]))108    leak = [b for b in range(3, 42, 2) if any(at(digit_polynomial(b, (1, 3, 3, 1)), c + 3 * (b - 1) // 2 - b * d) for c in (-1, 0, 1) for d in (-2, 2))]109    print("no carry leaves abs(c) <= 1 for the full cube at odd b = 3..41: %s" % (not leak))110    print("mismatches=%s" % bad)111112# THE LETTER, SYMBOLIC113114@lru_cache(maxsize=None)115def symbolic_coefficient(k, c, d, parity):116    if (c + d + parity + 1 - k) % 2:117        return sp.Integer(0), 2118    t = (c + 3 * (N - 1) - (2 * N - 1) * d - k) / 2119    total, start = sp.Integer(0), 2120    for size in range(4):121        for sub in combinations(range(3), size):122            x = sp.expand(t - sum(N - (1 if i < k else 0) for i in sub))123            poly = sp.Poly(x, N)124            alpha, beta = poly.coeff_monomial(N), poly.coeff_monomial(1)125            if alpha > 0:126                total += (-1) ** size * (x + 2) * (x + 1) / 2127                start = max(start, int(sp.ceiling(-beta / alpha)))128            else:129                start = max(start, int(sp.floor(-beta / alpha)) + 1)130    return sp.expand(total), start131132def symbolic_block(w, parity):133    coef, start = {}, 2134    for name, (c, d) in FIVE.items():135        total = sp.Integer(0)136        for k in range(4):137            if w[k]:138                x, s0 = symbolic_coefficient(k, c, d, parity)139                total += w[k] * x140                start = max(start, s0)141        coef[name] = sp.expand(total)142    s = sp.Matrix([[coef["g"], coef["g+b"]], [2 * coef["g+1"], coef["g+b-1"] + coef["g+b+1"]]])143    return s, start144145def symbolic_fill(w):146    return sp.expand(sum(w[k] * N ** (3 - k) * (N - 1) ** k for k in range(4)))147148def in_b(expr):149    return sp.expand(expr.subs(N, (B + 1) / 2))150151def parity_masses(w, parity):152    even, odd = R(w[0] + w[2], 4), R(w[1] + w[3], 4)153    return (odd, even) if parity == 0 else (even, odd)154155def limit_letter(u, v):156    return sp.Matrix([[3 * u / 4, v / 8], [3 * v / 2, u / 4]])157158def perron(m):159    lam = max(m.eigenvals(), key=lambda e: sp.N(e))160    right = (m - lam * sp.eye(2)).nullspace()[0]161    left = (m.T - lam * sp.eye(2)).nullspace()[0]162    return sp.radsimp(lam), right, left163164def first_order(a, r):165    lam, right, left = perron(a)166    return sp.radsimp(sp.simplify((left.T * r * right)[0] / (lam * (left.T * right)[0]))), lam167168CLASS = {0: "b = 3 mod 4", 1: "b = 1 mod 4"}169170def expansion(parity):171    s, start = symbolic_block(CODE, parity)172    a = s.applyfunc(lambda e: sp.Poly(in_b(e), B).coeff_monomial(B**2))173    first = s.applyfunc(lambda e: sp.Poly(in_b(e), B).coeff_monomial(B))174    second = s.applyfunc(lambda e: sp.Poly(in_b(e), B).coeff_monomial(1))175    return s, start, a, first, second176177def verb_letter():178    print("LETTER: the carry block of one letter over b^2 as an exact quadratic in t = 1/b, and its limit")179    for parity in (0, 1):180        s, start, a, first, second = expansion(parity)181        bad = []182        for b in range(3, 102, 2):183            n = (b + 1) // 2184            if n % 2 != parity:185                continue186            if n >= start and sp.Matrix(s.subs(N, n)) != sp.Matrix(block(b)):187                bad.append(b)188            if sp.Matrix(spectra_block(b)) != sp.Matrix(block(b)):189                bad.append(("spectra", b))190        print("%s: S_b/b^2 = A + B t + B_2 t^2 = %s + t %s + t^2 %s, extraction stable from b >= %d, mismatches at odd b = 3..101: %s" % (CLASS[parity], a.tolist(), first.tolist(), second.tolist(), 2 * start - 1, bad))191        u, w = parity_masses(CODE, parity)192        q = {0: R(1, 4), 1: R(3, 4)}193        f = {-1: R(1, 8), 0: R(3, 4), 1: R(1, 8)}194        three = all(sp.Poly(in_b(sum(CODE[k] * symbolic_coefficient(k, c, d, parity)[0] for k in range(4))), B).coeff_monomial(B**2) == f[d] * q[(c + d + 1 - parity) % 2] for c in (-1, 0, 1) for d in (-1, 0, 1))195        fb = sp.Poly(in_b(symbolic_fill(CODE)), B)196        print("  M_b[c, c']/b^2 -> f(c') q(c + c' + g) on all nine carries, f = 1/8, 3/4, 1/8, q = 1/4 even, 3/4 odd: %s; fill/b^3 = %s (1 + (%s)/b + O(b^-2))" % (three, fb.coeff_monomial(B**3), fb.coeff_monomial(B**2) / fb.coeff_monomial(B**3)))197        lam, _, _ = perron(a)198        print("  limit = [[3u/4, v/8], [3v/2, u/4]] with u = %s, v = %s: %s; trace %s, det %s, lambda = %s = %.6f" % (u, w, a == limit_letter(u, w), a.trace(), a.det(), lam, float(lam)))199        t = in_b(s.trace())200        dt = in_b(s.det())201        root = sp.factor(t / 2) + sp.sqrt(sp.factor(t**2 / 4 - dt))202        print("  exact root rho_b = %s, disc/16 = %s" % (sp.simplify(root).subs(sp.Abs(B - 1), B - 1), sp.factor(t**2 - 4 * dt) / 16))203        mu, _ = first_order(a, first)204        print("  rho_b/b^2 = lambda (1 + mu/b + O(b^-2)), mu = %s" % mu)205        print("  (log b)(log_b rho_b - log_b fill + 1) = log(2 lambda) + (mu - 3/2)/b + O(b^-2) = %.6f + (%s)/b, sign %s" % (math.log(2 * float(lam)), sp.radsimp(mu - R(3, 2)), "+" if 2 * lam > 1 else "-"))206    a3, a1 = expansion(0)[2], expansion(1)[2]207    lam, _, _ = perron(a1 * a3)208    print("pair: lambda(A_1 A_3) = %s, per letter sqrt = %.6f; [A_3, A_1] = %s" % (lam, math.sqrt(float(lam)), (a3 * a1 - a1 * a3).tolist()))209210# THE DRIFT211212def ink_logs(sides, marks, rows):213    v, scale, side2, out = [1.0, 0.0], 0.0, 0.0, {}214    for i, b in enumerate(sides, 1):215        s = spectra_block(b)216        x = [[s[0][0] / b**2, s[0][1] / b**2], [s[1][0] / b**2, s[1][1] / b**2]]217        if rows:218            v = [v[0] * x[0][0] + v[1] * x[1][0], v[0] * x[0][1] + v[1] * x[1][1]]219        else:220            v = [x[0][0] * v[0] + x[0][1] * v[1], x[1][0] * v[0] + x[1][1] * v[1]]221        top = max(v)222        v = [y / top for y in v]223        scale += math.log(top)224        side2 += 2 * math.log(b)225        if i in marks:226            out[i] = scale + math.log(v[0]) - math.log(0.75 + 0.25 * math.exp(-side2))227    return out228229WORDS = {230    "sides 3, 7, 11, .. coarsest first": (lambda k: 4 * k - 1, False),231    "sides 5, 9, 13, .. coarsest first": (lambda k: 4 * k + 1, False),232    "sides 3, 5, 7, .. coarsest first": (lambda k: 2 * k + 1, False),233    "sides 3, 5, 7, .. finest first": (lambda k: 2 * k + 1, True),234}235236def exponents():237    a3, b3 = expansion(0)[2:4]238    a1, b1 = expansion(1)[2:4]239    mu3, l3 = first_order(a3, b3)240    mu1, l1 = first_order(a1, b1)241    nu, lp = first_order(a1 * a3, b1 * a3 + a1 * b3)242    nu2, _ = first_order(a3 * a1, b3 * a1 + a3 * b1)243    return {"3": (l3, mu3 / 4), "1": (l1, mu1 / 4), "pair": (sp.sqrt(lp), nu / 4), "pair reversed": (sp.sqrt(lp), nu2 / 4)}, (a3, a1, b3, b1)244245def verb_drift():246    print("DRIFT: first-order perturbation of the limiting letter, then local fits at two lengths")247    ex, (a3, a1, b3, b1) = exponents()248    print("first-order term identical in both classes: %s = (3/8) [[1, 1/2], [2, 1]]: %s" % (b3.tolist(), b3 == b1 == R(3, 8) * sp.Matrix([[1, R(1, 2)], [2, 1]])))249    for key, (lam, gamma) in ex.items():250        print("%s: lambda = %s, gamma = %s = %.9f" % (key, sp.nsimplify(lam), sp.radsimp(gamma), float(gamma)))251    marks = [4000, 8000, 16000, 32000]252    keys = ["3", "1", "pair", "pair"]253    for (name, (side, rows)), key in zip(WORDS.items(), keys):254        lam, gamma = map(float, ex[key])255        out = ink_logs([side(k) for k in range(1, marks[-1] + 1)], set(marks), rows)256        est = [(out[2 * m] - out[m] - m * math.log(lam)) / math.log(2) for m in marks[:-1]]257        rich = [2 * y - x for x, y in zip(est, est[1:])]258        print("%s: derived gamma %.6f, local fits L -> 2L at L = 4000, 8000, 16000: %s; Richardson at 8000, 16000: %s" % (name, gamma, ", ".join("%.6f" % e for e in est), ", ".join("%.6f" % e for e in rich)))259260# THE CONSTANTS261262def neville(xs, ys):263    p = list(ys)264    for k in range(1, len(xs)):265        for i in range(len(xs) - 1, k - 1, -1):266            p[i] = (p[i] * xs[i - k] - p[i - 1] * xs[i]) / (xs[i - k] - xs[i])267    return p[-1]268269def constants_along(side, rows, marks, lam, gamma):270    v, scale, side2, out = [mp.mpf(1), mp.mpf(0)], mp.mpf(0), mp.mpf(0), {}271    for i in range(1, max(marks) + 1):272        b = side(i)273        s = spectra_block(b)274        bb = mp.mpf(b) ** 2275        x = [[s[0][0] / bb, s[0][1] / bb], [s[1][0] / bb, s[1][1] / bb]]276        if rows:277            v = [v[0] * x[0][0] + v[1] * x[1][0], v[0] * x[0][1] + v[1] * x[1][1]]278        else:279            v = [x[0][0] * v[0] + x[0][1] * v[1], x[1][0] * v[0] + x[1][1] * v[1]]280        top = max(v)281        v = [y / top for y in v]282        scale += mp.log(top)283        side2 += 2 * mp.log(b)284        if i in marks:285            out[i] = mp.exp(scale + mp.log(v[0]) - mp.log(mp.mpf(3) / 4 + mp.exp(-side2) / 4) - i * mp.log(lam) - gamma * mp.log(i))286    return out287288def blink(a3, a1):289    lam, right, _ = perron(a1 * a3)290    coarse = sp.radsimp((a3 * right)[0] / (sp.sqrt(lam) * right[0]))291    lam, _, left = perron(a3 * a1)292    fine = sp.radsimp((left.T * a3)[0] / (sp.sqrt(lam) * left[0]))293    return coarse, fine294295def verb_constants():296    print("CONSTANTS: ink / (lambda^L L^gamma) at L = 2^j, Neville-extrapolated in 1/L to two depths")297    mp.mp.dps = 40298    ex, (a3, a1, _, _) = exponents()299    keys = ["3", "1", "pair", "pair"]300    found = {}301    for (name, (side, rows)), key in zip(WORDS.items(), keys):302        lam, gamma = (mp.mpf(sp.N(x, 50)) for x in ex[key])303        lengths = [256 * 2**j for j in range(7)]304        for shift, tag in ((0, "even L"), (1, "odd L")):305            if key != "pair" and shift:306                continue307            marks = [m + shift for m in lengths]308            out = constants_along(side, rows, set(marks), lam, gamma)309            xs = [mp.mpf(1) / m for m in marks]310            ys = [out[m] for m in marks]311            deep, shallow = neville(xs, ys), neville(xs[:-1], ys[:-1])312            digits = int(-mp.log10(abs(deep - shallow) / deep)) if deep != shallow else 40313            found[name, tag] = deep314            print("%s, %s: C(L=%d) = %s; extrapolated %s and %s, agreeing to %d digits" % (name, tag, marks[-1], mp.nstr(ys[-1], 12), mp.nstr(deep, 18), mp.nstr(shallow, 18), digits))315    coarse, fine = blink(a3, a1)316    for name, exact in (("sides 3, 5, 7, .. coarsest first", coarse), ("sides 3, 5, 7, .. finest first", fine)):317        ratio = found[name, "odd L"] / found[name, "even L"]318        print("%s: blink C_odd/C_even = %s, derived %s = %s" % (name, mp.nstr(ratio, 15), sp.nsimplify(exact), mp.nstr(mp.mpf(sp.N(exact, 50)), 15)))319    print("row word, finest-first over coarsest-first: even L %s, odd L %s" % (mp.nstr(found["sides 3, 5, 7, .. finest first", "even L"] / found["sides 3, 5, 7, .. coarsest first", "even L"], 15), mp.nstr(found["sides 3, 5, 7, .. finest first", "odd L"] / found["sides 3, 5, 7, .. coarsest first", "odd L"], 15)))320321# THE DESIGNS322323def exact_sign(s, q):324    if s[0][1] * s[1][0] == 0:325        if s[0][0] == 0:326            return None327        return (s[0][0] > q) - (s[0][0] < q)328    t = s[0][0] + s[1][1]329    chi = q * q - t * q + s[0][0] * s[1][1] - s[0][1] * s[1][0]330    if 2 * q < t or chi < 0:331        return 1332    return 0 if chi == 0 else -1333334def last_root(p):335    p = sp.Poly(sp.expand(p), N)336    if p.is_zero or p.degree() == 0:337        return 0338    ends = [hi for (lo, hi), _ in p.intervals()]339    return int(math.ceil(max(ends))) + 1 if ends else 0340341def eventual(w, parity):342    s, start = symbolic_block(w, parity)343    f = symbolic_fill(w)344    if s[0, 1] == 0:345        polys, guard = [(2 * N - 1) * s[0, 0] - f], s[0, 0]346    else:347        t, d = s.trace(), s.det()348        polys, guard = [t * (2 * N - 1) - 2 * f, f**2 - t * f * (2 * N - 1) + d * (2 * N - 1) ** 2], s[0, 1] * s[1, 0]349    roots = max(last_root(p) for p in polys + [guard])350    zeros = [r for r in sp.Poly(sp.expand(guard), N).real_roots()]351    top = max(start, roots)352    n = top + 2 + (top + parity) % 2353    b = 2 * n - 1354    return exact_sign([[int(x) for x in row] for row in s.subs(N, n).tolist()], Fraction(int(f.subs(N, n)), b)), b, start, roots, max(zeros) if zeros else None355356def verb_designs():357    print("DESIGNS: all 256 dim 3 parity designs at infinite side, and the sign of the census exponent against solid dimension minus 1 at every odd base")358    classes = {}359    for code in range(1, 256):360        classes.setdefault(weights(code), []).append(code)361    tally, mism, late, law, reach, flat, guards, blind = {}, [], [], 0, [0, 0, 0], [], {}, set()362    for w, codes in sorted(classes.items()):363        e, o = w[0] + w[2], w[1] + w[3]364        row = []365        for parity in (0, 1):366            u, v = parity_masses(w, parity)367            lam = (2 * u + sp.sqrt(u**2 + 3 * v**2)) / 4368            s, start = symbolic_block(w, parity)369            for b in range(2 * start - 1, 62, 2):370                if ((b + 1) // 2) % 2 == parity and [[int(x) for x in r] for r in s.subs(N, (b + 1) // 2).tolist()] != block(b, w):371                    mism.append((w, b, "form"))372            if s.applyfunc(lambda x: sp.Poly(in_b(x), B).coeff_monomial(B**2)) != limit_letter(u, v):373                mism.append((w, parity, "limit"))374            sign, b_star, start, roots, zero = eventual(w, parity)375            if zero is not None:376                guards[w, parity] = zero377            if s[0, 0] == 0 and s[1, 1] == 0:378                blind.update(codes)379            reach = [max(reach[0], start), max(reach[1], b_star), max(reach[2], roots)]380            predicted = int(sp.sign(u - v))381            if sign != predicted or predicted != int(sp.sign(o - e)) * (1 - 2 * parity):382                mism.append((w, parity, "sign"))383            for b in range(3, b_star + 1, 2):384                if ((b + 1) // 2) % 2 == parity:385                    x = exact_sign(block(b, w), Fraction(fill(b, w), b))386                    if x != sign:387                        late.append((w, b, x, sign))388                    if x != sign and not (x is None and sign == -1):389                        law += len(codes)390            row.append((sign, lam))391        if e == o:392            for b in range(3, 42, 2):393                m = automaton(b, w)394                if any(sum(r) * b != fill(b, w) for r in m):395                    flat.append((w, b))396        key = (row[0][0], row[1][0])397        tally[key] = tally.get(key, 0) + len(codes)398        print("weights=%s e=%d o=%d designs=%d: 8 lambda/(e + o) = %.6f at 3 mod 4, %.6f at 1 mod 4, sign %s" % (w, e, o, len(codes), float(8 * row[0][1] / (e + o)), float(8 * row[1][1] / (e + o)), key))399    print("tally of signs (3 mod 4, 1 mod 4) over codes 1..255: %s" % dict(sorted(tally.items())))400    print("symbolic forms stable from b >= %d; every real root of every sign polynomial below n = %d, b = %d; every odd base up to b = %d checked exactly" % (2 * reach[0] - 1, reach[2], 2 * reach[2] - 1, reach[1]))401    print("largest real root of the case guard, s00 on a triangular class and s01 s10 otherwise: %s at %s" % ((max(guards.values()), max(guards, key=guards.get)) if guards else "none"))402    print("codes whose block has a zero diagonal on a class, count zero at every odd level there: %d" % len(blind))403    print("odd bases where the exact sign differs from the large-side sign: %s (None is an empty slice)" % late)404    print("designs violating sign = sgn(o - e) at 3 mod 4 and sgn(e - o) at 1 mod 4, an empty slice counted below: %d" % law)405    print("e = o designs with an automaton row sum other than fill/b at odd b = 3..41: %s" % flat)406    print("mismatches=%s" % mism)407408# THE CLOSED FORM AT SIDES 3 MOD 4409410def verb_closed():411    print("CLOSED: sides 3, 7, 11, .. coarsest first, the generating function and the constant")412    k, z = sp.symbols("k z")413    a = sp.Matrix([[18, 1], [12, 6]])414    c = sp.Matrix([[-6, 1], [0, 0]])415    b = 4 * k - 1416    letter = sp.Matrix([[3 * (b + 1) * (3 * b - 1) / 16, (b + 1) * (b + 5) / 32], [3 * (b + 1) ** 2 / 8, 3 * (b + 1) ** 2 / 16]])417    print("S_(4k-1) = (k/2)(k K_1 + K_0) with K_1 = %s, K_0 = %s: %s" % (a.tolist(), c.tolist(), sp.expand(letter - k / 2 * (k * a + c)) == sp.zeros(2)))418    d = 1 - 24 * z + 96 * z**2419    p = sp.Matrix([1 - 6 * z, 12 * z])420    lhs = (sp.eye(2) - z * a) * (d * p.diff(z) - R(3, 4) * d.diff(z) * p)421    rhs = d * (a + c) * p422    print("W = (1 - 24z + 96z^2)^(-3/4) (1 - 6z, 12z) solves (1 - z K_1) W' = (K_1 + K_0) W, W(0) = e_0: %s" % (sp.expand(lhs - rhs) == sp.zeros(2, 1)))423    y = [Fraction(1), Fraction(18)]424    for m in range(1, 60):425        y.append(((24 * m + 18) * y[m] - (96 * m + 48) * y[m - 1]) / (m + 1))426    v, bad = [1, 0], []427    for L in range(1, 61):428        s = spectra_block(4 * L - 1)429        v = [s[0][0] * v[0] + s[0][1] * v[1], s[1][0] * v[0] + s[1][1] * v[1]]430        series = y[L] - 6 * y[L - 1]431        if Fraction(math.factorial(L) ** 2, 2**L) * series != v[0]:432            bad.append(L)433        if L <= 4:434            print("L=%d count=%d" % (L, v[0]))435    print("count_L = (L!)^2 2^(-L) [z^L] (1 - 6z)(1 - 24z + 96z^2)^(-3/4) at L = 1..60: mismatches=%s" % bad)436    mp.mp.dps = 40437    exact = mp.gamma(mp.mpf(3) / 4) * (1 + mp.sqrt(3)) / (3 * (mp.sqrt(3) - 1) ** (mp.mpf(3) / 4))438    lam = (3 + mp.sqrt(3)) / 8439    marks = [256 * 2**j for j in range(7)]440    out = constants_along(lambda i: 4 * i - 1, False, set(marks), lam, mp.mpf(1) / 4)441    print("C_3 = Gamma(3/4)(1 + sqrt(3))/(3 (sqrt(3) - 1)^(3/4)) = %s; extrapolated %s" % (mp.nstr(exact, 25), mp.nstr(neville([mp.mpf(1) / m for m in marks], [out[m] for m in marks]), 20)))442443# RUN444445VERBS = {"words": verb_words, "letter": verb_letter, "drift": verb_drift, "constants": verb_constants, "closed": verb_closed, "designs": verb_designs}446447if __name__ == "__main__":448    for name in sys.argv[1:] or list(VERBS):449        start = time.time()450        VERBS[name]()451        print("%s: %.1fs" % (name, time.time() - start))