census.py
44.3 kB · python · 1003 lines
1import argparse2import cmath3import json4import math5import os6import time78LN2 = math.log(2.0)9PERIOD = 2.0 * math.pi / LN210ROUNDING = 1e-1411HERE = os.path.dirname(os.path.abspath(__file__))12CONTROL = os.path.join(HERE, "control.json")1314# ENGINE1516def roots(poly):17 n = len(poly) - 118 while n > 0 and abs(poly[n]) < 1e-14:19 n -= 120 c = [poly[i] / poly[n] for i in range(n + 1)]21 z = [(0.4 + 0.9j) ** i for i in range(n)]22 for _ in range(400):23 move = 0.024 for i in range(n):25 num = sum(c[j] * z[i] ** j for j in range(n + 1))26 den = 1.0 + 0j27 for j in range(n):28 if j != i:29 den *= z[i] - z[j]30 step = num / den31 z[i] -= step32 move = max(move, abs(step))33 if move < 1e-15:34 break35 return z363738class Rule:39 def __init__(self, width, code, peel):40 self.width, self.code, self.peel = width, code, peel41 self.size = 1 << (width - 1)42 mask = self.size - 143 self.T = [[0.0] * self.size for _ in range(self.size)]44 self.G1 = [[0.0] * self.size for _ in range(self.size)]45 for u in range(self.size):46 for a in (0, 1):47 w = (u << 1) | a48 if (code >> w) & 1:49 self.T[w & mask][u] += 1.050 self.G1[w & mask][u] += float(a)51 self.faddeev()52 self.degree = max(i for i in range(self.size + 1) if abs(self.det[i]) > 0.5)53 raw = [1.0 / x for x in roots(self.det[:self.degree + 1])] if self.degree else []54 self.eigen = [complex(z.real, 0.0) if abs(z.imag) < 1e-9 * abs(z) else z for z in raw]55 self.rho = max((abs(z) for z in self.eigen), default=0.0)56 self.alpha = math.log(self.rho) / LN2 if self.rho else float("-inf")57 self.low = [n for n in range(1, 1 << (peel - 1)) if self.accepts(n)]58 self.mid = [[] for _ in range(self.size)]59 for n in range(1 << (peel - 1), 1 << peel):60 if self.accepts(n):61 self.mid[n & mask].append(n)62 self.loglow = [math.log(n) for n in self.low]63 self.logmid = [[math.log(n) for n in row] for row in self.mid]64 self.counts = [0] * 6465 for j in range(1, 64):66 self.counts[j] = self.words(j)6768 def accepts(self, n):69 bits = bin(n)[2:]70 for i in range(len(bits) - self.width + 1):71 if not (self.code >> int(bits[i:i + self.width], 2)) & 1:72 return False73 return True7475 def words(self, j):76 if j < self.width:77 return 1 << (j - 1)78 live = [1.0 if u >> (self.width - 2) else 0.0 for u in range(self.size)]79 for _ in range(j - self.width + 1):80 live = [sum(self.T[v][u] * live[u] for u in range(self.size)) for v in range(self.size)]81 return sum(live)8283 def faddeev(self):84 n = self.size85 eye = [[1.0 if i == j else 0.0 for j in range(n)] for i in range(n)]86 m = [row[:] for row in eye]87 self.adj, c = [m], [1.0]88 for k in range(1, n + 1):89 tm = [[sum(self.T[i][t] * m[t][j] for t in range(n)) for j in range(n)] for i in range(n)]90 ck = -sum(tm[i][i] for i in range(n)) / k91 c.append(ck)92 m = [[tm[i][j] + (ck if i == j else 0.0) for j in range(n)] for i in range(n)]93 if k < n:94 self.adj.append(m)95 self.charpoly = list(reversed(c))96 self.det = c[:]9798 def deter(self, x):99 return sum(self.det[i] * x ** i for i in range(self.size + 1))100101 def apply_adj(self, x, v):102 out = [0j] * self.size103 p = 1.0 + 0j104 for k in range(self.size):105 for i in range(self.size):106 out[i] += p * sum(self.adj[k][i][j] * v[j] for j in range(self.size))107 p *= x108 return out109110 def adj_abs(self, x):111 a = abs(x)112 return [[sum(abs(self.adj[k][i][j]) * a ** k for k in range(self.size))113 for j in range(self.size)] for i in range(self.size)]114115 def tail(self, sigma):116 if sigma <= self.alpha + 1e-9:117 return float("inf")118 total = 0.0119 for j in range(self.peel, 64):120 total += self.counts[j] * math.exp(-(j - 1) * sigma * LN2)121 return total122123 def combs(self, m):124 out = []125 for z in self.eigen:126 re = math.log(abs(z)) / LN2 - m127 out.append((re, cmath.phase(z) / LN2))128 return out129130131def two_pow(w):132 return cmath.exp(-w * LN2)133134135def ladder(rule, s, deep, shift, cut):136 n = rule.size137 levels = max(1, int(math.ceil(shift - s.real)))138 rows = levels + cut + 2139 g = [[0j] * n for _ in range(rows)]140 bound = [[0.0] * n for _ in range(rows)]141 for j in range(levels, rows):142 t = rule.tail(s.real + j)143 bound[j] = [t] * n144 num, bnum = None, None145 for j in range(levels - 1, -1, -1):146 w = s + j147 sig, aw = w.real, abs(w)148 acc = [sum(cmath.exp(-w * L) for L in rule.logmid[u]) for u in range(n)]149 bacc = [0.0] * n150 c, cb = 1.0 + 0j, 1.0151 for l in range(1, cut + 1):152 c = c * (-w - (l - 1)) / l153 cb = cb * (aw + l - 1) / l154 weight = c * two_pow(w + l)155 wb = cb * math.exp(-(sig + l) * LN2)156 for i in range(n):157 acc[i] += weight * sum(rule.G1[i][t] * g[j + l][t] for t in range(n))158 bacc[i] += wb * sum(rule.G1[i][t] * bound[j + l][t] for t in range(n))159 ratio = ((aw + cut + 1) / (cut + 2)) * 2.0 ** (-rule.peel)160 if ratio >= 0.5:161 raise ValueError(f"the l-cut does not close at s = {s}, ratio {ratio}")162 rest = cb * (aw + cut) / (cut + 1) * math.exp(-(sig + cut + 1) * LN2)163 rest *= rule.tail(sig + cut + 1) / (1.0 - ratio)164 for i in range(n):165 bacc[i] += rest * sum(rule.G1[i][t] for t in range(n))166 if j == 0:167 num, bnum = acc, bacc168 if not deep:169 break170 x = two_pow(w)171 det = rule.deter(x)172 g[j] = [z / det for z in rule.apply_adj(x, acc)]173 ia = rule.adj_abs(x)174 bound[j] = [sum(ia[i][t] * bacc[t] for t in range(n)) / abs(det) for i in range(n)]175 return g[0], bound[0], num, bnum176177178class Engine:179 def __init__(self, rule, shift, cut, add=()):180 self.rule, self.shift, self.cut = rule, shift, cut181 self.add = tuple(sorted(add))182 self.logadd = [math.log(n) for n in self.add]183184 def poly_low(self, s):185 return (sum(cmath.exp(-s * L) for L in self.rule.loglow)186 + sum(cmath.exp(-s * L) for L in self.logadd))187188 def poly_abs(self, sigma):189 return (sum(math.exp(-sigma * L) for L in self.rule.loglow)190 + sum(math.exp(-sigma * L) for L in self.logadd))191192 def cofactor(self, s):193 _, _, num, bnum = ladder(self.rule, s, False, self.shift, self.cut)194 x = two_pow(s)195 det = self.rule.deter(x)196 top = sum(self.rule.apply_adj(x, num))197 ia = self.rule.adj_abs(x)198 n = self.rule.size199 lift = [sum(ia[i][j] for i in range(n)) for j in range(n)]200 bound = sum(lift[j] * bnum[j] for j in range(n))201 scale = abs(det) * self.poly_abs(s.real) + sum(lift[j] * abs(num[j]) for j in range(n))202 return det * self.poly_low(s) + top, bound + ROUNDING * scale203204 def zeta(self, s):205 g, b, _, _ = ladder(self.rule, s, True, self.shift, self.cut)206 value = self.poly_low(s) + sum(g)207 return value, sum(b) + ROUNDING * (self.poly_abs(s.real) + sum(abs(z) for z in g))208209 def residue(self, s0):210 value, bound = self.cofactor(s0)211 x = two_pow(s0)212 slope = sum(i * self.rule.det[i] * x ** (i - 1) for i in range(1, self.rule.size + 1))213 den = -x * LN2 * slope214 return value / den, bound / abs(den)215216# BOX217218def geometry(rule, left, right, height):219 combs, poles = [], []220 for re, off in rule.combs(0):221 pts = [complex(re, off + PERIOD * j) for j in range(-2, int(height / PERIOD) + 3)]222 combs.append((re, [p for p in pts if 0.0 < p.imag < height]))223 m = 1224 while re - m > left:225 poles += [complex(re - m, p.imag) for p in pts if 0.0 < p.imag < height]226 m += 1227 combs.sort(key=lambda c: -c[0])228 lines = sorted({round(c[0], 12) for c in combs} | {round(p.real, 12) for p in poles})229 lines = [v for v in lines if left < v < right]230 edges = [left] + [0.5 * (lines[i] + lines[i + 1]) for i in range(len(lines) - 1)] + [right]231 wide = []232 for i in range(len(edges) - 1):233 span = edges[i + 1] - edges[i]234 parts = max(1, int(math.ceil(span / 0.75)))235 wide += [edges[i] + span * t / parts for t in range(parts)]236 edges = wide + [right]237 heights = sorted({round(p.imag, 9) for c in combs for p in c[1]}238 | {round(p.imag, 9) for p in poles})239 rows = [0.02] + [0.5 * (heights[i] + heights[i + 1]) for i in range(len(heights) - 1)]240 rows = [r for r in rows if r < height] + [height]241 return combs, poles, edges, rows242243244def show(z, w=12):245 return f"{z.real:.{w}f}{z.imag:+.{w}f}i"246247248class Census:249 def __init__(self, engine, combs, poles, edges, rows, seed):250 self.engine, self.combs, self.poles = engine, combs, poles251 self.edges, self.rows, self.seed = edges, rows, seed252 self.cache, self.worst, self.bound, self.calls = {}, 0.0, 0.0, 0253 self.probes = 0254255 def value(self, s):256 key = (round(s.real, 12), round(s.imag, 12))257 hit = self.cache.get(key)258 if hit is None:259 hit, b = self.engine.cofactor(s)260 self.cache[key] = hit261 self.bound = max(self.bound, b)262 self.calls += 1263 return hit264265 def segment(self, a, b, fa, fb, depth):266 step = cmath.phase(fb / fa)267 if abs(step) <= 1.0 or depth >= 18:268 self.worst = max(self.worst, abs(step))269 return step270 m = 0.5 * (a + b)271 fm = self.value(m)272 return self.segment(a, m, fa, fm, depth + 1) + self.segment(m, b, fm, fb, depth + 1)273274 def side(self, a, b):275 steps = max(4, int(math.ceil(abs(b - a) / self.seed)))276 total, prev, fprev = 0.0, a, self.value(a)277 for i in range(1, steps + 1):278 nxt = a + (b - a) * i / steps279 fnext = self.value(nxt)280 total += self.segment(prev, nxt, fprev, fnext, 0)281 prev, fprev = nxt, fnext282 return total283284 def winding(self, box):285 x0, x1, y0, y1 = box286 corner = [complex(x0, y0), complex(x1, y0), complex(x1, y1), complex(x0, y1)]287 return sum(self.side(corner[i], corner[(i + 1) % 4]) for i in range(4)) / (2.0 * math.pi)288289 def count(self, box):290 x0, x1, y0, y1 = box291 turn = self.winding(box)292 held = sum(1 for p in self.poles if x0 < p.real < x1 and y0 < p.imag < y1)293 return int(round(turn + held)), turn, held294295 def split(self, box):296 x0, x1, y0, y1 = box297 marks = self.poles + [p for c in self.combs for p in c[1]]298 if x1 - x0 >= y1 - y0:299 cut = 0.5 * (x0 + x1)300 while any(abs(p.real - cut) < 1e-4 for p in marks):301 cut += 1.7e-3302 return (x0, cut, y0, y1), (cut, x1, y0, y1)303 cut = 0.5 * (y0 + y1)304 while any(abs(p.imag - cut) < 1e-4 for p in marks):305 cut += 1.7e-3306 return (x0, x1, y0, cut), (x0, x1, cut, y1)307308 def hunt(self, box, n, depth):309 if n <= 0:310 return []311 x0, x1, y0, y1 = box312 if (x1 - x0 < 0.01 and y1 - y0 < 0.01) or depth > 44:313 return [complex(0.5 * (x0 + x1), 0.5 * (y0 + y1))] * n314 left, right = self.split(box)315 nl = self.count(left)[0]316 return self.hunt(left, nl, depth + 1) + self.hunt(right, n - nl, depth + 1)317318 def polish(self, z):319 a, b = z, z + 1e-4320 fa, fb = self.value(a), self.value(b)321 for _ in range(60):322 if abs(fb - fa) < 1e-300:323 break324 c = b - fb * (b - a) / (fb - fa)325 if abs(c - b) < 1e-13:326 b = c327 break328 a, fa, b, fb = b, fb, c, self.value(c)329 return b, abs(self.value(b))330331 def circle(self, s0, eps, n=48):332 total, bound = 0j, 0.0333 for i in range(n):334 u = cmath.exp(2j * math.pi * i / n)335 v, b = self.engine.cofactor(s0 + eps * u)336 total += v * u337 bound = max(bound, b)338 self.probes += 1339 return eps * total / n, eps * bound340341 def arc(self, c, rad, a, b, fa, fb, depth):342 step = cmath.phase(fb / fa)343 if abs(step) <= 1.0 or depth >= 18:344 self.worst = max(self.worst, abs(step))345 return step346 m = 0.5 * (a + b)347 fm = self.value(c + rad * cmath.exp(1j * m))348 return self.arc(c, rad, a, m, fa, fm, depth + 1) + self.arc(c, rad, m, b, fm, fb, depth + 1)349350 def ring(self, c, rad, n=40):351 f0 = self.value(c + rad)352 total, prev, fprev = 0.0, 0.0, f0353 for i in range(1, n + 1):354 th = 2.0 * math.pi * i / n355 fnext = f0 if i == n else self.value(c + rad * cmath.exp(1j * th))356 total += self.arc(c, rad, prev, th, fprev, fnext, 0)357 prev, fprev = th, fnext358 return total / (2.0 * math.pi)359360 def disc(self, c, rad):361 held = sum(1 for p in self.poles if abs(p - c) < rad)362 return int(round(self.ring(c, rad))) + held363364 def guarded(self, c, rad, guard):365 near = min((abs(abs(p - c) - rad) for p in self.poles), default=float("inf"))366 return self.disc(c, rad - guard), self.disc(c, rad + guard), near367368 def regular(self, s0, eps, n=64):369 total = 0j370 for i in range(n):371 u = eps * cmath.exp(2j * math.pi * i / n)372 total += self.value(s0 + u) / self.engine.rule.deter(two_pow(s0 + u))373 return total / n374375# VERBS376377def census(args):378 start = time.perf_counter()379 rule = Rule(args.width, args.code, args.peel)380 engine = Engine(rule, args.shift, args.cut)381 combs, poles, edges, rows = geometry(rule, args.left, args.right, args.height)382 height = rows[-1]383 run = Census(engine, combs, poles, edges, rows, args.seed)384 print(f"code {args.code}, k = {args.width}, D = 1, base 2, states {rule.size},"385 f" peel {args.peel}, shift {args.shift}, cut {args.cut}")386 print(f"det(I - x T) coefficients {[round(c, 12) for c in rule.det]}")387 print(f"eigenvalues {' '.join(show(z, 12) for z in rule.eigen)}, alpha {rule.alpha:.12f}")388 for i, (re, pts) in enumerate(combs):389 print(f"comb {i} Re s = {re:.12f} teeth {len(pts)} first"390 f" {pts[0].imag if pts else float('nan'):.12f} spacing {PERIOD:.12f}")391 print(f"box Re s in [{args.left}, {args.right}], Im s in [{rows[0]}, {height:.12f}]")392 print(f"cofactor poles in the box {len(poles)}"393 f" on Re s = {sorted({round(p.real, 9) for p in poles})}")394 ranked = sorted(poles, key=lambda q: q.imag)395 for p in ranked:396 r1, b1 = run.circle(p, args.ring)397 r2, _ = run.circle(p, 0.4 * args.ring)398 lift, _ = engine.cofactor(p + 1e-5)399 print(f" pole Im {p.imag:10.6f} residue {show(r1, 12)} bound {b1:.2e},"400 f" at radius {0.4 * args.ring:.3f} gap {abs(r1 - r2):.2e},"401 f" simple to {abs(1e-5 * lift - r1) / abs(r1):.2e}")402 if len(ranked) > 1:403 blank = complex(ranked[0].real, 0.5 * (ranked[0].imag + ranked[1].imag))404 rb, _ = run.circle(blank, args.ring)405 print(f" blank point on the same line Im {blank.imag:10.6f}"406 f" circle mean {abs(rb):.2e}")407 cells, total = [], 0408 for i in range(len(edges) - 1):409 for j in range(len(rows) - 1):410 box = (edges[i], edges[i + 1], rows[j], rows[j + 1])411 n, turn, held = run.count(box)412 cells.append((box, n))413 total += n414 if n or held:415 print(f"cell Re [{box[0]:7.4f},{box[1]:7.4f}] Im [{box[2]:9.5f},{box[3]:9.5f}]"416 f" winding {turn:+9.5f} poles {held} zeros {n}")417 print(f"zeros in the box {total}, cells {len(cells)}")418 zeros = []419 for box, n in cells:420 for z in run.hunt(box, n, 0):421 zeros.append(run.polish(z))422 zeros.sort(key=lambda p: p[0].imag)423 print(f"located {len(zeros)} of {total},"424 f" largest residual {max((r for _, r in zeros), default=0.0):.3e}")425 for z, _ in zeros:426 gaps = " ".join(f"c{i} {min((abs(z - t) for t in pts), default=float('inf')):.9f}"427 for i, (_, pts) in enumerate(combs))428 print(f" zero {show(z)} {gaps}")429 on_edge = [z for z, _ in zeros430 if min(abs(z.real - args.left), abs(z.real - args.right),431 abs(z.imag - rows[0]), abs(z.imag - height)) < 0.02]432 on_pole = [p for p in poles433 if min(abs(p.real - args.left), abs(p.real - args.right),434 abs(p.imag - rows[0]), abs(p.imag - height)) < 0.02]435 print(f"zeros within 0.02 of a box edge {len(on_edge)} {[show(z, 9) for z in on_edge]}")436 print(f"poles within 0.02 of a box edge {len(on_pole)}")437 inner = [min([abs(z.real - e) for e in edges[1:-1]] + [abs(z.imag - r) for r in rows[1:-1]]438 + [float("inf")]) for z, _ in zeros]439 tight = min(range(len(zeros)), key=lambda i: inner[i])440 print(f"tightest clearance to an internal cell edge {inner[tight]:.6f}"441 f" at {show(zeros[tight][0])}")442 held = set()443 for i, (re, pts) in enumerate(combs):444 near = [min(abs(z - t) for z, _ in zeros) for t in pts]445 for t in pts:446 for z, _ in zeros:447 if abs(z - t) < args.rho:448 held.add((round(z.real, 9), round(z.imag, 9)))449 print(f"comb {i} Re s = {re:.12f} teeth {len(pts)}"450 f" occupied at radius {args.rho} {sum(1 for d in near if d < args.rho)}")451 print(f" distances {' '.join(f'{d:.9f}' for d in near)}")452 for t in pts:453 r, rb = engine.residue(t)454 reg = run.regular(t, args.eps)455 u1 = -r / reg456 best = min(zeros, key=lambda p: abs(p[0] - t))457 disc = sum(1 for z, _ in zeros if abs(z - t) < args.rho)458 print(f" tooth Im {t.imag:10.6f} r {show(r, 9)} bound {rb:.2e}"459 f" R {show(reg, 9)} u1 {show(u1, 9)} abs {abs(u1):.9f}"460 f" zero at {abs(best[0] - t):.9f}"461 f" modulus miss {abs(abs(best[0] - t) - abs(u1)):.9f}"462 f" vector miss {abs(best[0] - t - u1):.9f}"463 f" zeros in the disc {disc}")464 family = [z for z, _ in zeros if (round(z.real, 9), round(z.imag, 9)) not in held]465 print(f"tooth zeros {len(zeros) - len(family)}, second family {len(family)}")466 if family:467 print(f"second family Re s in [{min(z.real for z in family):.12f},"468 f" {max(z.real for z in family):.12f}],"469 f" mean {sum(z.real for z in family) / len(family):.12f}")470 for i, (re, _) in enumerate(combs):471 hug = [z for z, _ in zeros if abs(z.real - re) < args.hug]472 print(f"zeros within {args.hug} of comb {i}'s line: {len(hug)}"473 f" {[f'{show(z, 9)} at {abs(z.real - re):.9f}' for z in hug]}")474 axis = 0.5 * (combs[0][0] + combs[-1][0])475 mates = sum(1 for i, (z, _) in enumerate(zeros) for j, (w, _) in enumerate(zeros) if j != i476 and abs(w.imag - z.imag) < 0.05 and abs(w.real + z.real - 2.0 * axis) < 0.05)477 selfish = [z for z, _ in zeros if abs(2.0 * (z.real - axis)) < 0.05]478 near = min(abs(2.0 * (z.real - axis)) for z, _ in zeros)479 print(f"reflection partners other than the zero itself about Re s = {axis:.12f},"480 f" the axis of the comb lines: {mates} of {len(zeros)}")481 print(f"zeros whose own reflection lands within 0.05 of them {len(selfish)},"482 f" nearest self-reflection {near:.12f} {[show(z, 12) for z in selfish]}")483 guard, gbound = engine.zeta(complex(args.right, 0.0))484 print(f"zeta_W({args.right}) = {guard.real:.12f} with a_min = 1, bound {gbound:.3e}")485 print(f"cuts: Im s in [0, {rows[0]}] uncounted, Im s > {height:.12f} uncounted,"486 f" Re s < {args.left} uncounted")487 print(f"largest surviving phase step {run.worst:.6f} rad, largest bound {run.bound:.3e},"488 f" evaluations {run.calls} on the contour and the hunt,"489 f" {run.probes} off it on the residue circles, runtime {time.perf_counter() - start:.2f} s")490491492# CLASSES493494def window_maps(k):495 tables = set()496 for flip in (0, 1):497 for rev in (False, True):498 t = []499 for w in range(1 << k):500 ds = [((w >> (k - 1 - j)) & 1) ^ flip for j in range(k)]501 if rev:502 ds = ds[::-1]503 v = 0504 for c in ds:505 v = (v << 1) | c506 t.append(v)507 tables.add(tuple(t))508 return sorted(tables)509510511def classes(k):512 nbits = 1 << k513 maps = []514 for t in window_maps(k):515 arr = [0] * (1 << nbits)516 for code in range(1, 1 << nbits):517 low = code & -code518 arr[code] = arr[code ^ low] | (1 << t[low.bit_length() - 1])519 maps.append(arr)520 seen = bytearray(1 << nbits)521 out = []522 for c in range(1 << nbits):523 if seen[c]:524 continue525 orb = {m[c] for m in maps}526 for x in orb:527 seen[x] = 1528 out.append((c, len(orb)))529 return out530531532def spectrum(rule):533 order = sorted(rule.eigen, key=lambda z: -abs(z))534 mods = []535 for z in order:536 if not any(abs(abs(z) - m) < 1e-9 for m in mods):537 mods.append(abs(z))538 gap = min((abs(order[i] - order[j]) for i in range(len(order))539 for j in range(i + 1, len(order))), default=float("inf"))540 return order, mods, gap541542# TEETH543544def flank(rule, walls):545 left, right, low, high = walls546 out = []547 for re, off in rule.combs(0):548 pts = [off + PERIOD * j for j in range(-2, int(high / PERIOD) + 3)]549 m = 1550 while re - m > left - 1.0:551 out += [complex(re - m, y) for y in pts if low - 1.0 < y < high + 1.0]552 m += 1553 return out554555556def wall(p, walls):557 left, right, low, high = walls558 dx = max(left - p.real, p.real - right, 0.0)559 dy = max(low - p.imag, p.imag - high, 0.0)560 if dx == 0.0 and dy == 0.0:561 return min(p.real - left, right - p.real, p.imag - low, high - p.imag)562 return math.hypot(dx, dy)563564565def sweep(engine, rule, args):566 combs, poles, edges, rows = geometry(rule, args.left, args.right, args.height)567 run = Census(engine, combs, poles, edges, rows, args.seed)568 total = 0569 cells = []570 for i in range(len(edges) - 1):571 for j in range(len(rows) - 1):572 box = (edges[i], edges[i + 1], rows[j], rows[j + 1])573 n = run.count(box)[0]574 cells.append((box, n))575 total += n576 zeros = []577 for box, n in cells:578 for z in run.hunt(box, n, 0):579 zeros.append(run.polish(z))580 zeros.sort(key=lambda p: p[0].imag)581 return combs, poles, edges, rows, cells, total, zeros, run582583584def teeth(args):585 start = time.perf_counter()586 reps = classes(args.width)587 print(f"D = 1, base 2, width {args.width}, classes under G_(1,k) {len(reps)},"588 f" peel {args.peel}, shift {args.shift}, cut {args.cut}")589 print(f"box Re s in [{args.left}, {args.right}], Im s in [0.02, {args.height}],"590 f" contour seed {args.seed}, occupancy radius {args.rho}")591 picked, skipped = [], []592 for code, size in reps:593 rule = Rule(args.width, code, args.peel)594 if not rule.eigen:595 continue596 order, mods, gap = spectrum(rule)597 if len(mods) != 2:598 continue599 if gap < 1e-9:600 skipped.append(code)601 continue602 picked.append((code, size, rule, order, mods))603 picked.sort(key=lambda r: -r[4][1] / r[4][0])604 print(f"two-line classes {len(picked)} of {len(reps)},"605 f" dropped for a repeated eigenvalue {len(skipped)} {skipped}")606 table, law, guard = [], [], []607 for code, size, rule, order, mods in picked:608 engine = Engine(rule, args.shift, args.cut)609 combs, poles, edges, rows, cells, total, zeros, run = sweep(engine, rule, args)610 height = rows[-1]611 print(f"code {code}, k = {args.width}, class size {size}, states {rule.size},"612 f" det {[round(c, 9) for c in rule.det]}")613 print(f" eigenvalues {' '.join(show(z, 9) for z in order)}, rho {mods[0]:.12f},"614 f" abs(l2)/rho {mods[1] / mods[0]:.12f}")615 print(f" cells {len(cells)} poles {len(poles)} zeros {total} located {len(zeros)}"616 f" largest residual {max((r for _, r in zeros), default=0.0):.3e}"617 f" phase step {run.worst:.6f} bound {run.bound:.3e} calls {run.calls}")618 held, row = set(), []619 for j, m in enumerate(mods):620 pts = sorted([p for re, ps in combs if abs(2.0 ** re - m) < 1e-9 for p in ps],621 key=lambda p: p.imag)622 args_of = sorted({round(cmath.phase(z), 9) for z in order if abs(abs(z) - m) < 1e-9})623 near = [min((abs(z - t) for z, _ in zeros), default=float("inf")) for t in pts]624 for t, d in zip(pts, near):625 if d < args.rho:626 for z, _ in zeros:627 if abs(z - t) < args.rho:628 held.add((round(z.real, 9), round(z.imag, 9)))629 hit = sum(1 for d in near if d < args.rho)630 print(f" line {j} Re s = {math.log(m) / LN2:.12f} arg"631 f" {' '.join(f'{a:+.9f}' for a in args_of)} teeth {len(pts)}"632 f" occupied {hit} least distance {min(near, default=float('nan')):.9f}")633 print(f" distances {' '.join(f'{d:.9f}' for d in near)}")634 for t, d in zip(pts, near):635 r, rb = engine.residue(t)636 reg = run.regular(t, args.eps)637 u1 = -r / reg638 best = min(zeros, key=lambda p: abs(p[0] - t), default=None)639 vec = abs(best[0] - t - u1) if best else float("inf")640 disc = sum(1 for z, _ in zeros if abs(z - t) < args.rho)641 law.append((code, j, abs(u1), d, d < args.rho, vec, disc))642 print(f" tooth Im {t.imag:10.6f} r {show(r, 9)} bound {rb:.2e}"643 f" R {show(reg, 9)} u1 {show(u1, 9)} abs {abs(u1):.9f}"644 f" zero at {d:.9f} modulus miss {abs(d - abs(u1)):.9f}"645 f" vector miss {vec:.9f} zeros in the disc {disc}")646 row.append((len(pts), hit, args_of))647 off = [z for z, _ in zeros if (round(z.real, 9), round(z.imag, 9)) not in held]648 edge = [z for z, _ in zeros649 if min(abs(z.real - args.left), abs(z.real - args.right),650 abs(z.imag - rows[0]), abs(z.imag - height)) < 0.02]651 walls = (args.left, args.right, rows[0], height)652 reach = [wall(p, walls) for p in flank(rule, walls)]653 on_pole = [d for d in reach if d < 0.02]654 clear = min(reach + [float("inf")])655 guard.append((code, len(edge), len(on_pole), clear))656 print(f" zeros off every tooth {len(off)} of {total},"657 f" within 0.02 of a contour {len(edge)} {[show(z, 9) for z in edge]},"658 f" poles within 0.02 of a contour {len(on_pole)},"659 f" least pole clearance {clear:.12f}")660 reach = mods[1] and math.log(mods[1]) / LN2 - args.rho661 print(f" second line disc reaches Re s = {reach:.12f},"662 f" inside the box {reach > args.left}")663 table.append((code, size, mods[1] / mods[0], row, total, len(off)))664 print("table rule class ratio arg2 teeth1 occupied1 teeth2 occupied2 zeros off")665 for code, size, ratio, row, total, off in table:666 a2 = " ".join(f"{a:+.6f}" for a in row[1][2])667 print(f" {code:4d} {size:3d} {ratio:.9f} [{a2}]"668 f" {row[0][0]:3d} {row[0][1]:3d} {row[1][0]:3d} {row[1][1]:3d} {total:3d} {off:3d}")669 full = [r for r in table if r[3][1][0] and r[3][1][1] == r[3][1][0]]670 empty = [r for r in table if r[3][1][0] and r[3][1][1] == 0]671 print(f"second line full {len(full)} {[r[0] for r in full]},"672 f" empty {len(empty)} {[r[0] for r in empty]},"673 f" partial {len(table) - len(full) - len(empty)}")674 inside = [t for t in law if t[2] < args.eps]675 outside = [t for t in law if t[2] >= args.rho]676 print(f"teeth {len(law)}, occupied {sum(1 for t in law if t[4])};"677 f" first-order prediction abs(u1) below {args.eps}, inside the disc that builds R:"678 f" {len(inside)} teeth, occupied {sum(1 for t in inside if t[4])}")679 print(f"prediction abs(u1) at or above {args.rho}: {len(outside)} teeth,"680 f" occupied {sum(1 for t in outside if t[4])}")681 band = [t for t in law if args.eps <= t[2] < args.rho]682 print(f"prediction in [{args.eps}, {args.rho}): {len(band)} teeth,"683 f" occupied {sum(1 for t in band if t[4])} {[(t[0], round(t[2], 6), round(t[3], 6)) for t in band]}")684 worst = max(law, key=lambda t: abs(t[3] - t[2]), default=None)685 print(f"largest modulus miss of the first-order law {abs(worst[3] - worst[2]):.9f}"686 f" at code {worst[0]} line {worst[1]}, abs(u1) {worst[2]:.9f} zero at {worst[3]:.9f}")687 vworst = max(law, key=lambda t: t[5], default=None)688 print(f"largest vector miss of the first-order law {vworst[5]:.9f}"689 f" at code {vworst[0]} line {vworst[1]}, abs(u1) {vworst[2]:.9f} zero at {vworst[3]:.9f}")690 doubles = [t for t in law if t[6] > 1]691 print(f"teeth holding more than one zero inside {args.rho} {len(doubles)}"692 f" {[(t[0], t[1], t[6]) for t in doubles]}, zeros inside a tooth disc"693 f" {sum(t[6] for t in law)} against occupied teeth"694 f" {sum(1 for t in law if t[4])}")695 tight = min(guard, key=lambda g: g[3], default=None)696 print(f"contour guard: zeros within 0.02 of a contour {sum(g[1] for g in guard)},"697 f" poles within 0.02 of a contour {sum(g[2] for g in guard)},"698 f" least pole clearance {tight[3]:.12f} at code {tight[0]}")699 print(f"runtime {time.perf_counter() - start:.2f} s")700701702# BRIDGE703704def bridge(args):705 start = time.perf_counter()706 a, b = Rule(2, 7, args.peel), Rule(3, 55, args.peel)707 top = 1 << args.span708 extra = [n for n in range(1, top) if b.accepts(n) and not a.accepts(n)]709 lost = [n for n in range(1, top) if a.accepts(n) and not b.accepts(n)]710 print(f"code 7 at k = 2 against code 55 at k = 3 over 1 .. {top - 1}:"711 f" in 55 and not 7 {extra}, in 7 and not 55 {lost}")712 print(f"det(I - x T) code 7 {[round(c, 12) for c in a.det]},"713 f" code 55 {[round(c, 12) for c in b.det]}, states {a.size} against {b.size}")714 ea, eb = Engine(a, args.shift, args.cut), Engine(b, args.shift, args.cut)715 worst, over = 0.0, 0716 alpha = math.log((1.0 + math.sqrt(5.0)) / 2.0) / LN2717 for s in [complex(2.0, 0.0), complex(0.8, 0.0), complex(-0.95, 20.0),718 complex(alpha, PERIOD), complex(-alpha, 0.5 * PERIOD),719 complex(-alpha, 1.5 * PERIOD), complex(-0.442302243578, 4.612546440182)]:720 va, ba = ea.cofactor(s)721 vb, bb = eb.cofactor(s)722 want = b.deter(two_pow(s)) * cmath.exp(-s * math.log(3.0))723 gap = abs(vb - va - want)724 worst = max(worst, gap)725 over += gap > ba + bb726 print(f" s {show(s, 6)} Z7 {show(va, 12)} Z55 {show(vb, 12)}"727 f" gap {gap:.3e} bound {ba + bb:.3e} {'in' if gap <= ba + bb else 'OUT'}")728 print(f"largest gap {worst:.3e}, outside its bound {over} of 7,"729 f" runtime {time.perf_counter() - start:.2f} s")730731732def control(args):733 rule = Rule(2, 7, args.peel)734 engine = Engine(rule, args.shift, args.cut)735 data = json.load(open(CONTROL))736 print(f"control {data['source']}, dps {data['dps']}, peel {data['peel']},"737 f" shift {data['shift']}, cut {data['cut']}")738 worst, over = 0.0, 0739 for row in data["rows"]:740 s = complex(float(row["re"]), float(row["im"]))741 got, bound = (engine.cofactor(s) if row["kind"] == "Z" else742 engine.zeta(s) if row["kind"] == "zeta" else engine.residue(s))743 want = complex(float(row["value"][0]), float(row["value"][1]))744 gap = abs(got - want)745 worst = max(worst, gap)746 over += gap > bound747 print(f" {row['kind']:8s} {show(s, 6)} value {show(got, 12)} gap {gap:.3e}"748 f" bound {bound:.3e} {'in' if gap <= bound else 'OUT'}")749 print(f"largest gap {worst:.3e}, outside its bound {over} of {len(data['rows'])}")750751752# DIAL753754def lines_of(combs):755 out = []756 for re, pts in combs:757 key = round(re, 9)758 for i, (r, ps) in enumerate(out):759 if abs(r - key) < 1e-9:760 out[i] = (r, ps + list(pts))761 break762 else:763 out.append((key, list(pts)))764 out.sort(key=lambda c: -c[0])765 return [(r, sorted(ps, key=lambda p: p.imag)) for r, ps in out]766767768def exact_min(reg, t, cand, cap):769 import numpy as np770 half = len(cand) // 2771 left, right = cand[:half], cand[half:]772773 def table(part):774 sums = np.array([0.0 + 0.0j])775 size = np.array([0])776 for n in part:777 v = cmath.exp(-t * math.log(n))778 sums = np.concatenate([sums, sums + v])779 size = np.concatenate([size, size + 1])780 return sums, size781782 ls, lk = table(left)783 rs, rk = table(right)784 best, pick = float("inf"), None785 for j in range(ls.size):786 room = cap - int(lk[j])787 if room < 0:788 continue789 ok = np.flatnonzero(rk <= room)790 d = np.abs(reg + ls[j] + rs[ok])791 w = int(np.argmin(d))792 if float(d[w]) < best:793 best = float(d[w])794 idx = int(ok[w])795 pick = ([left[e] for e in range(len(left)) if j >> e & 1]796 + [right[e] for e in range(len(right)) if idx >> e & 1])797 return best, sorted(pick)798799800def dial(args):801 start = time.perf_counter()802 rule = Rule(args.width, args.code, args.peel)803 base = Engine(rule, args.shift, args.cut)804 combs, poles, edges, rows = geometry(rule, args.left, args.right, args.height)805 lines = lines_of(combs)806 height = rows[-1]807 cand = [n for n in range(2, args.top + 1) if not rule.accepts(n)]808 print(f"code {args.code}, k = {args.width}, D = 1, base 2, states {rule.size},"809 f" peel {args.peel}, shift {args.shift}, cut {args.cut}")810 print(f"det(I - x T) coefficients {[round(c, 12) for c in rule.det]}")811 print(f"eigenvalues {' '.join(show(z, 12) for z in rule.eigen)}, alpha {rule.alpha:.12f}")812 print(f"box Im s in [{rows[0]}, {height:.12f}], Re s in [{args.left}, {args.right}],"813 f" occupancy radius {args.rho}, ring samples 40, edge guard {args.guard},"814 f" R radius {args.eps}")815 for i, (re, pts) in enumerate(lines):816 print(f"line {i} Re s = {re:.12f} teeth {len(pts)}"817 f" Im {' '.join(f'{p.imag:.6f}' for p in pts)}")818 walls = (args.left, args.right, rows[0], height)819 wide = flank(rule, walls)820 reach = min(t.real for _, pts in lines for t in pts) - args.rho - args.guard821 outer = max([p.real for p in wide if p.real <= args.left] + [float("-inf")])822 print(f"cofactor poles in the box {len(poles)}"823 f" on Re s = {sorted({round(p.real, 9) for p in poles})}, padded pole list"824 f" {len(wide)} on Re s = {sorted({round(p.real, 9) for p in wide})}")825 print(f"deepest disc reach Re s = {reach:.12f}, nearest pole line left of the box"826 f" {outer:.12f}, clear {reach > outer}")827 print(f"candidates, the integers of 2 .. {args.top} outside S_W, {len(cand)}: {cand}")828 ref = Census(base, combs, wide, edges, rows, args.seed)829 print("falsification, the perturbation is entire so every pole and residue must hold")830 worst_res, worst_val, worst_off, worst_add = 0.0, 0.0, 0.0, 0.0831 for probe in ([cand[0]], cand[:3], [cand[-1]]):832 eng = Engine(rule, args.shift, args.cut, probe)833 run = Census(eng, combs, wide, edges, rows, args.seed)834 for p in sorted(poles, key=lambda q: q.imag):835 r0, _ = ref.circle(p, args.ring)836 r1, _ = run.circle(p, args.ring)837 worst_res = max(worst_res, abs(r1 - r0))838 gaps = []839 for _, pts in lines:840 for t in pts:841 v0, _ = base.cofactor(t)842 v1, _ = eng.cofactor(t)843 gaps.append(abs(v1 - v0))844 worst_val = max(worst_val, abs(v1 - v0))845 off, add = 0.0, 0.0846 for _, pts in lines:847 for t in pts:848 s0 = t + args.rho849 part = rule.deter(two_pow(s0)) * sum(850 cmath.exp(-s0 * math.log(n)) for n in probe)851 v0, _ = base.cofactor(s0)852 v1, _ = eng.cofactor(s0)853 off = max(off, abs(v1 - v0 - part))854 add = max(add, abs(part))855 worst_off, worst_add = max(worst_off, off), max(worst_add, add)856 print(f" F {probe} poles {len(poles)} largest residue gap"857 f" {worst_res:.3e}, teeth {len(gaps)} largest Z gap {max(gaps):.3e},"858 f" off-tooth identity largest miss {off:.3e} against an added part of"859 f" {add:.6f}")860 print(f"largest residue gap over every probe {worst_res:.3e},"861 f" largest tooth value gap {worst_val:.3e}, largest off-tooth identity miss"862 f" {worst_off:.3e} against an added part of up to {worst_add:.6f}")863 print("the residue and tooth probes are blind to the size of an entire addition,"864 " the circle mean annihilating it and the determinant vanishing at a tooth;"865 " the off-tooth identity is the one that measures it")866 rows_out = []867 for i, (re, pts) in enumerate(lines):868 for t in pts:869 r, rb = base.residue(t)870 reg = ref.regular(t, args.eps)871 u1 = -r / reg872 lo, hi, near = ref.guarded(t, args.rho, args.guard)873 rows_out.append((i, t, r, reg, u1, lo, hi, near))874 print(f" tooth line {i} Im {t.imag:10.6f} r {show(r, 9)} bound {rb:.2e}"875 f" R {show(reg, 9)} u1 {show(u1, 9)} abs {abs(u1):.9f}"876 f" zeros in the disc {lo}/{hi} pole clearance {near:.6f}")877 print(f"baseline occupied {sum(1 for q in rows_out if q[5] > 0)} of {len(rows_out)} teeth,"878 f" edge guard splits {sum(1 for q in rows_out if q[5] != q[6])}")879 grid, split = {}, 0880 for n in cand:881 eng = Engine(rule, args.shift, args.cut, [n])882 run = Census(eng, combs, wide, edges, rows, args.seed)883 for i, t, r, reg, u1, lo, hi, near in rows_out:884 a, b, _ = run.guarded(t, args.rho, args.guard)885 split += a != b886 grid[(round(t.imag, 6), n)] = (a, b)887 agree, cells, flips = 0, 0, []888 perline = {}889 print("grid, one row per tooth, one column per candidate in the printed order")890 for i, t, r, reg, u1, lo, hi, near in rows_out:891 meas, pred, bad = "", "", []892 for n in cand:893 a, b = grid[(round(t.imag, 6), n)]894 p = abs(-r / (reg + cmath.exp(-t * math.log(n)))) < args.rho895 meas += "1" if a else "0"896 pred += "1" if p else "0"897 cells += 1898 hit, tot, one, und, hitb, oneb = perline.get(i, (0, 0, 0, 0, 0, 0))899 perline[i] = (hit + (bool(a) == p), tot + 1, one + bool(a),900 und + (not a and b > 0), hitb + (bool(b) == p),901 oneb + bool(b))902 if bool(a) == p:903 agree += 1904 else:905 bad.append(n)906 got = [n for n, c in zip(cand, meas) if c == "1"]907 miss = [n for n, c in zip(cand, meas) if c == "0"]908 print(f" tooth line {i} Im {t.imag:10.6f} base {1 if lo else 0}"909 f" measured {meas} occupied {len(got)} of {len(cand)}")910 print(f" predicted {pred} mismatch {bad}")911 seam = [n for n in cand if grid[(round(t.imag, 6), n)][0]912 != grid[(round(t.imag, 6), n)][1]]913 print(f" last candidate reading empty {max(miss) if miss else None},"914 f" empty candidates {miss}")915 print(f" candidates whose disc holds a zero within {args.guard} of the"916 f" occupancy circle {seam}")917 out_got = [n for n in cand if grid[(round(t.imag, 6), n)][1]]918 out_miss = [n for n in cand if not grid[(round(t.imag, 6), n)][1]]919 print(f" outer reading occupied {len(out_got)} of {len(cand)},"920 f" smallest singleton {out_got[0] if out_got else None},"921 f" last candidate reading empty"922 f" {max(out_miss) if out_miss else None}")923 first = got[0] if got else None924 if first is None:925 print(f" no singleton at or below {args.top} occupies it")926 else:927 ph = (t.imag * math.log(first) / (2.0 * math.pi)) % 1.0928 print(f" smallest singleton {{{first}}} phase {ph:.6f}"929 f" predicted abs(u1) {abs(-r / (reg + cmath.exp(-t * math.log(first)))):.9f}"930 f" on the seam {first in seam}")931 if lo and len(got) < len(cand):932 flips.append((round(t.imag, 6), miss))933 print(f"grid cells {cells}, first-order law agrees {agree},"934 f" disagrees {cells - agree}, edge guard splits {split}")935 for i in sorted(perline):936 hit, tot, one, und, hitb, oneb = perline[i]937 print(f" line {i} cells {tot} measured occupied {one},"938 f" first-order law agrees {hit}, the constant occupied predictor agrees {one}")939 print(f" occupancy undetermined {und}, outer reading occupied {oneb},"940 f" first-order law agrees {hitb},"941 f" the constant occupied predictor agrees {oneb}")942 print(f"occupied teeth emptied by a singleton {len(flips)} {flips}")943 print(f"greedy against the regular part, the F that drives R + P_F to zero")944 for i, t, r, reg, u1, lo, hi, near in rows_out:945 if not lo:946 continue947 pool, pick, done = list(cand), [], None948 for _ in range(args.deep):949 best = min(pool, key=lambda n: abs(reg + sum(950 cmath.exp(-t * math.log(m)) for m in pick + [n])))951 pick.append(best)952 pool.remove(best)953 eng = Engine(rule, args.shift, args.cut, pick)954 run = Census(eng, combs, wide, edges, rows, args.seed)955 a, b, _ = run.guarded(t, args.rho, args.guard)956 pv = abs(reg + sum(cmath.exp(-t * math.log(m)) for m in pick))957 print(f" tooth line {i} Im {t.imag:10.6f} F {pick} abs(R + P_F) {pv:.9f}"958 f" predicted abs(u1) {abs(r) / pv:.9f} zeros in the disc {a}/{b}")959 if a == 0 and b == 0:960 done = list(pick)961 break962 print(f" tooth line {i} Im {t.imag:10.6f} emptied by {done}"963 f" size {len(done) if done else None}")964 if args.deep and len(cand) <= 24:965 floor, best = exact_min(reg, t, cand, args.deep)966 eng = Engine(rule, args.shift, args.cut, best)967 run = Census(eng, combs, wide, edges, rows, args.seed)968 a, b, _ = run.guarded(t, args.rho, args.guard)969 print(f" tooth line {i} Im {t.imag:10.6f} exact minimiser over the"970 f" {sum(math.comb(len(cand), e) for e in range(args.deep + 1))}"971 f" subsets of size at most {args.deep} {best}"972 f" abs(R + P_F) {floor:.9f} predicted abs(u1)"973 f" {abs(r) / floor:.9f} emptying threshold abs(r)/rho"974 f" {abs(r) / args.rho:.9f} zeros in the disc {a}/{b}")975 print(f"runtime {time.perf_counter() - start:.2f} s")976977978def main():979 parser = argparse.ArgumentParser()980 parser.add_argument("--width", type=int, default=2)981 parser.add_argument("--code", type=int, default=7)982 parser.add_argument("--peel", type=int, default=8)983 parser.add_argument("--shift", type=float, default=12.0)984 parser.add_argument("--cut", type=int, default=24)985 parser.add_argument("--left", type=float, default=-0.95)986 parser.add_argument("--right", type=float, default=2.0)987 parser.add_argument("--height", type=float, default=43.1)988 parser.add_argument("--seed", type=float, default=0.05)989 parser.add_argument("--rho", type=float, default=0.45)990 parser.add_argument("--eps", type=float, default=0.3)991 parser.add_argument("--hug", type=float, default=0.05)992 parser.add_argument("--ring", type=float, default=0.05)993 parser.add_argument("--span", type=int, default=18)994 parser.add_argument("--top", type=int, default=40)995 parser.add_argument("--deep", type=int, default=8)996 parser.add_argument("--guard", type=float, default=0.02)997 parser.add_argument("verb", choices=["census", "control", "teeth", "bridge", "dial"])998 args = parser.parse_args()999 {"census": census, "control": control, "teeth": teeth,1000 "bridge": bridge, "dial": dial}[args.verb](args)100110021003main()