question.py
31.6 kB · python · 853 lines
1import argparse2import math3import time4import urllib.error5import urllib.request6from fractions import Fraction78import numpy as np9from scipy.special import zeta as hurwitz1011OEIS = "https://raw.githubusercontent.com/oeis/oeisdata/main/seq/A002/A002487.seq"12E2 = "0.5312805062772051416"13LN2 = math.log(2.0)14PHI = (1.0 + 5.0 ** 0.5) / 2.015HOLDER = LN2 / (2.0 * math.log(PHI))16M0 = ((1, 1), (0, 1))17M1 = ((1, 0), (1, 1))1819# STERN2021def stern(top):22 s = [0, 1]23 for n in range(2, top + 1):24 s.append(s[n // 2] if n % 2 == 0 else s[n // 2] + s[n // 2 + 1])25 return s262728def mul(v, m):29 return (v[0] * m[0][0] + v[1] * m[1][0], v[0] * m[0][1] + v[1] * m[1][1])303132def mm(a, b):33 return tuple(tuple(sum(a[i][t] * b[t][j] for t in range(2)) for j in range(2)) for i in range(2))343536def oeis():37 request = urllib.request.Request(OEIS, headers={"User-Agent": "curl/8"})38 try:39 with urllib.request.urlopen(request, timeout=10) as handle:40 text = handle.read().decode()41 except (urllib.error.URLError, OSError, TimeoutError):42 return None43 terms = []44 for line in text.splitlines():45 if line[:2] in ("%S", "%T", "%U"):46 terms += [int(x) for x in line.split(" ", 2)[2].strip().strip(",").split(",")]47 return terms484950def brocot_row(depth, unit):51 frontier = [((0, 1), (1, 1) if unit else (1, 0))]52 for _ in range(depth + 1):53 nxt = []54 row = []55 for lo, hi in frontier:56 mid = (lo[0] + hi[0], lo[1] + hi[1])57 row.append(mid)58 nxt += [(lo, mid), (mid, hi)]59 frontier = nxt60 return [Fraction(p, q) for p, q in row]616263def carry(span, depth):64 t0 = time.time()65 s = stern((1 << (span + 1)) + 2)66 bad = 067 for n in range(1 << span):68 v = (s[n], s[n + 1])69 if mul(v, M0) != (s[2 * n], s[2 * n + 1]) or mul(v, M1) != (s[2 * n + 1], s[2 * n + 2]):70 bad += 171 print(f"carry: v(2n+b) = v(n) M_b mismatches below 2^{span}: {bad}")72 bad = 073 for n in range(1, 1 << span):74 v = (0, 1)75 for bit in bin(n)[2:]:76 v = mul(v, M1 if bit == "1" else M0)77 if v != (s[n], s[n + 1]):78 bad += 179 print(f"carry: v(n) = (0,1) M_(d_1) ... M_(d_L) mismatches below 2^{span}: {bad}")80 ref = oeis()81 if ref is None:82 print("carry: A002487 fetch skipped, offline")83 else:84 print(f"carry: A002487 listed terms {len(ref)}, agree: {ref == s[:len(ref)]}")85 print(f"carry: M_0 M_1 = {mm(M0, M1)}, M_1 M_0 = {mm(M1, M0)}, equal: {mm(M0, M1) == mm(M1, M0)}")86 print(f"carry: s(3), s(5), s(6) = {s[3]}, {s[5]}, {s[6]} on the words 011, 101, 110")87 for lead, low in (("without", 1), ("with", 0)):88 least = None89 for L in range(2, 12):90 lo = (1 << (L - 1)) if low else 091 for n in range(lo, 1 << L):92 for m in range(n + 1, 1 << L):93 if bin(n).count("1") == bin(m).count("1") and s[n] != s[m]:94 least = (L, format(n, f"0{L}b"), format(m, f"0{L}b"), s[n], s[m])95 break96 if least:97 break98 if least:99 break100 print(f"carry: least pair of one length and one weight with different s, {lead} leading zeros: {least}")101 print(f"carry: level 2 cells read {s[:4]}")102 worst = 0103 for d in range(depth + 1):104 row = brocot_row(d, False)105 cw = [Fraction(s[n], s[n + 1]) for n in range((1 << d), (1 << (d + 1)))]106 rev = [int(format(i, f"0{d}b")[::-1], 2) if d else 0 for i in range(1 << d)]107 worst = max(worst, sum(1 for i in range(1 << d) if row[i] != cw[rev[i]]))108 print(f"carry: Stern-Brocot row against Calkin-Wilf under bit reversal, depths 0..{depth}, mismatches: {worst}")109 print(f"carry: {time.time() - t0:.2f} s")110111112# QUESTION MARK113114def quotients(x):115 out = []116 while x:117 q = 1 / x118 a = q.numerator // q.denominator119 out.append(a)120 x = q - a121 return out122123124def denjoy(x):125 y = Fraction(0)126 total = 0127 for i, a in enumerate(quotients(x)):128 total += a129 y += (-1) ** i * Fraction(1, 2 ** (total - 1))130 return y131132133def farey(x):134 return x / (1 - x) if x <= Fraction(1, 2) else (1 - x) / x135136137def tent(y):138 return 2 * y if y <= Fraction(1, 2) else 2 - 2 * y139140141def runs_word(qs, length):142 word = "0" * (qs[0] - 1)143 bit = "1"144 for a in qs[1:]:145 word += bit * a146 bit = "1" if bit == "0" else "0"147 return word[:length] if len(word) >= length else None148149150def question(depth, level, alphabets):151 t0 = time.time()152 q = {Fraction(0): Fraction(0), Fraction(1): Fraction(1)}153 frontier = [(Fraction(0), Fraction(1))]154 for _ in range(depth + 1):155 nxt = []156 for a, b in frontier:157 m = Fraction(a.numerator + b.numerator, a.denominator + b.denominator)158 q[m] = (q[a] + q[b]) / 2159 nxt += [(a, m), (m, b)]160 frontier = nxt161 bad_denjoy = sum(1 for x in q if denjoy(x) != q[x])162 print(f"question: mediant recursion against Denjoy on {len(q)} nodes to depth {depth}: mismatches {bad_denjoy}")163 bad_tent = sum(1 for x in q if 0 < x < 1 and denjoy(farey(x)) != tent(q[x]))164 print(f"question: ?(F(x)) = T(?(x)) on the same nodes: mismatches {bad_tent}")165 bad_branch = 0166 for x in q:167 for a in range(1, 6):168 if denjoy(1 / (a + x)) != Fraction(2 - q[x], 2 ** a):169 bad_branch += 1170 print(f"question: ?(1/(a+x)) = 2^(-a)(2 - ?(x)) for a = 1..5 on the same nodes: mismatches {bad_branch}")171 bad_addr = 0172 for d in range(depth + 1):173 for i, x in enumerate(brocot_row(d, True)):174 if q[x] != Fraction(2 * i + 1, 2 ** (d + 1)):175 bad_addr += 1176 print(f"question: ?(node i of row d) = (2i+1)/2^(d+1), depths 0..{depth}: mismatches {bad_addr}")177 for k, alpha in alphabets:178 code = code_of(k, alpha)179 accepted = {w for w in words(level) if accepts(code, k, w)}180 prefixes = prefixes_of(alpha, k, level)181 extra = accepted - prefixes182 missing = prefixes - accepted183 heads = sorted({(w[0], len(w) - len(w.lstrip(w[0]))) for w in extra})184 print(f"question: width {k} code {code} alphabet {name(alpha)}: accepted {len(accepted)}, prefixes of ?(E_A) {len(prefixes)}, missing {len(missing)}, extra {len(extra)} with leading runs {heads}")185 print(f"question: {time.time() - t0:.2f} s")186187188# RULES189190def words(length):191 return [format(i, f"0{length}b") for i in range(1 << length)]192193194def accepts(code, k, w):195 return all((code >> int(w[i:i + k], 2)) & 1 for i in range(len(w) - k + 1))196197198def forbidden(alpha, k):199 pats = []200 if alpha[0] == "finite":201 m = max(alpha[1])202 pats += ["0" * (m + 1), "1" * (m + 1)]203 gaps = [j for j in range(1, m) if j not in alpha[1]]204 else:205 gaps = sorted(alpha[1])206 for j in gaps:207 pats += ["1" + "0" * j + "1", "0" + "1" * j + "0"]208 assert all(len(p) <= k for p in pats)209 return pats210211212def code_of(k, alpha):213 pats = forbidden(alpha, k)214 code = 0215 for w in range(1 << k):216 ws = format(w, f"0{k}b")217 if not any(p in ws for p in pats):218 code |= 1 << w219 return code220221222def name(alpha):223 if alpha[0] == "finite":224 return "{" + ",".join(str(a) for a in alpha[1]) + "}"225 return "N" if not alpha[1] else "N\\{" + ",".join(str(a) for a in alpha[1]) + "}"226227228def alphabets_at(k):229 out = []230 for mask in range(1, 1 << (k - 1)):231 out.append(("finite", tuple(j + 1 for j in range(k - 1) if (mask >> j) & 1)))232 for mask in range(1 << (k - 2)):233 out.append(("cofinite", tuple(j + 1 for j in range(k - 2) if (mask >> j) & 1)))234 return out235236237def prefixes_of(alpha, k, level):238 if alpha[0] == "finite":239 letters = list(alpha[1])240 else:241 letters = [a for a in range(1, level + 2) if a not in alpha[1]]242 out = set()243 stack = [((), 0)]244 while stack:245 qs, total = stack.pop()246 if total >= level + 1 and len(qs) >= 2:247 w = runs_word(qs, level)248 if w is not None:249 out.add(w)250 continue251 for a in letters:252 stack.append((qs + (a,), total + a))253 return out254255256def transfer(code, k):257 n = 1 << (k - 1)258 m = np.zeros((n, n))259 for s in range(n):260 for c in range(2):261 w = s * 2 + c262 if (code >> w) & 1:263 m[s, w % n] = 1.0264 return m265266267def perron(code, k):268 ev = np.linalg.eigvals(transfer(code, k))269 return float(max(ev.real))270271272def nacci(alpha):273 if alpha[0] == "finite":274 m = max(alpha[1])275 p = [0] * (m + 1)276 p[m] = 1277 for a in alpha[1]:278 p[m - a] -= 1279 return p280 f = max(alpha[1]) if alpha[1] else 0281 p = [0] * (f + 2)282 p[f + 1] += 1283 p[f] -= 2284 for j in alpha[1]:285 p[f - j + 1] += 1286 p[f - j] -= 1287 return p288289290def largest_root(p):291 r = np.roots(p[::-1])292 return float(max(z.real for z in r if abs(z.imag) < 1e-9))293294295def charpoly_int(m):296 n = m.shape[0]297 a = np.array(np.rint(m).astype(int), dtype=object)298 cs = [1]299 mj = a.copy()300 ident = np.array([[1 if i == j else 0 for j in range(n)] for i in range(n)], dtype=object)301 for j in range(1, n + 1):302 tr = sum(mj[i, i] for i in range(n))303 c = -tr // j304 cs.append(c)305 mj = a.dot(mj + c * ident)306 return cs307308309def divides(p, q):310 q = list(q)311 p = list(p)312 while q and q[-1] == 0:313 q.pop()314 while len(q) >= len(p):315 f = q[-1]316 if f % p[-1]:317 return False318 f //= p[-1]319 d = len(q) - len(p)320 for i, c in enumerate(p):321 q[d + i] -= f * c322 while q and q[-1] == 0:323 q.pop()324 return not q325326327# TRANSFER OPERATOR328329def nodes(n):330 x = (1.0 - np.cos(np.pi * np.arange(n) / (n - 1))) / 2.0331 w = (-1.0) ** np.arange(n)332 w[0] /= 2.0333 w[-1] /= 2.0334 return x, w335336337def bary_rows(t, x, w):338 d = t[:, None] - x[None, :]339 hit = np.abs(d) < 1e-15340 r = w[None, :] / np.where(hit, 1.0, d)341 r = np.where(hit.any(axis=1)[:, None], hit.astype(float), r)342 return r / r.sum(axis=1, keepdims=True)343344345def diff_matrix(x, w):346 n = len(x)347 d = np.zeros((n, n))348 for j in range(n):349 for i in range(n):350 if i != j:351 d[j, i] = (w[i] / w[j]) / (x[j] - x[i])352 d[j, j] = -d[j].sum()353 return d354355356def operator(s, alpha, x, w, cut, taylor):357 n = len(x)358 m = np.zeros((n, n))359 if alpha[0] == "finite":360 letters = list(alpha[1])361 else:362 letters = [a for a in range(1, cut + 1) if a not in alpha[1]]363 for a in letters:364 t = 1.0 / (a + x)365 m += ((a + x) ** (-2.0 * s))[:, None] * bary_rows(t, x, w)366 if alpha[0] == "cofinite":367 d = diff_matrix(x, w)368 dk = np.eye(n)369 fact = 1.0370 for k in range(taylor + 1):371 m += hurwitz(2.0 * s + k, cut + 1.0 + x)[:, None] * (dk[0][None, :] / fact)372 dk = dk.dot(d)373 fact *= k + 1374 return m375376377def leading(s, alpha, x, w, cut, taylor):378 ev = np.linalg.eigvals(operator(s, alpha, x, w, cut, taylor))379 i = int(np.argmax(np.abs(ev)))380 return ev[i].real, abs(ev[i].imag)381382383def pressure_zero(alpha, modes, cut, taylor):384 x, w = nodes(modes)385 if alpha[0] == "finite" and len(alpha[1]) == 1:386 lam, _ = leading(0.0, alpha, x, w, cut, taylor)387 return 0.0, lam388 lo, hi = (0.0, 1.0) if alpha[0] == "finite" else (0.51, 1.1)389 g = lambda s: math.log(leading(s, alpha, x, w, cut, taylor)[0])390 glo, ghi = g(lo), g(hi)391 assert glo > 0 > ghi, (alpha, glo, ghi)392 for _ in range(70):393 mid = (lo + hi) / 2.0394 gm = g(mid)395 if gm > 0:396 lo = mid397 else:398 hi = mid399 if hi - lo < 1e-16:400 break401 s = (lo + hi) / 2.0402 return s, leading(s, alpha, x, w, cut, taylor)[0]403404405def control(modes, cut, taylor):406 a12 = ("finite", (1, 2))407 s, lam = pressure_zero(a12, modes, cut, taylor)408 ref = float(E2)409 digits = -math.log10(abs(s - ref)) if s != ref else 17410 print(f"control: A = {{1,2}} pressure zero {s:.16f} against {E2}, gap {abs(s - ref):.1e}, {digits:.1f} digits, leading eigenvalue {lam:.16f}")411 assert abs(s - ref) < 1e-10, "E_2 control failed"412 for m in (24, 32, 48, 56):413 s2, _ = pressure_zero(a12, m, cut, taylor)414 print(f"control: A = {{1,2}} at {m} modes: {s2:.16f}, gap to {modes} modes {abs(s2 - s):.1e}")415 sn, lam = pressure_zero(("cofinite", ()), modes, cut, taylor)416 print(f"control: A = N pressure zero {sn:.16f} against 1, gap {abs(sn - 1.0):.1e}, leading eigenvalue {lam:.16f}")417 assert abs(sn - 1.0) < 1e-10, "Gauss control failed"418 for c in (500, 1000, 4000):419 s3, _ = pressure_zero(("cofinite", (1,)), modes, c, taylor)420 print(f"control: A = N\\{{1}} at cut {c}: {s3:.16f}")421 s4, _ = pressure_zero(("cofinite", (1,)), modes, cut, taylor)422 print(f"control: A = N\\{{1}} at cut {cut}: {s4:.16f}")423 return s424425426def table(widths, modes, cut, taylor):427 t0 = time.time()428 control(modes, cut, taylor)429 named = {(2, 7): 1.618033988749, (3, 23): 1.465571231876, (3, 54): 1.324717957244, (3, 127): 1.839286755214}430 for (k, code), rho in named.items():431 assert abs(perron(code, k) - rho) < 1e-9, (k, code)432 print(f"table: transfer matrix convention pinned on the four named codes of the census: ok")433 print("| k | code | A | P_A | rho | log_2 rho | dim_CF | alpha log_2 rho |")434 print("|---|---|---|---|---|---|---|---|")435 rows = []436 for k in widths:437 for alpha in alphabets_at(k):438 code = code_of(k, alpha)439 rho = perron(code, k)440 p = nacci(alpha)441 r = largest_root(p)442 cp = charpoly_int(transfer(code, k))443 div = divides(p, cp[::-1])444 assert abs(rho - r) < 1e-9 and div, (k, code, rho, r, div)445 s, lam = pressure_zero(alpha, modes, cut, taylor)446 dy = math.log(rho) / LN2447 bound = HOLDER * dy448 assert s + 1e-12 >= bound, (k, code, s, bound)449 rows.append((k, code, name(alpha), poly_text(p), rho, dy, s, bound))450 print(f"| {k} | {code} | `{name(alpha)}` | `{poly_text(p)}` | {trunc(rho)} | {trunc(dy)} | {trunc(s)} | {trunc(bound)} |")451 seen = {}452 for k, code, a, p, rho, dy, s, b in rows:453 seen.setdefault(a, []).append((k, code, s))454 print(f"table: {len(rows)} codes on {len(seen)} alphabets; the same alphabet at two widths prints the same dim_CF: {all(max(v)[2] - min(v)[2] == 0 for v in seen.values())}")455 print(f"table: {time.time() - t0:.2f} s")456457458def trunc(x):459 return f"{math.floor(x * 1e12 + 1e-3) / 1e12:.12f}"460461462def poly_text(p):463 terms = []464 for i in range(len(p) - 1, -1, -1):465 c = p[i]466 if c == 0:467 continue468 mono = "" if i == 0 else ("x" if i == 1 else f"x^{i}")469 if abs(c) == 1 and i > 0:470 body = mono471 else:472 body = f"{abs(c)}{mono}" if i > 0 else f"{abs(c)}"473 terms.append(("- " if c < 0 else "+ ") + body)474 text = " ".join(terms)475 return text[2:] if text.startswith("+ ") else "-" + text[2:]476477478def symmetric_orphans(k):479 codes = {code_of(k, a) for a in alphabets_at(k)}480 total = 1 << (1 << k)481 sym = []482 for code in range(total):483 flip = 0484 rev = 0485 for w in range(1 << k):486 if (code >> w) & 1:487 flip |= 1 << ((1 << k) - 1 - w)488 rev |= 1 << int(format(w, f"0{k}b")[::-1], 2)489 if flip == code and rev == code:490 sym.append(code)491 orphans = [c for c in sym if c not in codes]492 live = [c for c in orphans if perron(c, k) > 1.0 + 1e-9]493 return codes, sym, orphans, live494495496def obstruction(widths):497 for k in widths:498 codes, sym, orphans, live = symmetric_orphans(k)499 total = 1 << (1 << k)500 print(f"obstruction: width {k}: {len(codes)} run-length codes of {total}, {len(sym)} codes fixed by both the digit flip and reversal, {len(orphans)} of those with no alphabet")501 print(f"obstruction: width {k}: {len(live)} of the {len(orphans)} carry rho > 1")502 for c in orphans if k <= 3 else live[:1]:503 allowed = [format(w, f"0{k}b") for w in range(1 << k) if (c >> w) & 1]504 forbid = [format(w, f"0{k}b") for w in range(1 << k) if not (c >> w) & 1]505 print(f"obstruction: width {k} code {c} allows {allowed if k <= 3 else len(allowed)} forbids {forbid}, rho {trunc(perron(c, k))}")506507508# GRAPH509510def graph_of(code, k):511 n = 1 << (k - 1)512 edges = {}513 tails = {}514 for u in range(n):515 b = 1 - (u & 1)516 s = u517 for a in range(1, k + 1):518 w = 2 * s + b519 if not (code >> w) & 1:520 break521 s = w % n522 if a < k:523 edges.setdefault((u, s), []).append(a)524 else:525 tails[u] = s526 return edges, tails527528529def components(n, edges, tails):530 succ = {u: set() for u in range(n)}531 for u, v in edges:532 succ[u].add(v)533 for u, v in tails.items():534 succ[u].add(v)535 reach = [[False] * n for _ in range(n)]536 for u in range(n):537 stack = list(succ[u])538 while stack:539 v = stack.pop()540 if not reach[u][v]:541 reach[u][v] = True542 stack += list(succ[v])543 comp = {}544 for u in range(n):545 if reach[u][u]:546 comp[u] = frozenset(v for v in range(n) if reach[u][v] and reach[v][u])547 return comp548549550def branch_table(x, w, cut):551 return np.array([bary_rows(1.0 / (a + x), x, w) for a in range(1, cut + 1)])552553554def taylor_rows(x, w, taylor):555 d = diff_matrix(x, w)556 dk = np.eye(len(x))557 fact = 1.0558 rows = []559 for j in range(taylor + 1):560 rows.append(dk[0] / fact)561 dk = dk.dot(d)562 fact *= j + 1563 return np.array(rows)564565566def graph_operator(s, k, comp, edges, tails, x, B, rows, cut):567 n = len(x)568 states = sorted(comp)569 idx = {u: i for i, u in enumerate(states)}570 m = np.zeros((n * len(states), n * len(states)))571 W = (np.arange(1, cut + 1)[:, None] + x[None, :]) ** (-2.0 * s)572 small = {a: W[a - 1][:, None] * B[a - 1] for a in range(1, k)}573 tail = None574 for (u, v), labels in edges.items():575 if u in comp and v in comp[u]:576 for a in labels:577 m[idx[u] * n:(idx[u] + 1) * n, idx[v] * n:(idx[v] + 1) * n] += small[a]578 for u, v in tails.items():579 if u in comp and v in comp[u]:580 if tail is None:581 tail = np.einsum("ai,aij->ij", W[k - 1:], B[k - 1:])582 for j in range(rows.shape[0]):583 tail += hurwitz(2.0 * s + j, cut + 1.0 + x)[:, None] * rows[j][None, :]584 m[idx[u] * n:(idx[u] + 1) * n, idx[v] * n:(idx[v] + 1) * n] += tail585 return m586587588def graph_shape(code, k):589 n = 1 << (k - 1)590 edges, tails = graph_of(code, k)591 comp = components(n, edges, tails)592 live = any(u in comp and v in comp[u] for u, v in tails.items())593 out = {u: 0 for u in comp}594 for (u, v), labels in edges.items():595 if u in comp and v in comp[u]:596 out[u] += len(labels)597 for u, v in tails.items():598 if u in comp and v in comp[u]:599 out[u] += 2600 return edges, tails, comp, live, max(out.values(), default=0)601602603def graph_zero(code, k, x, B, rows, cut):604 edges, tails, comp, live, fan = graph_shape(code, k)605 g = lambda s: math.log(max(np.linalg.eigvals(graph_operator(s, k, comp, edges, tails, x, B, rows, cut)).real))606 if fan <= 1:607 return 0.0, math.exp(g(0.0)) if comp else 0.0608 lo, hi = (0.51, 1.1) if live else (0.0, 1.0)609 glo, ghi = g(lo), g(hi)610 assert glo > 0 > ghi, (code, glo, ghi)611 for _ in range(70):612 mid = (lo + hi) / 2.0613 if g(mid) > 0:614 lo = mid615 else:616 hi = mid617 if hi - lo < 1e-16:618 break619 s = (lo + hi) / 2.0620 return s, math.exp(g(s))621622623def markov_prefixes(letters, pairs, level):624 out = set()625 stack = [((), 0)]626 while stack:627 qs, total = stack.pop()628 if total >= level + 1 and len(qs) >= 2:629 w = runs_word(qs, level)630 if w is not None:631 out.add(w)632 continue633 for a in letters:634 if qs and (qs[-1], a) in pairs:635 continue636 stack.append((qs + (a,), total + a))637 return out638639640def graph(widths, modes, cut, taylor, level):641 t0 = time.time()642 x, w = nodes(modes)643 B = branch_table(x, w, cut)644 rows = taylor_rows(x, w, taylor)645 worst = 0.0646 for k in widths:647 for alpha in alphabets_at(k):648 code = code_of(k, alpha)649 s_graph, _ = graph_zero(code, k, x, B, rows, cut)650 s_alpha, _ = pressure_zero(alpha, modes, cut, taylor)651 worst = max(worst, abs(s_graph - s_alpha))652 assert abs(s_graph - s_alpha) < 1e-12, (k, code, s_graph, s_alpha)653 print(f"graph: the graph form recovers the pressure zero of all 18 run-length codes, largest gap {worst:.1e}")654 k = 4655 code = 11892656 accepted = {v for v in words(level) if accepts(code, k, v)}657 prefixes = markov_prefixes((1, 2), {(2, 2)}, level)658 extra = accepted - prefixes659 heads = sorted({(v[0], len(v) - len(v.lstrip(v[0]))) for v in extra})660 print(f"graph: code {code} at level {level}: accepted {len(accepted)}, prefixes of ?(M) for M = {{1,2}} without the pair 22: {len(prefixes)}, missing {len(prefixes - accepted)}, extra {len(extra)} with leading runs {heads}")661 edges, tails, comp, live, fan = graph_shape(code, k)662 print(f"graph: code {code} graph: recurrent states {sorted(format(u, '03b') for u in comp)}, edges {sorted((format(u, '03b'), format(v, '03b'), tuple(l)) for (u, v), l in edges.items() if u in comp and v in comp[u])}, tail edges {live}")663 _, _, _, live_codes = symmetric_orphans(k)664 order = [code] + [c for c in live_codes if c != code]665 print("| code | forbids | states | tail | rho | log_2 rho | dim_CF | alpha log_2 rho |")666 print("|---|---|---|---|---|---|---|---|")667 for c in order:668 rho = perron(c, k)669 dy = math.log(rho) / LN2670 s, _ = graph_zero(c, k, x, B, rows, cut)671 bound = HOLDER * dy672 assert s + 1e-12 >= bound, (c, s, bound)673 edges, tails, comp, live, fan = graph_shape(c, k)674 forbid = " ".join(format(v, f"0{k}b") for v in range(1 << k) if not (c >> v) & 1)675 print(f"| {c} | `{forbid}` | {len(comp)} | {'yes' if live else 'no'} | {trunc(rho)} | {trunc(dy)} | {trunc(s)} | {trunc(bound)} |")676 print("| code | rule | CF constraint | rho | log_2 rho | dim_CF | alpha log_2 rho |")677 print("|---|---|---|---|---|---|---|")678 named = ((2, 7, "no 11", "even quotients 1"), (3, 23, "at most one 1 per 3", "even quotients 1, odd quotients at least 2 past the first"), (3, 54, "no 11, no 000", "even quotients 1, odd quotients 1 or 2 past the first"), (3, 127, "no 111", "even quotients 1 or 2"))679 for k2, c, rule, cf in named:680 rho = perron(c, k2)681 dy = math.log(rho) / LN2682 s, _ = graph_zero(c, k2, x, B, rows, cut)683 bound = HOLDER * dy684 assert s + 1e-12 >= bound, (c, s, bound)685 print(f"| {c} | {rule} | {cf} | {trunc(rho)} | {trunc(dy)} | {trunc(s)} | {trunc(bound)} |")686 print(f"graph: {time.time() - t0:.2f} s")687688689# SUBLEADING690691def complex_parts(modes, letters):692 x, w = nodes(modes)693 return [np.log(a + x) for a in letters], [bary_rows(1.0 / (a + x), x, w) for a in letters]694695696def determinant(s, logs, B):697 n = len(logs[0])698 m = np.zeros((n, n), dtype=complex)699 for lg, b in zip(logs, B):700 m += np.exp(-2.0 * s * lg)[:, None] * b701 return np.linalg.det(np.eye(n) - m), m702703704def newton(s, logs, B):705 h = 1e-6706 for _ in range(60):707 f, _ = determinant(s, logs, B)708 fp = (determinant(s + h, logs, B)[0] - determinant(s - h, logs, B)[0]) / (2 * h)709 step = f / fp710 s = s - step711 if abs(step) < 1e-15:712 break713 return s714715716def winding(logs, B, s0, s1, t0, t1, n=3000):717 pts = [complex(s0 + (s1 - s0) * k / n, t0) for k in range(n)]718 pts += [complex(s1, t0 + (t1 - t0) * k / n) for k in range(n)]719 pts += [complex(s1 - (s1 - s0) * k / n, t1) for k in range(n)]720 pts += [complex(s0, t1 - (t1 - t0) * k / n) for k in range(n)]721 ph = np.angle(np.array([determinant(p, logs, B)[0] for p in pts]))722 d = np.diff(np.concatenate([ph, ph[:1]]))723 d = (d + np.pi) % (2 * np.pi) - np.pi724 return int(round(d.sum() / (2 * np.pi)))725726727def zeros_in(logs, B, s0, s1, ds, t0, t1, dt):728 sig = np.arange(s0, s1 + ds / 2, ds)729 tau = np.arange(t0, t1 + dt / 2, dt)730 F = np.array([[abs(determinant(complex(a, b), logs, B)[0]) for b in tau] for a in sig])731 mins = sorted((F[i, j], sig[i], tau[j]) for i in range(1, len(sig) - 1) for j in range(1, len(tau) - 1) if F[i, j] == F[i - 1:i + 2, j - 1:j + 2].min())732 zeros = []733 for _, a, b in mins:734 z = newton(complex(a, b), logs, B)735 inside = s0 <= z.real <= s1 and t0 <= z.imag <= t1736 if inside and abs(determinant(z, logs, B)[0]) < 1e-10 and not any(abs(z - q) < 1e-8 for q in zeros):737 zeros.append(z)738 return sorted(zeros, key=lambda z: -z.real)739740741def census_of(m, jmax):742 import importlib.util743 import os744 path = os.path.join(os.path.dirname(os.path.abspath(__file__)), "..", "ford-horocycle", "ford_horocycle.py")745 spec = importlib.util.spec_from_file_location("ford_horocycle", path)746 module = importlib.util.module_from_spec(spec)747 spec.loader.exec_module(module)748 return module.census_walk(m, 2 ** jmax)749750751def fit_zero(cum, delta, z, jlo, jmax, per):752 js = np.arange(jlo * per, jmax * per + 1) / per753 Q = np.array([int(round(2.0 ** j)) for j in js])754 R = np.array([int(cum[q]) for q in Q], dtype=float) / Q ** (2 * delta)755 lnQ = np.log(Q)756 env = Q ** (2 * (z.real - delta))757 X0 = np.ones((len(Q), 1))758 X1 = np.column_stack([np.ones(len(Q)), env * np.cos(2 * z.imag * lnQ), env * np.sin(2 * z.imag * lnQ)])759 out = []760 for X in (X0, X1):761 coef, *_ = np.linalg.lstsq(X, R, rcond=None)762 res = R - X @ coef763 out.append((coef, math.sqrt(float(np.mean(res ** 2))), float(np.abs(res).max())))764 return Q, R, out765766767def subleading(modes, jmax, jlo, per):768 t0 = time.time()769 delta = float(E2)770 letters = (1, 2)771 logs, B = complex_parts(modes, letters)772 small = zeros_in(logs, B, 0.0, 0.53, 0.01, 0.2, 14.0, 0.05)773 assert len(small) == 1 and winding(logs, B, 0.0, 0.53, 0.2, 14.0) == 1, small774 z = small[0]775 _, m = determinant(z, logs, B)776 ev = np.linalg.eigvals(m)777 near = ev[int(np.argmin(np.abs(ev - 1.0)))]778 print(f"subleading: A = {{1,2}}, one zero of det(1 - L_s) in 0 <= sigma <= 0.53, 0.2 <= tau <= 14 at {modes} modes, the winding number agreeing: s_1 = {z.real:.12f} + {z.imag:.12f} i, eigenvalue nearest 1 there {near.real:.12f} {near.imag:+.1e} i")779 for m2 in (60, 100, 140):780 logs2, B2 = complex_parts(m2, letters)781 z2 = newton(z, logs2, B2)782 print(f"subleading: s_1 at {m2} modes {z2.real:.15f} + {z2.imag:.15f} i, gap to {modes} modes {abs(z2 - z):.1e}")783 logs2, B2 = complex_parts(100, letters)784 wide = zeros_in(logs2, B2, 0.0, 0.53, 0.01, 0.2, 80.0, 0.1)785 count = winding(logs2, B2, 0.0, 0.53, 0.2, 80.0)786 assert len(wide) == count and abs(wide[0] - z) < 1e-9, (len(wide), count)787 print(f"subleading: {len(wide)} zeros in 0 <= sigma <= 0.53, 0.2 <= tau <= 80 at 100 modes, the winding number agreeing; the three of largest real part: " + ", ".join(f"{q.real:.9f} + {q.imag:.9f} i" for q in wide[:3]))788 period = math.pi / (z.imag * LN2)789 ratio = 2.0 ** (2.0 * (z.real - delta))790 alias = 2.0 * z.imag * LN2 - 3.0 * math.pi791 print(f"subleading: delta_2 - sigma_1 = {delta - z.real:.6f}; a zero s = sigma + i tau puts Q^(2 sigma) cos(2 tau log Q) into N(Q), so the period in octaves is pi/(tau log 2) = {period:.4f} and the amplitude ratio per octave 2^(2(sigma - delta_2)) = {ratio:.4f}; per octave the phase advances 2 tau log 2 = 3 pi + {alias:.4f}, so at integer octaves the wave reads as a sign alternation under an envelope of period {2 * math.pi / alias:.1f} octaves")792 t1 = time.time()793 cum = census_of(2, jmax)794 print(f"subleading: E_2 census walk to 2^{jmax} by lab/py/ford-horocycle census_walk, N_2(2^{jmax}) = {int(cum[2 ** jmax])} ({time.time() - t1:.1f}s)")795 Q, R, (flat, one) = fit_zero(cum, delta, z, jlo, jmax, per)796 C, b1, b2 = one[0]797 print(f"subleading: N_2(Q)/Q^(2 delta_2) on {len(Q)} points, {per} per octave from 2^{jlo} to 2^{jmax}: constant fit C = {flat[0][0]:.6f} leaves rms {flat[1]:.2e}, max {flat[2]:.2e}; constant plus the wave of s_1 with only C, amplitude and phase fitted: C = {C:.6f}, amplitude {math.hypot(b1, b2):.6f}, phase {math.atan2(-b2, b1):.4f}, rms {one[1]:.2e}, max {one[2]:.2e}")798 fit = lambda q: C + q ** (2 * (z.real - delta)) * (b1 * math.cos(2 * z.imag * math.log(q)) + b2 * math.sin(2 * z.imag * math.log(q)))799 rows = []800 for j in range(jlo, jmax):801 q = 2 ** j802 obs = math.log2(int(cum[2 * q]) / int(cum[q])) - 2 * delta803 rows.append((j, obs, math.log2(fit(2 * q) / fit(q))))804 print("subleading: octave exponent minus 2 delta_2, observed against the fitted wave of s_1: " + ", ".join(f"j = {j}: {o:+.4f} / {p:+.4f}" for j, o, p in rows))805 assert all(abs(o - p) < 4e-3 for _, o, p in rows), rows806 from scipy.optimize import least_squares807 lnQ = np.log(Q)808 resid = lambda p: p[2] + Q ** (2 * (p[0] - delta)) * (p[3] * np.cos(2 * p[1] * lnQ) + p[4] * np.sin(2 * p[1] * lnQ)) - R809 sol = least_squares(resid, [z.real, z.imag, C, b1, b2])810 print(f"subleading: sigma and tau fitted freely on the same points: {sol.x[0]:.4f}, {sol.x[1]:.4f} against the operator's {z.real:.4f}, {z.imag:.4f}, rms {math.sqrt(float(np.mean(sol.fun ** 2))):.2e}")811 delta3 = 0.705660908028812 logs3, B3 = complex_parts(100, (1, 2, 3))813 three = zeros_in(logs3, B3, 0.0, 0.71, 0.01, 0.2, 80.0, 0.05)814 count3 = winding(logs3, B3, 0.0, 0.71, 0.2, 80.0)815 assert len(three) == count3, (len(three), count3)816 z3 = three[0]817 cum3 = census_of(3, 18)818 _, _, (flat3, one3) = fit_zero(cum3, delta3, z3, 10, 18, per)819 print(f"subleading: control A = {{1,2,3}}: {len(three)} zeros in 0 <= sigma <= 0.71, 0.2 <= tau <= 80 at 100 modes, the winding number agreeing, the largest real part at {z3.real:.12f} + {z3.imag:.12f} i, delta_3 - sigma = {delta3 - z3.real:.4f}, ratio per octave {2.0 ** (2.0 * (z3.real - delta3)):.4f}; N_3(Q)/Q^(2 delta_3) from 2^10 to 2^18: constant fit rms {flat3[1]:.2e}, with the wave {one3[1]:.2e}")820 print(f"subleading: {time.time() - t0:.2f} s")821822823def main():824 parser = argparse.ArgumentParser()825 parser.add_argument("verb", choices=["carry", "question", "table", "obstruction", "graph", "subleading", "all"])826 parser.add_argument("--span", type=int, default=16)827 parser.add_argument("--depth", type=int, default=12)828 parser.add_argument("--level", type=int, default=14)829 parser.add_argument("--modes", type=int, default=40)830 parser.add_argument("--cut", type=int, default=2000)831 parser.add_argument("--taylor", type=int, default=4)832 parser.add_argument("--jmax", type=int, default=24)833 parser.add_argument("--jlo", type=int, default=12)834 parser.add_argument("--per", type=int, default=16)835 args = parser.parse_args()836 widths = (2, 3, 4)837 if args.verb in ("carry", "all"):838 carry(args.span, args.depth)839 if args.verb in ("question", "all"):840 alphas = [(k, a) for k in widths for a in alphabets_at(k)]841 question(args.depth, args.level, alphas)842 if args.verb in ("table", "all"):843 table(widths, args.modes, args.cut, args.taylor)844 if args.verb in ("obstruction", "all"):845 obstruction(widths)846 if args.verb in ("graph", "all"):847 graph(widths, args.modes, args.cut, args.taylor, args.level)848 if args.verb in ("subleading", "all"):849 subleading(args.modes, args.jmax, args.jlo, args.per)850851852if __name__ == "__main__":853 main()