returns.py
52.3 kB · python · 1215 lines
1import argparse2import time3from fractions import Fraction4from itertools import combinations5from math import comb, exp, gcd, log67MINPOLY = ([-1, 1], [-2, 1], [-3, 1], [-6, 1], [45, -15, 1], [90, -26, 1],8 [-3402, 945, -63, 1], [1134, -99, 1], [1299078, -293787, 16065, -255, 1],9 [-96228, 17469, -392, 1])1011def up(x, d):12 f = 10 ** d13 n = int(x * f)14 return (n + 1 if n < x * f else n) / f1516def down(x, d):17 f = 10 ** d18 n = int(x * f)19 return (n - 1 if n > x * f else n) / f2021# BAND AUTOMATON2223def succ(j, z1, z2):24 if j % 3 == 0:25 return (j // 3, (j + z1) // 3)26 if (j - z2) % 3 == 0:27 return ((j - z2) // 3,)28 return ()2930def first_return(z1, z2):31 start = z1 // 332 seen = {start}33 frontier = [start]34 n = 135 while frontier:36 n += 137 nxt = []38 for j in frontier:39 for t in succ(j, z1, z2):40 if t == 0:41 return n42 if t not in seen:43 seen.add(t)44 nxt.append(t)45 frontier = nxt46 return 04748def band(z1, z2):49 return -((z2 - 1) // 2), (z1 - 1) // 25051def degree_profile(z1, z2):52 lo, hi = band(z1, z2)53 cnt = [0, 0, 0]54 edges = 055 for j in range(lo, hi + 1):56 d = len(succ(j, z1, z2))57 cnt[d] += 158 edges += d59 for t in succ(j, z1, z2):60 if not lo <= t <= hi:61 raise ValueError("the band is not invariant at " + str((z1, z2, j)))62 return cnt, edges, hi - lo + 16364# COLUMN TRANSFER6566def slot_profile(k, n):67 return [(n - r + k - 1) // k for r in range(k)]6869def column_count(k, n, prim=True, twist=1):70 if n < k:71 return 072 w = (3 ** k - 1) // 273 slots = slot_profile(k, n)74 top = sum(s * 3 ** r for r, s in enumerate(slots))75 total = 076 for j in range(top // w + 1):77 M = j * w78 cur = {0: 1}79 for r in range(k):80 mr = M // 3 ** r % 381 lo = 1 if prim and r == 0 else 082 nxt = {}83 for t, v in cur.items():84 for c in range(lo, slots[r] + 1):85 if (t + c - mr) % 3:86 continue87 u = (t + c - mr) // 388 nxt[u] = nxt.get(u, 0) + v * comb(slots[r] - lo, c - lo) * twist ** c89 cur = nxt90 total += cur.get(M // 3 ** k, 0)91 return total9293def lift_count(k, n, prim=True):94 return column_count(k, n, prim) - (0 if prim else 1)9596def model_returns(k, N, cap=40):97 w = (3 ** k - 1) // 298 out = []99 lo = hi = 0.0100 prev = 0101 for n in range(k, min(N, cap * k) + 1):102 a = column_count(k, n, True, 2)103 d = a - prev104 prev = a105 if d:106 lo += d * w / 3 ** n107 hi += d * w / 3 ** (n - 1)108 out.append((n, lo, hi))109 n = out[-1][0]110 while n < N:111 n += 1112 lo += 4.0 / 9113 hi += 4.0 / 3114 out.append((n, lo, hi))115 return out116117def exact_L(k):118 acc = 0.0119 for s in range(1 << (k - 1)):120 a = sum(3 ** (i + 1) for i in range(k - 1) if s >> i & 1)121 acc += 1.0 / (1 + 2 * a)122 return acc123124def brute_lift_count(k, n, prim=True):125 w = (3 ** k - 1) // 2126 pw = [3 ** i for i in range(n)]127 out = 0128 for s in range(1, 1 << n):129 if prim and not s & 1:130 continue131 v = sum(pw[i] for i in range(n) if s >> i & 1)132 if v % w == 0:133 out += 1134 return out135136# THE FIXED MATRIX AT n = bk137138def block_head(b, j, K=64):139 M = j * (3 ** K - 1) // 2140 d = [M // 3 ** r % 3 for r in range(K)]141 r0 = K - 1142 while r0 > 1 and d[r0 - 1] == d[K - 1]:143 r0 -= 1144 return d[:r0], d[K - 1], j // 2 - 1 if j and not j % 2 else j // 2145146def block_pieces(b, j):147 head, bulk, fin = block_head(b, j)148 S = b // 2 + 1149 cur = {0: 1}150 for r, mr in enumerate(head):151 lo = 1 if r == 0 else 0152 nxt = {}153 for t, v in cur.items():154 for c in range(lo, b + 1):155 if (t + c - mr) % 3:156 continue157 u = (t + c - mr) // 3158 nxt[u] = nxt.get(u, 0) + v * comb(b - lo, c - lo)159 cur = nxt160 A = [[0] * S for _ in range(S)]161 for t in range(S):162 for c in range(b + 1):163 if (t + c - bulk) % 3:164 continue165 u = (t + c - bulk) // 3166 if u < S:167 A[u][t] += comb(b, c)168 elif comb(b, c):169 raise ValueError("the carry leaves the band at " + str((b, j, t, c)))170 return A, [cur.get(t, 0) for t in range(S)], fin, len(head)171172def block_count(b, k):173 tot = 0174 for j in range(b + 1):175 A, h, f, r0 = block_pieces(b, j)176 if k < r0:177 return None178 for _ in range(k - r0):179 h = [sum(A[i][t] * h[t] for t in range(len(h))) for i in range(len(A))]180 if f < len(h):181 tot += h[f]182 return tot183184def block_closure(A, s):185 seen, front = {s}, [s]186 while front:187 nxt = []188 for t in front:189 for u in range(len(A)):190 if A[u][t] and u not in seen:191 seen.add(u)192 nxt.append(u)193 front = nxt194 return seen195196def block_rho(A, steps=600):197 v = [1.0] * len(A)198 g = 0.0199 for i in range(steps):200 v = [sum(A[u][t] * v[t] for t in range(len(A))) for u in range(len(A))]201 s = sum(v)202 if s <= 0:203 return 0.0204 v = [x / s for x in v]205 if i >= steps // 2:206 g += log(s)207 return exp(g / (steps - steps // 2))208209def block_bracket(b):210 dev, r0s, rhos = set(), [], []211 for j in range(b + 1):212 A, h, f, r0 = block_pieces(b, j)213 r0s.append(r0)214 rhos.append(block_rho(A))215 for t in range(len(A)):216 dev.add(Fraction(3 * sum(A[u][t] for u in range(len(A))) - 2 ** b, 3))217 par = {Fraction(-1, 3), Fraction(2, 3)} if not b % 2 else {Fraction(-2, 3), Fraction(1, 3)}218 if not dev <= par:219 raise ValueError("a column sum leaves the parity pair at b = " + str(b))220 A, h, f, _ = block_pieces(b, 1)221 C = block_closure(A, 0)222 if f or not h[0] or not A[0][0]:223 raise ValueError("the j = 1 sector does not start and end at the carry 0 at b = " + str(b))224 for t in sorted(C):225 if sum(A[u][t] for u in C) != sum(A[u][t] for u in range(len(A))):226 raise ValueError("the closure of the carry 0 is not closed at b = " + str(b))227 if 0 not in block_closure(A, t):228 raise ValueError("the carry 0 is not reachable from " + str(t) + " at b = " + str(b))229 cs = [Fraction(3 * sum(A[u][t] for u in C) - 2 ** b, 3) for t in sorted(C)]230 return sorted(dev), min(cs), max(cs), len(C), r0s, max(rhos)231232# EXACT ALGEBRA233234def poly_div(P, Q):235 P = list(P)236 n, m = len(P) - 1, len(Q) - 1237 if n < m:238 return None239 out = [0] * (n - m + 1)240 for i in range(n - m, -1, -1):241 if P[i + m] % Q[m]:242 return None243 out[i] = P[i + m] // Q[m]244 for t in range(m + 1):245 P[i + t] -= out[i] * Q[t]246 return out if not any(P[:m]) else None247248def poly_at(P, x):249 v = 0250 for c in reversed(P):251 v = v * x + c252 return v253254def poly_show(P):255 out = ""256 for i in range(len(P) - 1, -1, -1):257 if not P[i]:258 continue259 s = "-" if P[i] < 0 else "+"260 a = abs(P[i])261 t = str(a) if a != 1 or not i else ""262 t += "x^" + str(i) if i > 1 else ("x" if i == 1 else "")263 out += (t if not out and s == "+" else " " + s + " " + t)264 return out or "0"265266def poly_roots(P, R):267 n = len(P) - 1268 q = [P[i] / R ** (n - i) for i in range(n + 1)]269 z = [(0.4 + 0.9j) ** i for i in range(n)]270 for _ in range(4000):271 mv = 0.0272 for i in range(n):273 d = 1.0 + 0j274 for t in range(n):275 if t != i:276 d *= z[i] - z[t]277 v = 0j278 for c in reversed(q):279 v = v * z[i] + c280 s = v / d281 z[i] -= s282 mv = max(mv, abs(s))283 if mv < 1e-16:284 break285 return [w * R for w in z]286287def poly_minfactor(P, rs, lam):288 rest = [r for r in rs if r is not lam]289 for d in range(1, len(P)):290 for S in combinations(rest, d - 1):291 f = [1.0 + 0j]292 for r in (lam,) + S:293 g = [0j] * (len(f) + 1)294 for i, c in enumerate(f):295 g[i] -= c * r296 g[i + 1] += c297 f = g298 if max(abs(c.imag) for c in f) > 0.4:299 continue300 g = [round(c.real) for c in f]301 if max(abs(g[i] - f[i].real) for i in range(len(g))) > 0.45:302 continue303 if poly_div(P, g) is not None:304 return g305 return None306307def poly_newton(P, x, steps=5, bits=256):308 D = [i * P[i] for i in range(1, len(P))]309 v = Fraction(x)310 for _ in range(steps):311 d = poly_at(D, v)312 if not d:313 return v314 v -= Fraction(poly_at(P, v)) / d315 v = Fraction(round(v * 2 ** bits), 2 ** bits)316 return v317318def poly_isolate(P, x, digits):319 w = Fraction(1, 10 ** 6)320 cap = Fraction(abs(x)) / 1000 + 1321 while w < cap and poly_at(P, Fraction(x) - w) * poly_at(P, Fraction(x) + w) > 0:322 w *= 2323 if w >= cap:324 return None325 lo, hi = Fraction(x) - w, Fraction(x) + w326 eps = Fraction(1, 10 ** (digits + 3))327 while hi - lo > eps:328 mid = (lo + hi) / 2329 if poly_at(P, lo) * poly_at(P, mid) <= 0:330 hi = mid331 else:332 lo = mid333 return lo, hi334335def poly_radical(P):336 if len(P) == 2:337 return str(-P[0])338 if len(P) != 3:339 return "-"340 p, q = -P[1], P[0]341 d = p * p - 4 * q342 f = 1343 t = 2344 while t * t <= d:345 while d % (t * t) == 0:346 d //= t * t347 f *= t348 t += 1349 g = gcd(gcd(p, f), 2)350 num = str(p // g) + " + " + (str(f // g) + " " if f // g != 1 else "") + "sqrt " + str(d)351 return num if 2 // g == 1 else "(" + num + ") / 2"352353def linear_solve(rows, rhs):354 n = len(rows[0])355 A = [[Fraction(v) for v in r] + [Fraction(rhs[i])] for i, r in enumerate(rows)]356 piv, r = [], 0357 for c in range(n):358 p = next((i for i in range(r, len(A)) if A[i][c]), None)359 if p is None:360 continue361 A[r], A[p] = A[p], A[r]362 A[r] = [v / A[r][c] for v in A[r]]363 for i in range(len(A)):364 if i != r and A[i][c]:365 f = A[i][c]366 A[i] = [A[i][t] - f * A[r][t] for t in range(n + 1)]367 piv.append(c)368 r += 1369 if r == len(A):370 break371 for i in range(r, len(A)):372 if A[i][n] and not any(A[i][:n]):373 return None374 x = [Fraction(0)] * n375 for i, c in enumerate(piv):376 x[c] = A[i][n]377 return x378379def min_recurrence(seq, maxord):380 for L in range(1, maxord + 1):381 if len(seq) < 2 * L + 6:382 return None383 c = linear_solve([[seq[i + t] for t in range(L)] for i in range(L + 2)],384 [seq[i + L] for i in range(L + 2)])385 if c is None:386 continue387 if all(sum(c[t] * seq[i + t] for t in range(L)) == seq[i + L] for i in range(len(seq) - L)):388 return [-int(v) for v in c] + [1]389 return None390391def lam_block(b, terms=0, maxord=0, kbig=160):392 terms = terms or 4 * b + 14393 maxord = maxord or 2 * b + 4394 seq = [block_count(b, k) for k in range(3, 3 + terms)]395 P = min_recurrence(seq, maxord)396 if P is None:397 return None398 rs = poly_roots(P, 2.0 ** b / 3 + 1)399 lam = max(rs, key=lambda z: z.real)400 ratio = block_count(b, kbig + 1) / block_count(b, kbig)401 if abs(lam - ratio) > 1e-2 * ratio:402 raise ValueError("the dominant root and the block ratio disagree at b = " + str(b))403 Q = poly_minfactor(P, rs, lam)404 lam2 = max(abs(z) for z in rs if z is not lam) if len(rs) > 1 else 0.0405 return P, Q, poly_isolate(Q or P, poly_newton(Q or P, lam.real), 12), ratio, lam2 / lam.real406407# THE LIFT COUNT408409def lift_support(word):410 k = len(word) + 1411 T = [i + 1 for i, c in enumerate(word) if c]412 m = 1 + 2 * sum(3 ** i for i in T)413 supp = [p for p in range(k) if p not in T] + [k + p for p in T]414 return k, m, sorted(supp)415416def submask_count(word):417 k, m, supp = lift_support(word)418 h = len(supp) // 2419 lo, hi = supp[:h], supp[h:]420 d = {}421 for s in range(1 << len(lo)):422 v = sum(pow(3, lo[i], m) for i in range(len(lo)) if s >> i & 1) % m423 d[v] = d.get(v, 0) + 1424 out = 0425 for s in range(1 << len(hi)):426 v = sum(pow(3, hi[i], m) for i in range(len(hi)) if s >> i & 1) % m427 out += d.get((-v) % m, 0)428 return out429430def submask_set(word):431 k, m, supp = lift_support(word)432 n = len(supp)433 return [s for s in range(1 << n) if sum(3 ** supp[i] for i in range(n) if s >> i & 1) % m == 0], n434435def closure_faults(D, n):436 S = set(D)437 full = (1 << n) - 1438 bad = 0439 for a in D:440 if (full ^ a) not in S:441 bad += 1442 for b in D:443 if not a & b and (a | b) not in S:444 bad += 1445 if a & b == b and (a ^ b) not in S:446 bad += 1447 return bad448449def irreducibles(D):450 S = set(D)451 out = []452 for a in D:453 if not a:454 continue455 if not any(b and b != a and b & a == b and (a ^ b) in S for b in D):456 out.append(a)457 return out458459def lift_census(k):460 tot = triv = pw = dis = 0461 ito = imax = 0462 cls = {}463 for w in range(1 << (k - 1)):464 word = tuple((w >> i) & 1 for i in range(k - 1))465 D, n = submask_set(word)466 I = irreducibles(D)467 tot += len(D)468 ito += len(I)469 imax = max(imax, len(I))470 triv += len(D) == 2471 pw += not len(D) & (len(D) - 1)472 dis += all(not I[a] & I[b] for a in range(len(I)) for b in range(a))473 r = lift_support(word)[1] % 9474 c = cls.setdefault(r, [0, 0, 0])475 c[0] += 1476 c[1] += len(D)477 c[2] += len(I)478 return tot, ito, imax, triv, pw, dis, cls479480def rank_exact(rows):481 rows = [[Fraction(x) for x in r] for r in rows]482 n = len(rows[0])483 r = 0484 for c in range(n):485 p = next((i for i in range(r, len(rows)) if rows[i][c]), None)486 if p is None:487 continue488 rows[r], rows[p] = rows[p], rows[r]489 pv = rows[r][c]490 for i in range(len(rows)):491 if i != r and rows[i][c]:492 f = rows[i][c] / pv493 for j in range(c, n):494 rows[i][j] -= f * rows[r][j]495 r += 1496 if r == len(rows):497 break498 return r499500def all_words(n):501 return [tuple((s >> i) & 1 for i in range(L)) for L in range(n + 1) for s in range(1 << L)]502503def hankel_rank(p, q, fn=None):504 fn = fn or submask_count505 U, V = all_words(p), all_words(q)506 memo = {}507 H = []508 for u in U:509 row = []510 for v in V:511 w = u + v512 if w not in memo:513 memo[w] = fn(w)514 row.append(memo[w])515 H.append(row)516 return rank_exact(H), len(U)517518def irreducible_count(word):519 D, _ = submask_set(word)520 return len(irreducibles(D))521522# THE FLOOR AND THE MODEL523524def floor_count(k):525 w = (3 ** k - 1) // 2526 pw = [3 ** i for i in range(k)]527 out = 0528 for s in range(1, (1 << k) - 1):529 v = sum(pw[i] for i in range(k) if s >> i & 1)530 if gcd(v, w) == 1:531 out += 1532 return out533534def model_by_top(k):535 out = [0.0] * k536 for t in range(1, k):537 base = 1 + 2 * 3 ** t538 acc = 0.0539 for s in range(1 << (t - 1)):540 a = sum(3 ** (i + 1) for i in range(t - 1) if s >> i & 1)541 acc += 1.0 / (base + 2 * a)542 out[t] = acc * 2 ** k543 return out544545# THE FITS546547TQ = (12.706, 4.303, 3.182, 2.776, 2.571, 2.447, 2.365, 2.306, 2.262, 2.228,548 2.201, 2.179, 2.16, 2.145, 2.131, 2.12, 2.11, 2.101, 2.093, 2.086,549 2.08, 2.074, 2.069, 2.064, 2.06, 2.056, 2.052, 2.048, 2.045, 2.042)550551def ols(xs, ys):552 n = len(xs)553 mx, my = sum(xs) / n, sum(ys) / n554 sxx = sum((x - mx) ** 2 for x in xs)555 sl = sum((xs[i] - mx) * (ys[i] - my) for i in range(n)) / sxx556 a = my - sl * mx557 se = (sum((ys[i] - a - sl * xs[i]) ** 2 for i in range(n)) / (n - 2) / sxx) ** 0.5558 t = TQ[n - 3] if 3 <= n <= 32 else 1.96559 return sl, se, sl - t * se, sl + t * se, n560561def ols_diag(xs, ys):562 n = len(xs)563 mx = sum(xs) / n564 sxx = sum((x - mx) ** 2 for x in xs)565 sl = sum((xs[i] - mx) * (ys[i] - sum(ys) / n) for i in range(n)) / sxx566 a = sum(ys) / n - sl * mx567 r = [ys[i] - a - sl * xs[i] for i in range(n)]568 s2 = sum(v * v for v in r) / (n - 2)569 h = [1 / n + (xs[i] - mx) ** 2 / sxx for i in range(n)]570 dw = sum((r[i] - r[i - 1]) ** 2 for i in range(1, n)) / sum(v * v for v in r)571 m = n // 2572 v1 = sum(v * v for v in r[:m]) / m573 v2 = sum(v * v for v in r[m:]) / (n - m)574 stu = max(abs(r[i]) / (s2 * (1 - h[i])) ** 0.5 for i in range(n))575 cook = max(r[i] ** 2 * h[i] / (2 * s2 * (1 - h[i]) ** 2) for i in range(n))576 return dw, max(v1, v2) / min(v1, v2), stu, cook577578def ols2(x1, x2, ys):579 n = len(ys)580 d = [sum(x1) / n, sum(x2) / n, sum(ys) / n]581 u = [v - d[0] for v in x1]582 w = [v - d[1] for v in x2]583 y = [v - d[2] for v in ys]584 a, b, c = sum(v * v for v in u), sum(u[i] * w[i] for i in range(n)), sum(v * v for v in w)585 e, f = sum(u[i] * y[i] for i in range(n)), sum(w[i] * y[i] for i in range(n))586 det = a * c - b * b587 return (c * e - b * f) / det, (a * f - b * e) / det588589def survival(hist, k):590 s, b = [], 2591 while True:592 v = sum(x for d, x in hist.items() if d > b * k)593 if not v:594 break595 s.append(v)596 b += 1597 return s598599def geom_fit(s):600 n = [s[i] - s[i + 1] for i in range(len(s) - 1)] + [s[-1]]601 tot = sum(n)602 r = 1 - tot / sum((i + 1) * n[i] for i in range(len(n)))603 e = [tot * (1 - r) * r ** i for i in range(len(n))]604 raw = sum((n[i] - e[i]) ** 2 / e[i] for i in range(len(n)))605 o, x, co, ce = [], [], 0, 0.0606 for i in range(len(n)):607 co += n[i]608 ce += e[i]609 if ce >= 5:610 o.append(co)611 x.append(ce)612 co, ce = 0, 0.0613 if co:614 o[-1] += co615 x[-1] += ce616 return r, raw, len(n) - 2, sum((o[i] - x[i]) ** 2 / x[i] for i in range(len(o))), len(o) - 2617618def dec_down(x, d):619 n = x.numerator * 10 ** d // x.denominator620 return str(n // 10 ** d) + "." + str(n % 10 ** d).rjust(d, "0")621622def dec_up(x, d):623 n = -((-x.numerator * 10 ** d) // x.denominator)624 return str(n // 10 ** d) + "." + str(n % 10 ** d).rjust(d, "0")625626# SWEEP627628def sweep(k):629 w = (3 ** k - 1) // 2630 hist = {}631 cand = 0632 for z1 in range(3, w, 3):633 if gcd(z1, w) != 1:634 continue635 cand += 1636 d = first_return(z1, w - z1)637 if d:638 hist[d] = hist.get(d, 0) + 2639 return hist, 2 * cand640641# VERBS642643def cmd_automaton(args):644 print("the band automaton of the direction (z1, z2), z1 + z2 = w, 3 | z1: states are the integers j in "645 "[-(z2-1)//2, (z1-1)//2]; from j the moves are j -> (j + a)/3 over the increments a in {0, z1, -z2} "646 "that keep the quotient integral, so out-degree 2 on j = 0 mod 3, 1 on j = z2 mod 3, 0 on the third class; "647 "the walk leaves 0 by the forced increment z1 and a return is a walk back to 0, of length n exactly when "648 "the multiplier m it spells has m w binary of base-3 length n")649 print("w z1 z2 states edges mean deg0 deg1 deg2 firstreturn")650 for w in args.weights:651 for z1 in range(3, w, 3):652 z2 = w - z1653 if gcd(z1, z2) != 1:654 continue655 cnt, edges, states = degree_profile(z1, z2)656 if z1 % 30 and w > 40:657 continue658 print(w, z1, z2, states, edges, down(edges / states, 6), cnt[0], cnt[1], cnt[2],659 first_return(z1, z2))660 print("criticality over every coprime direction of these weights; the mean out-degree is exactly "661 "1 + (n0 - n2) / N for the counts n0, n1, n2 of the residue classes 0, z2 and the third among the N "662 "band states, so |mean - 1| <= 1 / N, and mean = 1 exactly when n0 = n2, which 3 | N gives and does "663 "not exhaust")664 every = flat = split = 0665 for w in args.weights:666 worst = 0.0667 tot = 0668 for z1 in range(3, w, 3):669 z2 = w - z1670 if gcd(z1, z2) != 1:671 continue672 lo, hi = band(z1, z2)673 n0 = sum(1 for j in range(lo, hi + 1) if not j % 3)674 n2 = sum(1 for j in range(lo, hi + 1) if (j - z2) % 3 and j % 3)675 cnt, edges, states = degree_profile(z1, z2)676 if 2 * cnt[2] + cnt[1] != edges or edges != states + n0 - n2 or abs(n0 - n2) > 1:677 raise ValueError("the out-degree is not carried by the residue classes at " + str((w, z1)))678 worst = max(worst, abs(edges / states - 1))679 tot += 1680 every += 1681 if n0 == n2:682 flat += 1683 split += states % 3 != 0684 print("w", w, "directions", tot, "max |mean out-degree - 1|", up(worst, 6))685 print("over", every, "coprime directions of these weights the mean out-degree is exactly 1 on", flat,686 "of them, and", split, "of those have 3 not dividing the band size, so 3 | N is sufficient for exact "687 "criticality and not necessary")688 print("the two smallest witnesses, read off the moves themselves:")689 for z1, z2 in ((3, 1), (3, 2)):690 cnt, edges, states = degree_profile(z1, z2)691 lo, hi = band(z1, z2)692 print(" (z1, z2) =", (z1, z2), "states", list(range(lo, hi + 1)), "degrees",693 [len(succ(j, z1, z2)) for j in range(lo, hi + 1)], "N", states, "mean out-degree",694 down(edges / states, 6), "| 3 divides N:", not states % 3)695 print("so (3, 1) is the smallest non-critical direction, mean 3/2, and (3, 2) is critical at N = 2 with 3 "696 "not dividing N: criticality is a statement about the residue split of the band and never an identity")697698def cmd_returns(args):699 print("L(k, n) = #{m >= 1 : m R_k binary in base 3 and below 3^n} is the return count of the weight R_k at "700 "horizon n, the number of distinct returns the band automata of weight R_k can spell; it is computed by "701 "the column transfer, never by enumeration; free = 2^n / R_k is the count a random binary string of "702 "length n would give and rho = L R_k / 2^n the excess over it")703 print("depth b, L(k, bk) at k = " + str(args.kmin) + "..; lam = L(k+1, b(k+1)) / L(k, bk) at two large k, "704 "the growth of the return count per block, against the free rate 2^b / 3; these ratios are a reading "705 "and not the rate, the second root of the recurrence sitting within a percent of lam_b at even b, so "706 "the verb ladder is where lam_b is computed exactly")707 for b in range(1, args.bmax + 1):708 row = [lift_count(k, b * k) for k in range(args.kmin, args.kmin + 6)]709 r1 = lift_count(args.kbig + 1, b * (args.kbig + 1)) / lift_count(args.kbig, b * args.kbig)710 r2 = lift_count(2 * args.kbig + 1, b * (2 * args.kbig + 1)) / lift_count(2 * args.kbig, b * 2 * args.kbig)711 free = 2.0 ** b / 3.0712 print("b", b, row, "lam", down(r1, 6), down(r2, 6), "free", down(free, 6),713 "excess", down(r2 * 3 / 2.0 ** b, 6))714 print("the raw counts, every m including 3 | m:")715 for b in range(1, args.bmax + 1):716 row = [lift_count(k, b * k, False) for k in range(args.kmin, args.kmin + 6)]717 print("b", b, row)718 print("rho = L(k, bk) R_k / 2^(bk - 1), the excess over the free model:")719 for b in range(2, args.bmax + 1):720 row = []721 for k in range(args.kmin, args.kmax + 1):722 w = (3 ** k - 1) // 2723 row.append(down(lift_count(k, b * k) * w / 2.0 ** (b * k - 1), 4))724 print("b", b, row)725 k = args.kfine726 print("the fine return count at k =", k, ": n, L(k, n), new = L(k, n) - L(k, n-1)")727 prev = 0728 gaps = []729 for n in range(k, args.bmax * k + 1):730 cur = lift_count(k, n)731 if cur == prev:732 gaps.append(n)733 prev = cur734 print("n with no return at all, k <", k, "* bmax:", gaps)735 for b in range(1, args.bmax + 1):736 print("b", b, "L", lift_count(k, b * k), "L at bk+1", lift_count(k, b * k + 1))737738def cmd_hist(args):739 print("d(z) is the first return time of the band automaton of (z, R_k - z), the base-3 length of the shortest "740 "binary lift m R_k that carries z; Z = #{z : d(z) < infinity}, Phi = #{d = k} the coprime submask floor, "741 "U = #{d <= 2k}, V = #{d <= 3k}, tail = Z - U")742 print("k R_k cand Z Phi U V tail dmax dmax/sqrt(R_k) secs")743 rows = []744 for k in range(args.kmin, args.kmax + 1):745 t = time.time()746 hist, cand = sweep(k)747 Z = sum(hist.values())748 phi = sum(v for d, v in hist.items() if d <= k)749 U = sum(v for d, v in hist.items() if d <= 2 * k)750 V = sum(v for d, v in hist.items() if d <= 3 * k)751 mx = max(hist) if hist else 0752 w = (3 ** k - 1) // 2753 print(k, w, cand, Z, phi, U, V, Z - U, mx, down(mx / w ** 0.5, 4), round(time.time() - t, 1))754 rows.append((k, hist, Z, phi, U, V))755 print("the fine histogram in the bulk, h(t) = #{z : d(z) = k + t + 1} against the model "756 "E(t) = Sum_{max T = t} 2^k / m_T, the model computed from the multipliers alone:")757 for k, hist, Z, phi, U, V in rows:758 if k < args.bulkmin:759 continue760 mod = model_by_top(k)761 h = [hist.get(k + t + 1, 0) for t in range(k)]762 print("k", k, "h", h[1:], "model", [down(mod[t], 1) for t in range(1, k)],763 "ratio", [down(h[t] / mod[t], 4) if mod[t] else "-" for t in range(1, k)])764 print("the depth histogram of the tail, b -> #{z : d(z) in ((b-1)k, bk]}:")765 for k, hist, Z, phi, U, V in rows:766 if Z == U:767 continue768 dh = {}769 for d, v in hist.items():770 if d > 2 * k:771 dh[-(-d // k)] = dh.get(-(-d // k), 0) + v772 print("k", k, sorted(dh.items()))773 print("the support of the first return time inside [k, dmax]: the lengths reached, the lengths missing, and "774 "the first holes past the proved gap at k + 1:")775 for k, hist, Z, phi, U, V in rows:776 ds = sorted(hist)777 miss = [n for n in range(k, ds[-1] + 1) if n not in hist]778 print("k", k, "dmax", ds[-1], "distinct lengths", len(ds), "missing", len(miss),779 "first four holes", miss[:4], "first six lengths", ds[:6])780 print("the tail survival S(b) = #{z : d(z) > bk}, in z values and in distinct directions, a direction being "781 "the pair (z, R_k - z) which the sweep counts twice, with the sample size behind every exponent and a "782 "maximum-likelihood geometric fitted to the depth counts:")783 for k, hist, Z, phi, U, V in rows:784 if Z == U:785 continue786 s = survival(hist, k)787 if any(v % 2 for v in s):788 raise ValueError("a survival count is odd at k = " + str(k))789 print("k", k, "S(2..) z values", s)790 print(" S(2..) directions", [v // 2 for v in s], "ratios",791 [down(s[i + 1] / s[i], 4) for i in range(len(s) - 1)])792 e = []793 b = 2794 while 2 * b < len(s) + 2:795 e.append((b, down(-log(s[2 * b - 2] / s[b - 2]) / log(2), 3) if s[2 * b - 2] else "-",796 s[b - 2] // 2, s[2 * b - 2] // 2))797 b *= 2798 print(" local exponent -log2(S(2b)/S(b)) at (b, exponent, directions behind S(b), behind S(2b))", e)799 r, raw, rdf, pool, pdf = geom_fit(s)800 print(" one geometric fitted by maximum likelihood to the depth counts: ratio", down(r, 4), ", chi2",801 up(raw, 1), "on", rdf, "df over the unpooled depths and", up(pool, 1) if pdf > 0 else "-", "on", pdf,802 "df with the bins pooled to expectation 5; the unpooled depths fall below expectation 5 in the "803 "tail, where the statistic is not valid, so neither verdict is carried by the counts")804 print("Phi against the independent floor count, and the depth-1 identity:")805 for k, hist, Z, phi, U, V in rows:806 f = floor_count(k)807 if f != phi:808 raise ValueError("the floor " + str(f) + " against the depth-1 count " + str(phi) + " at k = " + str(k))809 print("k", k, "Phi", f, "matched")810811def brute_depths(k, n):812 w = (3 ** k - 1) // 2813 out = {}814 for c in range(1, 1 << n):815 K = sum(3 ** i for i in range(n) if c >> i & 1)816 if K % w:817 continue818 m = K // w819 d = max(i for i in range(n) if c >> i & 1) + 1820 sup = [i for i in range(n) if c >> i & 1]821 for t in range(1, 1 << len(sup)):822 A = sum(3 ** sup[i] for i in range(len(sup)) if t >> i & 1)823 if A % m:824 continue825 z = A // m826 if 0 < z < w and gcd(z, w) == 1 and (out.get(z, n + 1) > d):827 out[z] = d828 return out829830def sweep_weight(w):831 hist = {}832 for z1 in range(3, w, 3):833 if gcd(z1, w) != 1:834 continue835 d = first_return(z1, w - z1)836 if d:837 hist[d] = hist.get(d, 0) + 1838 return hist839840def cmd_critical(args):841 print("the deepest first return of the band automaton over every coprime direction of a weight w prime to 3, "842 "against the scale sqrt(w) a critical walk on a band of w/2 states predicts; Z is the occupied count, "843 "med the median first return and q9 its ninth decile")844 print("w Z dmax med q9 dmax/sqrt(w) med/sqrt(w) q9/sqrt(w) secs")845 rows = []846 w = args.wmin847 while w <= args.wmax:848 if w % 3:849 t = time.time()850 hist = sweep_weight(w)851 if hist:852 Z = sum(hist.values())853 ds = []854 for d in sorted(hist):855 ds += [d] * hist[d]856 med, q9 = ds[len(ds) // 2], ds[9 * len(ds) // 10]857 rows.append((w, Z, max(hist), med, q9))858 print(w, Z, max(hist), med, q9, down(max(hist) / w ** 0.5, 4), down(med / w ** 0.5, 4),859 down(q9 / w ** 0.5, 4), round(time.time() - t, 1))860 w = int(w * args.step) + 1861 band = [r[2] / r[0] ** 0.5 for r in rows]862 print("dmax / sqrt(w) over", len(rows), "weights: min", down(min(band), 4), "max", up(max(band), 4),863 "mean", down(sum(band) / len(band), 4))864 print("the exponent the ladder carries: log d on log w by ordinary least squares, with the 95 percent "865 "interval and the t statistic against the critical exponent 1/2:")866 for name, sel, col in (("dmax, every weight", lambda r: True, 2), ("dmax, Z >= 8", lambda r: r[1] >= 8, 2),867 ("median, every weight", lambda r: True, 3), ("median, Z >= 4", lambda r: r[1] >= 4, 3),868 ("median, Z >= 8", lambda r: r[1] >= 8, 3), ("median, Z >= 16", lambda r: r[1] >= 16, 3),869 ("ninth decile", lambda r: True, 4)):870 sub = [r for r in rows if sel(r)]871 if len(sub) < 4:872 continue873 sl, se, lo, hi, n = ols([log(r[0]) for r in sub], [log(r[col]) for r in sub])874 print(" ", name, "n", n, "exponent", down(sl, 4), "95 percent [", down(lo, 4), ",", up(hi, 4),875 "] t against 1/2", down((sl - 0.5) / se, 2), "| excludes 1/2:", not lo <= 0.5 <= hi)876 print("the median is not one number: a weight with Z <= 2 has its median equal to its dmax, and", 877 sum(1 for r in rows if r[1] <= 2), "of the", len(rows), "weights have Z <= 2 and",878 sum(1 for r in rows if r[1] <= 4), "have Z <= 4, so the median exponent moves with the cut and only "879 "its sign below 1/2 survives every cut")880 loo = []881 for i in range(len(rows)):882 sub = rows[:i] + rows[i + 1:]883 sl, se, lo, hi, n = ols([log(r[0]) for r in sub], [log(r[2]) for r in sub])884 loo.append((sl, lo, hi, rows[i][0]))885 a, z = min(loo), max(loo)886 print("leave one out on the dmax fit: the slope ranges over [", down(a[0], 4), ",", up(z[0], 4),887 "], the floor deleting w =", a[3], "and the ceiling deleting w =", z[3], ";",888 sum(1 for v in loo if v[1] <= 0.5 <= v[2]), "of the", len(loo),889 "leave-one-out intervals cover 1/2 at the weights", " ".join(str(v[3]) for v in loo if v[1] <= 0.5 <= v[2]),890 ", the widest upper endpoint", up(max(v[2] for v in loo), 4),891 "- so the exclusion of 1/2 rests on single weights and is not a property of the ladder; what survives "892 "every deletion is the sign, dmax ~ w^0.41 with the exponent below 1/2")893 c1, c2 = ols2([log(r[0]) for r in rows], [log(r[1]) for r in rows], [log(r[2]) for r in rows])894 print("two predictors, log dmax on log w and log Z:", down(c1, 4), "on log w and", down(c2, 4),895 "on log Z, both positive, so controlling for the sample size lowers the exponent on the weight: the "896 "one-predictor slope overstates and the drift below 1/2 is understated, never overstated")897 dw, vr, stu, cook = ols_diag([log(r[0]) for r in rows], [log(r[2]) for r in rows])898 print("the dmax fit's diagnostics: Durbin-Watson", down(dw, 2), ", residual variance ratio across the "899 "halves of the ladder", up(vr, 2), ", largest studentised residual", up(stu, 2),900 ", largest Cook distance", up(cook, 3), "- ordinary least squares is not what breaks here; the "901 "leverage of single weights and the Z <= 2 weights are")902 sl, se, lo, hi, n = ols([log(r[0]) for r in rows], [log(r[1]) for r in rows])903 print(" log Z on log w: exponent", down(sl, 4), "95 percent [", down(lo, 4), ",", up(hi, 4),904 "], so the sample a weight offers grows with the weight")905906def cmd_model(args):907 print("the model return count D(k, N) = Sum over the primitive lifts K of length at most N of 2^L(K) / m(K), "908 "the equidistribution count of (lift, direction) pairs the band automata of weight R_k spell inside "909 "horizon N; 2^L(K) is summed by the same column transfer and 1/m(K) is sandwiched by the length of K, "910 "so lo and hi bracket D with lo the safe lower and hi the safe upper; D(k, 2k) must reproduce 2^k L_k "911 "with L_k = Sum_T 1 / m_T, the aggregate the lift half already carries")912 print("k L_k 2^k L_k D(k,2k) lo hi ratio")913 for k in range(args.kmin, args.kmax + 1):914 rows = model_returns(k, 2 * k)915 lo, hi = rows[-1][1], rows[-1][2]916 Lk = exact_L(k)917 print(k, down(Lk, 5), down(2 ** k * Lk, 2), "[", down(lo, 2), up(hi, 2), "]",918 down(lo / (2 ** k * Lk), 4), up(hi / (2 ** k * Lk), 4))919 print("D(k, bk) / 2^k at every depth b, the aggregate return count the model gives with no enumeration:")920 for k in range(args.kmin, args.kmax + 1):921 rows = model_returns(k, args.bmax * k)922 band = []923 for b in range(2, args.bmax + 1):924 n, lo, hi = rows[b * k - k]925 band.append((b, down(lo / 2 ** k, 3), up(hi / 2 ** k, 3)))926 print("k", k, band)927 print("the deep part at the critical cutoff N = floor(sqrt(R_k)), the scale a critical walk on the band "928 "predicts for the deepest return; beyond a depth of 40 blocks the per-digit increment of D is exactly "929 "4/9 and the row extends by it")930 print("share is the percentage of the deep part carried by the 4/9 extrapolation past 40 blocks rather than "931 "by the transfer; both endpoints of the extrapolated stretch are built on the free per-digit increment "932 "4/9 and its triple 4/3, and rho > 1 at every depth puts the true increment above 4/9, so the "933 "extrapolated part of the band leans low at both ends and hi stays an upper bound only while rho < 3")934 print("k N D(k,2k)/2^k D(k,N)/2^k deep = (D(k,N) - D(k,2k))/2^k share")935 for k in range(args.kmin, args.kmax + 1):936 w = (3 ** k - 1) // 2937 N = int(w ** 0.5)938 rows = model_returns(k, N)939 base = rows[k][1], rows[k][2]940 cap = rows[min(N, 40 * k) - k][1]941 top = rows[-1][1], rows[-1][2]942 deep = top[0] - base[0]943 print(k, N, "[", down(base[0] / 2 ** k, 4), up(base[1] / 2 ** k, 4), "]",944 "[", down(top[0] / 2 ** k, 4), up(top[1] / 2 ** k, 4), "]",945 "[", down((top[0] - base[0]) / 2 ** k, 4), up((top[1] - base[1]) / 2 ** k, 4), "]",946 up(100 * (top[0] - cap) / deep, 0) if deep else 0)947948def cmd_ladder(args):949 print("at n = bk the slot profile is uniform, s_r = b for every column r, so the column transfer of Q.2 is a "950 "matrix fixed in k: L(k, bk) = Sum_{j = 0..b} w_j B(b, j)^(k - r0(j)) h_j with B(b, j), h_j, w_j and "951 "the head length r0(j) all independent of k, hence L(k, bk) obeys a constant-coefficient linear "952 "recurrence in k whose dominant root is lam_b; the carry state space is [0, b // 2], every column sum "953 "of every B(b, j) is Sum_{c = a mod 3} binom(b, c) = (2^b + 2 cos(pi (b - 2a) / 3)) / 3, so the "954 "deviation from 2^b / 3 takes exactly two values and they are set by the parity of b, and a "955 "nonnegative matrix has its spectral radius between its least and its greatest column sum")956 print("the parity bracket, asserted per b below: 2^b / 3 - 1/3 <= lam_b <= 2^b / 3 + 2/3 at even b and "957 "2^b / 3 - 2/3 <= lam_b <= 2^b / 3 + 1/3 at odd b, so 3 lam_b - 2^b lies in [-1, 2] at even b and in "958 "[-2, 1] at odd b and the headline is the two-sided |3 lam_b / 2^b - 1| <= 2^(1 - b) at every b; the "959 "parity refines which edge is which, not the two-sided rate")960 print("the fixed matrix is asserted against the general transfer at every k = 3.." + str(args.kcheck))961 exc = []962 for b in range(args.bmin, args.bmax + 1):963 for k in range(3, args.kcheck + 1):964 if set(slot_profile(k, b * k)) != {b}:965 raise ValueError("the slot profile is not uniform at " + str((k, b)))966 if block_count(b, k) != lift_count(k, b * k):967 raise ValueError("the fixed matrix misses the transfer at " + str((k, b)))968 dev, clo, chi, states, r0s, rmax = block_bracket(b)969 dlo, dhi = dev[0], dev[-1]970 exc.append((b, None))971 out = lam_block(b)972 if out is None:973 print("b", b, "no recurrence of order at most", 2 * b + 4)974 continue975 P, Q, iv, ratio, sub = out976 if iv is None:977 print("b", b, "order", len(P) - 1, "lam_b is not isolated by the exact bisection")978 continue979 e0, e1 = 3 * iv[0] - 2 ** b, 3 * iv[1] - 2 ** b980 sl = Fraction(1, 10 ** 9)981 if e0 < 3 * dlo - sl or e1 > 3 * dhi + sl:982 raise ValueError("lam_b leaves the column-sum bracket at b = " + str(b))983 if abs(Fraction(ratio) - iv[0]) > iv[1] / 100:984 raise ValueError("the isolated root misses L(k+1, b(k+1)) / L(k, bk) at b = " + str(b))985 print("b", b, "order", len(P) - 1, "charpoly", poly_show(P))986 if Q is None:987 print(" the minimal polynomial of lam_b is not identified; lam_b is a root of the polynomial above")988 else:989 print(" deg", len(Q) - 1, "minpoly of lam_b", poly_show(Q), "| exact",990 poly_radical(Q) if len(Q) <= 3 else991 "in radicals, degree " + str(len(Q) - 1) if len(Q) <= 5 else992 "no radical form, degree " + str(len(Q) - 1))993 print(" lam_b in [", dec_down(iv[0], 9), ",", dec_up(iv[1], 9), "] 3 lam_b - 2^b in [",994 dec_down(e0, 6), ",", dec_up(e1, 6), "] column sum - 2^b / 3 in {",995 ", ".join(str(v) for v in dev), "} saturated:", e0 <= 3 * dlo + sl or e1 >= 3 * dhi - sl)996 print(" head lengths r0(j), j = 0..b:", " ".join(str(v) for v in r0s), " so the exponent k - r0(j) is "997 "k - 2 for", r0s.count(2), "of the", b + 1, "sectors")998 print(" the j = 1 closure carries", states, "of", b // 2 + 1, "carry states and its column sums lie in "999 "2^b / 3 + [", clo, ",", chi, "]; max_j rho(B(b, j)) reads", down(rmax, 6), "against lam_b, gap",1000 up(abs(rmax - float(iv[0])) / float(iv[0]), 9))1001 if abs(rmax - float(iv[0])) > 1e-6 * float(iv[0]):1002 raise ValueError("a block outgrows lam_b at b = " + str(b))1003 exc[-1] = (b, e0)1004 print(" the block ratio L(k + 1, b(k + 1)) / L(k, bk) at k = 160 reads", down(ratio, 6),1005 "and the second root of the recurrence is", up(sub, 6), "of lam_b, so that ratio holds about",1006 int(-160 * log(sub, 10)) if 0 < sub < 1 else 0,1007 "correct digits at k = 160: a ratio is not how this rate is read")1008 good = [v for v in exc if v[1] is not None]1009 print("the excess 3 lam_b - 2^b, truncated down, at b =", good[0][0], "..", good[-1][0], ":",1010 ", ".join(dec_down(v[1], 6) for v in good))1011 print("every one of them is positive, so lam_b sits above the free rate 2^b / 3 at every depth printed and "1012 "the live edge of the parity bracket is the upper one, while the two-sided rate stays 2^(1 - b)")10131014def cmd_lift(args):1015 t0 = time.time()1016 print("the depth-2 lift census: K = m R_k with one position per column, N_K(m) the submasks of K "1017 "divisible by m, M_k = Sum_T N_K(m), the sweep exhaustive over all 2^(k-1) sets T inside [1, k-1]")1018 seq = []1019 band = []1020 for k in range(1, args.kbig + 1):1021 M = sum(submask_count(tuple((w >> i) & 1 for i in range(k - 1))) for w in range(1 << (k - 1)))1022 seq.append(M)1023 band.append(M / 2.0 ** k)1024 print("k", k, "M_k", M, "M_k / 2^k", repr(M / 2.0 ** k))1025 print("M_k / 2^k lies inside [" + str(down(min(band), 5)) + ", " + str(up(max(band), 5)) + "], the endpoints "1026 "rounded outward so the band holds the exact values")1027 n = (len(seq) + 1) // 21028 H = [[seq[i + j] for j in range(n)] for i in range(n)]1029 print("the Hankel matrix of M_k is", n, "by", n, "of rank", rank_exact(H), "reading every one of the",1030 len(seq), "terms k = 1.." + str(len(seq)) + ", so no linear recurrence of order at most", n - 1,1031 "holds on them")1032 print("")1033 print("k M_k/2^k Sum_T iota_T / 2^k max_T iota_T trivial T power-of-two N disjoint irreducibles m_T mod 9")1034 ib = []1035 for k in range(2, args.kmax + 1):1036 tot, ito, imax, triv, pw, dis, cls = lift_census(k)1037 pk = 2.0 ** k1038 ib.append(ito / pk)1039 print("k", k, repr(tot / pk), repr(ito / pk), imax, str(triv) + "/" + str(1 << (k - 1)),1040 str(pw) + "/" + str(1 << (k - 1)), str(dis) + "/" + str(1 << (k - 1)),1041 " ".join(str(r) + ":" + str(c[0]) + ":" + str(round(c[1] / pk, 4)) + ":" + str(round(c[2] / pk, 4))1042 for r, c in sorted(cls.items())))1043 print("Sum_T iota_T / 2^k lies inside [" + str(down(min(ib), 6)) + ", " + str(up(max(ib), 6)) + "], the "1044 "endpoints rounded outward; the even readings fall at every step from k = 6 and the odd ones do not")1045 print("the disjoint column counts the T whose irreducibles are pairwise disjoint, where N_K(m) = 2^iota; "1046 "the power-of-two column is larger, so a power-of-two count does not force disjointness")1047 print("the residue classes of the ladder are m_T mod 9, printed as class:sets:M share:iota share; T inside "1048 "[1, k-1] gives a_T = Sum 3^i with i >= 1, so a_T is 3 mod 9 when 1 is in T and 0 mod 9 otherwise, and "1049 "m_T = 1 + 2 a_T is 7 or 1 mod 9 and never 4, class 4 needing a_T = 6 mod 9")1050 print("")1051 top = 01052 for n in range(1, args.nmax + 1):1053 a, rows = hankel_rank(n, n)1054 top = max(top, a)1055 print("Hankel of N over the column word, p = q =", n, "rows", rows, "rank", a, "full rank", rows)1056 print("the reversed reading is the transpose of this matrix at p = q and carries no further information")1057 for n in range(1, min(args.nmax, 5) + 1):1058 a, rows = hankel_rank(n, n, irreducible_count)1059 print("Hankel of the irreducible count iota over the column word, p = q =", n, "rows", rows,1060 "rank", a, "full rank", rows)1061 print("a column transfer with a state set free of k is a linear representation of N as a series over the "1062 "column word, so its dimension is at most its state count and at least the Hankel rank; the rank "1063 "reaches", top, "on the words of length at most", 2 * args.nmax, "that is on k at most",1064 2 * args.nmax + 1, "so no such transfer with fewer than that many states exists there, and whether "1065 "the rank is unbounded is observed in the deficiency column and not proved here")1066 print("elapsed", round(time.time() - t0, 1), "s")10671068def cmd_check(args):1069 for k in range(2, 6):1070 for n in range(k, min(4 * k, 21) + 1):1071 for q in (True, False):1072 a = lift_count(k, n, q)1073 b = brute_lift_count(k, n, q)1074 if a != b:1075 raise ValueError("transfer " + str(a) + " against brute " + str(b) + " at k, n, prim = " + str((k, n, q)))1076 print("the column transfer matches the brute enumeration of binary multiples, primitive and raw, "1077 "at every k = 2..5, n = k..min(4k, 21)")1078 for k in range(2, 9):1079 if lift_count(k, k) != 1 or lift_count(k, k + 1) != 1:1080 raise ValueError("the length k+1 gap fails at k = " + str(k))1081 if lift_count(k, 2 * k) != 2 ** (k - 1) + 1:1082 raise ValueError("the depth-2 count fails at k = " + str(k))1083 if lift_count(k, 2 * k, False) != 2 ** k + 1:1084 raise ValueError("the raw depth-2 count fails at k = " + str(k))1085 if lift_count(k, 3 * k, False) != 2 * 3 ** k + 1:1086 raise ValueError("the raw depth-3 count fails at k = " + str(k))1087 if lift_count(k, 3 * k) != 3 ** k + 1:1088 raise ValueError("the depth-3 count fails at k = " + str(k))1089 print("L(k, k) = L(k, k+1) = 1, L(k, 2k) = 2^(k-1) + 1 and L(k, 3k) = 3^k + 1 at every k = 2..8, the raw "1090 "counts 2^k + 1 at 3^(2k) and 2 * 3^k + 1 at 3^(3k) reproducing the block ladder")1091 for k in range(2, 6):1092 n = 3 * k1093 ref = brute_depths(k, n)1094 w = (3 ** k - 1) // 21095 for z1 in range(3, w, 3):1096 if gcd(z1, w) != 1:1097 continue1098 d = first_return(z1, w - z1)1099 r = ref.get(z1, 0)1100 if r and d != r:1101 raise ValueError("the BFS return " + str(d) + " against the brute lift " + str(r) + " at " + str((k, z1)))1102 if not r and d and d <= n:1103 raise ValueError("the BFS returns at " + str(d) + " where no lift does, at " + str((k, z1)))1104 print("k", k, "brute lifts to 3^" + str(n), "depths matched on", len(ref), "directions")1105 print("the BFS first return equals the shortest binary lift on every direction the brute enumeration reaches")1106 for k in range(2, 10):1107 hist, _ = sweep(k)1108 if hist and min(hist) < k:1109 raise ValueError("a return shorter than k at k = " + str(k))1110 if any(v % 2 for v in hist.values()):1111 raise ValueError("a first-return count is odd at k = " + str(k))1112 print("no return is shorter than k at any k = 2..9, and every first-return count is even, the sweep "1113 "counting the direction (z, R_k - z) once at z and once at R_k - z")1114 for b in range(1, 9):1115 for k in range(3, 8):1116 if set(slot_profile(k, b * k)) != {b}:1117 raise ValueError("the slot profile is not uniform at " + str((k, b)))1118 if block_count(b, k) != lift_count(k, b * k):1119 raise ValueError("the fixed matrix misses the transfer at " + str((k, b)))1120 print("the slot profile is uniform and the fixed matrix reproduces the column transfer at every b = 1..8, "1121 "k = 3..7")1122 for b, mp in enumerate(MINPOLY, 1):1123 P, Q, iv, ratio, sub = lam_block(b)1124 if Q != mp:1125 raise ValueError("the minimal polynomial of lam_" + str(b) + " reads " + str(Q))1126 d = 3 * iv[0] - 2 ** b, 3 * iv[1] - 2 ** b1127 lo, hi = (-1, 2) if not b % 2 else (-2, 1)1128 if d[0] < lo or d[1] > hi:1129 raise ValueError("lam_" + str(b) + " leaves the parity bracket")1130 if lam_block(4)[0] != [90, -153, 77, -15, 1] or lam_block(5)[0] != [-45, 60, -16, 1]:1131 raise ValueError("the characteristic polynomial at b = 4 or b = 5 moved")1132 print("lam_b has the pinned minimal polynomial at every b = 1..10, lam_4 = 6 with characteristic polynomial "1133 "(x - 1)(x - 3)(x - 5)(x - 6) and lam_5 = 3 (5 + sqrt 5) / 2 with (x - 1)(x^2 - 15x + 45), and "1134 "3 lam_b - 2^b sits inside [-1, 2] at even b and inside [-2, 1] at odd b at every one")1135 for k in range(2, 10):1136 for w in range(1 << (k - 1)):1137 word = tuple((w >> i) & 1 for i in range(k - 1))1138 D, n = submask_set(word)1139 if len(D) != submask_count(word):1140 raise ValueError("the meet in the middle misses the brute submask count at " + str((k, w)))1141 if closure_faults(D, n):1142 raise ValueError("the solution set is not closed at " + str((k, w)))1143 if len(D) % 2:1144 raise ValueError("an odd submask count at " + str((k, w)))1145 print("the meet-in-the-middle submask count matches the brute enumeration, and the solution set is closed "1146 "under complement in K, under disjoint union and under nested difference, over every one of the "1147 "2^(k-1) sets T at every k = 2..9, so every count is even")1148 for k, ref in ((11, 5224), (12, 11852), (13, 20888), (14, 43364)):1149 M = sum(submask_count(tuple((w >> i) & 1 for i in range(k - 1))) for w in range(1 << (k - 1)))1150 if M != ref:1151 raise ValueError("M_" + str(k) + " reads " + str(M))1152 print("M_k reads 5224, 11852, 20888, 43364 at k = 11..14, the cut-free aggregate of the lift half")1153 D, n = submask_set((1, 0, 0, 0))1154 I = irreducibles(D)1155 if len(D) != 6 or sorted(I) != [5, 11, 20, 26]:1156 raise ValueError("the k = 5, T = {1} irreducibles moved: " + str((len(D), I)))1157 print("at k = 5, T = {1}, m = 7 the support {0, 2, 3, 4, 6} has four irreducibles and two decompositions "1158 "of the whole, 6 solutions against 7 packings, so the decomposition is not unique and the count is "1159 "the number of distinct unions and not the number of packings")1160 D, n = submask_set((0, 1, 0, 0, 0, 0))1161 I = irreducibles(D)1162 if len(D) != 8 or len(I) != 6 or all(not I[a] & I[b] for a in range(len(I)) for b in range(a)):1163 raise ValueError("the k = 7, T = {2} witness moved: " + str((len(D), len(I))))1164 print("at k = 7, T = {2}, m = 19 the count is 8, a power of two, with 6 irreducibles that overlap, so a "1165 "power-of-two count does not force pairwise disjointness")1166 r, raw, rdf, pool, pdf = geom_fit(survival(sweep(13)[0], 13))1167 if abs(r - 0.8868) > 5e-4 or abs(raw - 143.5) > 0.2 or rdf != 30 or abs(pool - 47.9) > 0.2 or pdf != 12:1168 raise ValueError("the geometric fit at k = 13 moved: " + str((r, raw, rdf, pool, pdf)))1169 print("the maximum-likelihood geometric at k = 13 keeps its ratio 0.8868 and its chi2 143.5 on 30 df "1170 "unpooled, 47.9 on 12 df pooled")11711172def main():1173 p = argparse.ArgumentParser()1174 s = p.add_subparsers(dest="cmd", required=True)1175 a = s.add_parser("automaton")1176 a.add_argument("--weights", type=int, nargs="+", default=[13, 40, 100, 101, 121, 257, 364, 1093])1177 a.set_defaults(fn=cmd_automaton)1178 b = s.add_parser("returns")1179 b.add_argument("--kmin", type=int, default=3)1180 b.add_argument("--kmax", type=int, default=24)1181 b.add_argument("--bmax", type=int, default=8)1182 b.add_argument("--kfine", type=int, default=8)1183 b.add_argument("--kbig", type=int, default=160)1184 b.set_defaults(fn=cmd_returns)1185 c = s.add_parser("hist")1186 c.add_argument("--kmin", type=int, default=2)1187 c.add_argument("--kmax", type=int, default=13)1188 c.add_argument("--bulkmin", type=int, default=9)1189 c.set_defaults(fn=cmd_hist)1190 e = s.add_parser("model")1191 e.add_argument("--kmin", type=int, default=6)1192 e.add_argument("--kmax", type=int, default=15)1193 e.add_argument("--bmax", type=int, default=8)1194 e.set_defaults(fn=cmd_model)1195 f = s.add_parser("critical")1196 f.add_argument("--wmin", type=int, default=2000)1197 f.add_argument("--wmax", type=int, default=400000)1198 f.add_argument("--step", type=float, default=1.35)1199 f.set_defaults(fn=cmd_critical)1200 g = s.add_parser("ladder")1201 g.add_argument("--bmin", type=int, default=1)1202 g.add_argument("--bmax", type=int, default=14)1203 g.add_argument("--kcheck", type=int, default=7)1204 g.set_defaults(fn=cmd_ladder)1205 h = s.add_parser("lift")1206 h.add_argument("--kmax", type=int, default=12)1207 h.add_argument("--kbig", type=int, default=17)1208 h.add_argument("--nmax", type=int, default=7)1209 h.set_defaults(fn=cmd_lift)1210 d = s.add_parser("check")1211 d.set_defaults(fn=cmd_check)1212 args = p.parse_args()1213 args.fn(args)12141215main()