""" pairs.py -- the pair automaton G_d: every two codings of one point of P_m, over every m at once, as the paths of one finite graph, and its live part as a carry-window automaton K_d. QUESTION. Two codings y, y' of one point start at s_0, s'_0 with equal first coordinate and the rest equal mod m; E_k = (s_k - s'_k)/m runs E_(k+1) = M^(-1)(E_k - (y_k - y'_k) v_0) from E_0 in {0} x Z^(d-1), integral and bounded, m gone. z is multiply coded iff one coding of z is the y-label of a non-diagonal infinite path of G_d (nodes (E, a, a'), a, a' the runs) from a start node ((0, e), 0, 0). What is G_d's live part, in one d-shape, and does it read sigma single and -sigma multiple, as fibre.py does? DESIGN. fibre.py's eigen-coordinates: E's contracting coordinate l lies in the polygon of b_l times a difference of legal tail sums, support h(u) + h(-u); the beta-coordinate is bounded by max(S0, |b_beta|/(beta - 1)), S0 the bound a first coordinate 0 puts on it at the start. G_d's live part by BFS in the box, in two stages. The RELAXED step (delta free in {-1, 0, 1}, no runs) on lattice points: BFS from the starts, the box's points with E_1 = 0, each successor kept iff inside; R_d, the reached points with an infinite continuation inside the reached set, is the box's relaxed live part reached from a start through live points (by hand: a continuation from a reached point is reached through it, and a path to a live point is live through it). The box's sides bound differences of legal tail sums; F checks each side's support over free delta is no larger, the beta-coordinate's bounds using only |delta| <= 1 and the sides, so the box holds every bounded relaxed path from a start and R_d is the relaxed live part outright. Then legality, the runs, on R_d: a legal pair's E's are relaxed live, so G_d's live part is the legal graph's over R_d, read on arrays, a node (E, a, a') one cell of a boolean mask over R_d x d x d, reached from the starts and peeled to its live part. Three cuts keep dead points out and lose no live one: the starts enumerated coordinate by coordinate, a partial vector kept while every row of the box's polytope can still be met by the rest (the slack by linear programming); joint cuts, since on a bounded relaxed path ell_l(E_t) = b_l sum_k delta_(t+k) lam_l^k, so Re sum_l u_l ell_l(E) <= sum_k |Re sum_l u_l b_l lam_l^k| for every complex u (300 random u); and a two-step lookahead, a point with no relaxed path of two steps inside being dead. A product with a lasso word u V^inf (y forced) adds the phase; a path is non-diagonal once it leaves E = 0 or starts off it, read as: the live product holds a start node with E != 0, or a node with E = 0 and y != y' with a live successor. A product run after G_d is seeded from G_d's live starts and kept in its live set, on its mask with a phase coordinate. The carry coordinates, derived by hand from the step E' = (E_2 - delta, E_3, ..., E_d, E_1 - E_2 - ... - E_d), delta = y - y'. With w = (E_2, ..., E_d, E_1 - E_2 - ... - E_d), so that E = (sum w, w_1, ..., w_(d-1)), the step is a shift: w' = (w_2, ..., w_d, c), c = 2 w_1 - sum w - delta, and sum w' = w_1 - delta. So a pair of codings is one integer carry sequence x, w its length-d windows, with delta_m = x_m - (x_(m+1) + ... + x_(m+d)): y - y' is, up to a term at the start, the carry 1^d 0 <-> 0^d 1 applied x_m times at position m. A start (0, e) is the window (e, -sum e), of sum 0. Setting y = y' = 0 wherever delta = 0 is never worse, so a continuation is legal iff delta's runs of +1 and of -1 stay below d, the first counted on from a or a'. K_d, the carry graph: nodes (w, a, a'), w in {-1, 0, 1}^d with sum in {-1, 0, 1}, those edges, starts (w, 0, 0) with sum w = 0: G_d's step on the E in {-1, 0, 1}^d with |E_1 - E_2 - ... - E_d| <= 1, with no box. Its live part reachable from its starts is G_d's iff every live E of G_d lies in {-1, 0, 1}^d (P1; the second bound then holds, an entry of each successor) and the box cuts no infinite path of K_d (a bounded path's contracting coordinates are tail sums, its beta-coordinate within the bound the recursion keeps). Every zero-sum w is a live start of K_d, at every d: delta = w_1 keeps the sum 0 and appends c = w_1, so the window rotates and the pair is y = [w = 1]^inf, y' = [w = -1]^inf, period d, legal as a nonzero zero-sum w holds both signs. The two words code one point by periodic.py's criterion: the E-path returns after d steps, so E_0 is the periodic solution, E_0 = -iota(X_y - X_y'), first coordinate 0. W, the triangle wave, by hand. Let x have period P = d - 1. Then x_(t+1) + ... + x_(t+d) is a period's sum S plus x_(t+d) = x_(t+1), so delta_t = x_t - x_(t+1) - S and the window at t sums to S + x_t. With S = 0 and x_0 = 0 the pair theorem asks only for steps of x of size <= 1. The wave climbs 0 -> a, falls a -> -a, climbs -a -> 0, then rests at 0 for P - 4a steps, a = floor(P/4); x_(t + 2a) = -x_t on the moving part, so S = 0. y = [delta = 1] = (0^a 1^(2a) 0^(P - 3a))^inf and y' = [delta = -1] = (1^a 0^(2a) 1^a 0^(P - 4a))^inf have runs of at most 2a < d, cyclically too; y has period P, so codes a point of some P_m (periodic codings), and y' codes it with carries reaching a. W1 reads each wave against the periodic-codings criterion (D = (M^P - I)^(-1) sum_k delta_k M^k v_0 integral, D_1 = 0) and walks its E-path; W2 reads the d = 10 pair's prefix difference D_n, D_n alpha's distance to Z^9 in high precision, against D_n + 1. PREDICTIONS, CONTROLS AND KILL (fixed before the runs). C1 CONTROL, read first, d = 4..8: the S-product has no non-diagonal live path from a start node (sigma single, as fibre.py reads), and the A-product has one, its start E_0 = +-e_3 (by hand: s_A - s_B = -4 iota(beta^2 - beta - 1) = -4 e_3). A mismatch voids the run. C2 CONTROL, d = 3, 4: every census class fibre.py labels multiple at m = 2..4 has its first coding paired, every single one not. PLANT: the box halved loses the A/B pair at some d. P1 the live E-set's entries lie in {-1, 0, 1} at every d read. P2 (TRANSPLANT from tilings, no confidence) 2(2^(d-1) - 1) nonzero E's among the live start nodes. KILL of the d-uniform route: P1 fails at d = 8 with an entry of size >= 3, or the live E count at d = 8 exceeds 10^4. Fixed after d = 3..7 were read as a fit, before d = 8 was run: K1 at d = 8, G_8's live set read in w is K_8's live part reachable from its starts, symmetric difference 0. K2 at d = 8, G_8's live starts are the zero-sum windows. Fixed after the carry bound failed at d = 9 and the wave was proved, before the BFS was built: B, the largest carry and window sum over G_d's live part is max(1, floor((d - 1)/4)) at d = 3..9; its control, the route reproducing the box enumeration's sets at d = 3..9 (7213 relaxed points, 67875 live nodes at d = 9). Z (TRANSPLANT from the carry graph), R_d is the carry graph's windows and +-e_1 at d = 3..8, e_1 the finite carry 0^d 1 <-> 1^d 0, which needs d ones in a row. Fixed after the run, from a review: F, each side's free support within its legal one at d = 3..9. B, C1 and F carried to d = 10 after the mask route had read d = 10 outside this script, its control the set route's live part at d = 5..10, cell for cell. CONTROL, read first: K_d with c's delta sign flipped misses, a nonzero symmetric difference at d = 3, 4, 5. KILL of the rule: a nonzero symmetric difference at d = 8. FINDINGS (entered after the runs; prints copied from them). - C1 holds at d = 4..10: the S-product has no non-diagonal live path on G_d's live part, the A-product's starts at E_0 = (0, 0, -1, 0, ..., 0) = -e_3. C2 holds at d = 3, 4, m = 2..4: paired iff labelled multiple, no exception. The PLANT loses the A pair at every d = 4..8. - KILL misses; P1 holds at d = 3..8: every live E has entries in {-1, 0, 1}; it fails at d = 9 (B below). The route at d = 3..9: starts 7, 19, 51, 165, 726, 4017, 27280 past the cuts, reached 15, 37, 95, 319, 1471, 8635, 60385, R_d 15, 37, 95, 255, 699, 1947, 7213; live nodes 33, 130, 463, 1566, 5117, 16402, 67875 on 13, 35, 93, 253, 697, 1945, 7211 E's; live nonzero starts 6, 18, 50, 140, 392, 1106 at d = 3..8. P2 holds at d = 3 alone (6) and fails from d = 4 (18 against 14). - W1 holds: every wave at d = 5..30 passes the criterion, its E-path closing, its largest carry a (7 at d = 30). W2: the d = 10 pair (001111000, 110000110), |D_n alpha| 0.00266, 8.16e-12, 4.97e-23 at n = 200, 800, 1600, control 0.5. - B holds at d = 3..10: largest carry and window sum 1 at d = 3..8, 2 at d = 9, 10, over G_d's live part, and largest carry the same over R_d. At d = 10: starts 198558, reached 450083 in 21 layers, R_d 20649, live nodes 215218 on 20647 E's. Its control, frozen as the box enumeration's sets at d = 3..9, was read before this route replaced the enumeration as less: the relaxed sets equal at d = 3..7, the counts at d = 9 (7213, 67875); K1 and Z read the sets at d = 3..8 against the carry graph. The starts and reached counts move with the joint cuts drawn; R_d and G_d cannot, the cuts losing no live point. F holds at d = 3..10: the largest excess of a side's free support over its legal one is -1.01e-07, the sides' own widening. Z holds at d = 3..8: R_d = the carry graph's windows and +-e_1, 15 against 13 up to 1947 against 1945. - The rule's KILL misses; K1 and K2 hold at d = 3..8, d = 8 the frozen test: G_d's live part is K_d's reachable live part, 33, 130, 463, 1566, 5117, 16402 nodes in both; the flipped control misses by 30, 126, 458. So the live starts are the zero-sum windows of {-1, 0, 1}^d, T_d - 1 nonzero ones, T_d the central trinomial coefficient: the 6, 18, 50, 140, 392, 1106 above. RUN RECORD: the bare run, d = 3..10 on masks, 74.9 s, peak commit 276.8 MB; d = 10's route 36.0 s, its products 5.8 s. Before the masks, d = 3..9 took 42.8 s, 414.4 MB (the route 24.7 s, d = 9 8.4 s; the products 1.5 s, seeded from G_d's live starts and kept in its live set; the Python BFS it replaced took 184.9 s at d = 8, the box enumeration 293 s at d = 9). The mask route at d = 11 took 393.2 s and 421.3 MB alone, its products 465.6 MB, so the bare run stops at 10. --plant d = 4..8 4.4 s, 24.0 MB. """ import cmath from fractions import Fraction import itertools from array import array import math import os import sys import time os.environ.setdefault("OPENBLAS_NUM_THREADS", "1") import numpy as np # noqa: E402 from fibre import DIRS, Minv_col, Num, census, eigen, support # noqa: E402 CHECKS = [] def check(name, ok, detail=""): CHECKS.append(bool(ok)) print(f" [{'PASS' if ok else 'FAIL'}] {name}" + (f" -- {detail}" if detail else ""), flush=True) class Box: """The box of E: polygons for the contracting coordinates, a bound on the beta-coordinate.""" def __init__(self, N, shrink=1.0): d = self.d = N.dim lam, b, rows = eigen(N) self.rows = rows ib = self.ib = max(range(d), key=lambda i: abs(lam[i])) beta = lam[ib].real dirs = [cmath.exp(2j * math.pi * t / DIRS) for t in range(DIRS)] self.poly = {i: [(u, shrink * (support(lam[i], b[i], d, u) + support(lam[i], b[i], d, -u)) + 1e-7) for u in dirs] for i in range(d) if i != ib} bnd = [0.0] * d for i, hs in self.poly.items(): bnd[i] = max(h for _, h in hs) / math.cos(math.pi / DIRS) S0 = sum(bnd[i] * abs(lam[i]) ** (d - 1) for i in self.poly) \ / beta ** (d - 1) bnd[ib] = shrink * max(S0, abs(b[ib]) / (beta - 1)) self.bnd = [x * (1 + 1e-9) + 1e-9 for x in bnd] self.lam, self.b = lam, b R = [[l ** (d - 1 - i) for l in lam] for i in range(d)] self.top = [int(sum(abs(R[j][l]) * self.bnd[l] for l in range(d))) + 1 for j in range(d)] def inside(self, s): for i in range(self.d): z = sum(self.rows[i][j] * s[j] for j in range(self.d)) if abs(z) > self.bnd[i]: return False if i != self.ib and any((u * z).real > h for u, h in self.poly[i]): return False return True def starts(self): """The lattice points (0, e) inside: e over the coordinate box that the eigen-bounds give, |s_j| <= sum_l |r_l,j| bnd_l.""" out = [] for e in itertools.product(*[range(-t, t + 1) for t in self.top[1:]]): s = (0,) + e if self.inside(s): out.append(s) return out class Nodes: """The nodes (E, a, a', phase) packed as integers, decoded on read: G_d at d = 8 holds about 10^6 of them.""" def __init__(self, d, rad, P): self.d, self.rad, self.P, self.codes = d, rad, P, [] def pack(self, nd): E, a, a2, ph = nd c = 0 for e in E: c = c * self.rad + e + self.rad // 2 return ((c * self.d + a) * self.d + a2) * self.P + ph def __getitem__(self, i): c = self.codes[i] c, ph = divmod(c, self.P) c, a2 = divmod(c, self.d) c, a = divmod(c, self.d) E = [] for _ in range(self.d): c, e = divmod(c, self.rad) E.append(e - self.rad // 2) return tuple(E[::-1]), a, a2, ph def __len__(self): return len(self.codes) def build(N, box, word=None): """G_d, or its product with the lasso word (u, V): nodes (E, a, a', phase); an edge labelled (y, y') to node t stored as 4t + 2y + y'. Returns nodes, succ, the number of starts.""" d, v0 = N.dim, N.v0 seq = None if word is not None: u, V = word seq = u + V def nxt(ph): return ph + 1 if ph + 1 < len(seq) else len(u) nodes = Nodes(d, 2 * max(box.top) + 1, len(seq) if seq else 1) idx = {} for E in box.starts(): c = nodes.pack((E, 0, 0, 0)) idx[c] = len(nodes.codes) nodes.codes.append(c) nst = len(nodes.codes) succ = [] i = 0 while i < len(nodes.codes): E, a, a2, ph = nodes[i] out = [] ys = (0, 1) if seq is None else (seq[ph],) for y in ys: for y2 in (0, 1): r1, r2 = (a + 1 if y else 0), (a2 + 1 if y2 else 0) if r1 >= d or r2 >= d: continue E1 = Minv_col(N, tuple(e - (y - y2) * w for e, w in zip(E, v0))) if not box.inside(E1): continue c = nodes.pack((E1, r1, r2, 0 if seq is None else nxt(ph))) t = idx.get(c) if t is None: t = idx[c] = len(nodes.codes) nodes.codes.append(c) out.append(4 * t + 2 * y + y2) succ.append(out) i += 1 return nodes, succ, nst def live_part(succ): """fibre.peel on flat arrays: the nodes with an infinite continuation (the predecessor lists of G_d at d = 8 pass 512 MB as lists).""" nn = len(succ) cnt = array("i", (len(o) for o in succ)) off = array("i", [0]) * (nn + 1) for o in succ: for x in o: off[(x >> 2) + 1] += 1 for u in range(nn): off[u + 1] += off[u] pred, fill = array("i", [0]) * off[nn], array("i", off[:nn]) for u, o in enumerate(succ): for x in o: t = x >> 2 pred[fill[t]] = u fill[t] += 1 del fill live = bytearray([1]) * nn stack = [u for u in range(nn) if cnt[u] == 0] while stack: u = stack.pop() if not live[u]: continue live[u] = 0 for k in range(off[u], off[u + 1]): q = pred[k] cnt[q] -= 1 if cnt[q] == 0 and live[q]: stack.append(q) return live def paired(nodes, succ, nst, live): """A non-diagonal live path from a start: a live start off E = 0, or a live node on E = 0 reachable along the diagonal with a live edge y != y'.""" zero = tuple([0] * len(nodes[0][0])) off = [i for i in range(nst) if live[i] and nodes[i][0] != zero] if off: return ("start", nodes[off[0]][0]) seen = [i for i in range(nst) if live[i] and nodes[i][0] == zero] stack, mark = list(seen), set(seen) while stack: u = stack.pop() for x in succ[u]: t, y, y2 = x >> 2, (x >> 1) & 1, x & 1 if not live[t]: continue if y != y2: return ("branch", nodes[u]) if t not in mark: mark.add(t) stack.append(t) return None def words(d): S = [1] + [0] * (d + 1) + [1] * (d - 3) + [0, 1, 0] A = ([0, 0] + [1] * (d - 1) + [0] * (d - 2) + [1, 0, 0, 0] + [1] * (d - 2) + [0, 1, 1, 0] + [1] * (d - 4) + [0, 0, 1]) return S, A def lasso(y, top=150): """y (400 digits) as (u, V) with y = u V^inf on its digits.""" for p in range(1, top): s = max((i + 1 for i in range(len(y) - p) if y[i] != y[i + p]), default=0) if s < len(y) - 2 * p - 50: return y[:s], y[s:s + p] return None def w_of(E): """The carry window of E: (E_2, ..., E_d, E_1 - E_2 - ... - E_d).""" return tuple(E[1:]) + (E[0] - sum(E[1:]),) def windows(E): """w_of on the rows of an array.""" return np.hstack([E[:, 1:], E[:, :1] - E[:, 1:].sum(axis=1, keepdims=True)]) def carry(d, flip=False): """K_d's live part reachable from its starts: nodes (w, a, a'), w in {-1, 0, 1}^d with sum in {-1, 0, 1}, the step w -> (w_2, ..., w_d, c), c = 2 w_1 - sum w - delta, sum w' = w_1 - delta; starts (w, 0, 0) with sum w = 0. With flip, c = 2 w_1 - sum w + delta (the control).""" idx, nodes, succ = {}, [], [] for w in itertools.product((-1, 0, 1), repeat=d): if sum(w) == 0: idx[(w, 0, 0)] = len(nodes) nodes.append((w, 0, 0)) i = 0 while i < len(nodes): w, a, a2 = nodes[i] s, out = sum(w), [] for y in (0, 1): for y2 in (0, 1): r1, r2 = (a + 1 if y else 0), (a2 + 1 if y2 else 0) dl = y - y2 c = 2 * w[0] - s + (dl if flip else -dl) if r1 >= d or r2 >= d or abs(c) > 1 or abs(w[0] - dl) > 1: continue nd = (w[1:] + (c,), r1, r2) t = idx.get(nd) if t is None: t = idx[nd] = len(nodes) nodes.append(nd) out.append(4 * t + 2 * y + y2) succ.append(out) i += 1 live = live_part(succ) return {nodes[i] for i in range(len(nodes)) if live[i]} class VBox: """Box's test on arrays of E, its linear constraints A E <= b, and joint cuts every bounded relaxed path obeys.""" def __init__(self, B, joint=0, seed=1): d = self.d = B.d self.rows = np.array(B.rows, dtype=complex) self.bnd = np.array(B.bnd) self.cs = [i for i in range(d) if i != B.ib] self.U = {i: np.array([u for u, _ in B.poly[i]]) for i in self.cs} self.H = {i: np.array([h for _, h in B.poly[i]]) for i in self.cs} self.top = np.array(B.top) A, b = [], [] for i in self.cs: for u, h in B.poly[i]: A.append((u * self.rows[i]).real) b.append(h) A.append(self.rows[B.ib].real) b.append(B.bnd[B.ib]) A.append(-self.rows[B.ib].real) b.append(B.bnd[B.ib]) self.J, self.h = self.joint(B, joint, seed) A, b = np.array(A + list(self.J)), np.array(b + list(self.h)) key = np.round(np.hstack([A, b[:, None]]), 9) _, keep = np.unique(key, axis=0, return_index=True) self.A, self.b = A[np.sort(keep)], b[np.sort(keep)] def slack(self, first): """slack[j]_i = -min A_i[j:] . E[j:] over the polytope A E <= b, E_1 = first: what the coordinates from j on can take off row i, by linear programming; and self.range, each coordinate's integer range there, kept as [-range, range]: a live point lies in -P as well, the box being symmetric and each joint cut holding for -u too.""" from scipy.optimize import linprog d, A, b = self.d, self.A, self.b eq = (np.eye(1, d), [first]) bd = [(None, None)] * d def lo(c): r = linprog(c, A_ub=A, b_ub=b, A_eq=eq[0], b_eq=eq[1], bounds=bd, method="highs") return r.fun self.range = np.array([int(np.floor(-lo(-np.eye(1, d, j)[0]) + 1e-9)) for j in range(d)]) out = [] for j in range(d + 1): c = np.zeros((len(A), d)) c[:, j:] = A[:, j:] out.append(np.array([-lo(row) if j < d else 0.0 for row in c])) return out def joint(self, B, n, seed): """n joint directions over the contracting coordinates: on a bounded relaxed path ell_l(E_t) = b_l sum_k delta_(t+k) lam_l^k, so Re sum_l u_l ell_l(E) <= sum_k |Re sum_l u_l b_l lam_l^k| for any complex u, as rows J E <= h.""" if not n: return np.zeros((0, self.d)), np.zeros(0) rng = np.random.default_rng(seed) lam = np.array([B.lam[i] for i in self.cs]) bb = np.array([B.b[i] for i in self.cs]) K = int(np.ceil(np.log(1e-13) / np.log(np.max(np.abs(lam))))) pw = bb[None, :] * lam[None, :] ** np.arange(K)[:, None] U = rng.normal(size=(n, len(lam))) + 1j * rng.normal(size=(n, len(lam))) J = (U @ self.rows[self.cs]).real h = np.abs((pw @ U.T).real).sum(axis=0) + 1e-7 return J, h def inside(self, E, chunk=20000): ok = np.empty(len(E), dtype=bool) for s in range(0, len(E), chunk): F = E[s:s + chunk].astype(float) Z = F @ self.rows.T o = np.all(np.abs(Z) <= self.bnd, axis=1) o &= np.all(F @ self.J.T <= self.h, axis=1) for i in self.cs: o &= np.all(np.outer(Z[:, i], self.U[i]).real <= self.H[i], axis=1) ok[s:s + chunk] = o return ok def starts(self, keep, chunk=4000): """The lattice points inside with E_1 = 0 that keep passes, as an int8 array: depth first by chunks, each partial vector kept while every row of A E <= b can still be met (slack).""" d, A, b = self.d, self.A, self.b slack = self.slack(0) top = self.range found = [] def grow(P, j): if j == d: P = P[self.inside(P)] found.append(P[keep(P)]) return vals = (np.array([0], dtype=np.int8) if j == 0 else np.arange(-top[j], top[j] + 1, dtype=np.int8)) step_ = max(1, chunk * 66 // len(A)) for s in range(0, len(P), step_): Q = P[s:s + step_] Q = np.hstack([np.repeat(Q, len(vals), axis=0), np.tile(vals, len(Q))[:, None]]) S = Q.astype(float) @ A[:, :j + 1].T grow(Q[np.all(S - slack[j + 1] <= b + 1e-9, axis=1)], j + 1) grow(np.zeros((1, 0), dtype=np.int8), 0) return np.vstack(found) def viable(V, Q, k, top): """Which rows of Q have a relaxed path of k steps inside the box: a point without one is dead, so never live.""" ok = np.zeros(len(Q), dtype=bool) if k == 0: return ~ok for dl in (-1, 0, 1): idx = np.flatnonzero(~ok) if not len(idx): break R = step(Q[idx], dl) inb = np.all(np.abs(R) <= top, axis=1) inb[inb] = V.inside(R[inb]) idx, R = idx[inb], R[inb] ok[idx[viable(V, R, k - 1, top)]] = True return ok def reach(V, chunk=500000, look=0): """The box's points reached from its starts (the slice E_1 = 0) by the relaxed step, each successor kept iff inside and with look steps more inside: BFS by layers, the seen set a sorted array of keys.""" top = V.top.astype(np.int16) S = V.starts(lambda P: viable(V, P.astype(np.int16), look, top)) k = keys(S, top) k, i = np.unique(k, return_index=True) seen, pts, front, layers = k, [S[i]], S[i], 0 while len(front): layers += 1 new = [] for s in range(0, len(front), chunk): for dl in (-1, 0, 1): Q = step(front[s:s + chunk], dl) Q = Q[np.all(np.abs(Q) <= top, axis=1)] Q = Q[V.inside(Q)] new.append(Q[viable(V, Q, look, top)].astype(np.int8)) Q = np.vstack(new) kq, i = np.unique(keys(Q, top), return_index=True) pos = np.minimum(np.searchsorted(seen, kq), len(seen) - 1) fresh = seen[pos] != kq front = Q[i[fresh]] pts.append(front) seen = np.union1d(seen, kq[fresh]) return np.vstack(pts), len(S), layers def keys(E, top): """E as int64 keys, coordinate j offset by top_j + 1.""" k = np.zeros(len(E), dtype=np.int64) for j in range(E.shape[1]): k = k * (2 * int(top[j]) + 3) + (E[:, j].astype(np.int64) + int(top[j]) + 1) return k def step(E, dl): """E' = (E_2 - delta, E_3, ..., E_d, E_1 - E_2 - ... - E_d).""" E = E.astype(np.int16) out = np.empty_like(E) out[:, 0] = E[:, 1] - dl out[:, 1:-1] = E[:, 2:] out[:, -1] = E[:, 0] - E[:, 1:].sum(axis=1) return out def relaxed_live(V, P, chunk=500000): """R_d: the relaxed step on the box's points P; the points with an infinite continuation, reachable from the starts (E_1 = 0). Chunked, small dtypes.""" top = V.top.astype(np.int16) k = keys(P, top) order = np.argsort(k) P, k = P[order], k[order] del order n = len(P) succ = np.full((n, 3), -1, dtype=np.int32) for s in range(0, n, chunk): for c, dl in enumerate((-1, 0, 1)): Q = step(P[s:s + chunk], dl) inb = np.all(np.abs(Q) <= top, axis=1) Q[~inb] = 0 kq = keys(Q, top) pos = np.minimum(np.searchsorted(k, kq), n - 1) hit = inb & (k[pos] == kq) succ[s:s + chunk, c][hit] = pos[hit] del Q, kq, pos, hit, inb del k has = succ >= 0 safe = np.where(has, succ, 0) live = np.ones(n, dtype=bool) rounds = 0 while True: new = live & (live[safe] & has).any(axis=1) rounds += 1 if (new == live).all(): break live = new del safe seen = live & (P[:, 0] == 0) front = np.flatnonzero(seen) while len(front): nx = succ[front].ravel() nx = nx[nx >= 0] nx = np.unique(nx[live[nx] & ~seen[nx]]) seen[nx] = True front = nx return P[seen], rounds def legal_dense(d, R): """The legal graph over R on arrays: node (i, a, a') as (i d + a) d + a', i a row of R; reached from the starts ((0, e), 0, 0), then peeled. Returns the live mask, shape (len(R), d, d).""" top = np.abs(R).max(axis=0).astype(np.int16) + 2 k = keys(R, top) order = np.argsort(k) R, k = R[order], k[order] n = len(R) nxt = {} for dl in (-1, 0, 1): Q = step(R, dl) inb = np.all(np.abs(Q) <= top, axis=1) kq = keys(np.where(inb[:, None], Q, 0), top) pos = np.minimum(np.searchsorted(k, kq), n - 1) nxt[dl] = np.where(inb & (k[pos] == kq), pos, -1).astype(np.int32) a = np.arange(d) moves = [] for y in (0, 1): for y2 in (0, 1): r1 = a + 1 if y else np.zeros(d, dtype=int) r2 = a + 1 if y2 else np.zeros(d, dtype=int) moves.append((nxt[y - y2], r1, r2)) def succ(M): """Each move's image of the node mask M.""" out = np.zeros_like(M) for t_, r1, r2 in moves: ok1, ok2 = r1 < d, r2 < d sub = M[:, ok1][:, :, ok2] i, x, z = np.nonzero(sub) j = t_[i] h = j >= 0 out[j[h], r1[ok1][x[h]], r2[ok2][z[h]]] = True return out reach = np.zeros((n, d, d), dtype=bool) reach[R[:, 0] == 0, 0, 0] = True front = reach.copy() while front.any(): new = succ(front) & ~reach reach |= new front = new live = reach.copy() while True: has = np.zeros_like(live) for t_, r1, r2 in moves: a1, a2 = np.flatnonzero(r1 < d), np.flatnonzero(r2 < d) tgt = live[np.where(t_ >= 0, t_, 0)][:, r1[a1]][:, :, r2[a2]] has[np.ix_(np.arange(n), a1, a2)] |= tgt & (t_ >= 0)[:, None, None] new = live & has if (new == live).all(): break live = new return R, live def paired_dense(d, R, L, V): """build and paired on arrays: the product of G_d's live part (R, L) with y = V^inf, phase ph, y = V[ph], y' free; a non-diagonal live path from a start, as paired reads it: ("start", E), ("branch", (E, a, a', ph)) or None.""" n, P = len(R), len(V) gid = np.full(n * d * d, -1, dtype=np.int32) flat = np.flatnonzero(L.ravel()) gid[flat] = np.arange(len(flat), dtype=np.int32) i, a, a2 = np.unravel_index(flat, (n, d, d)) top = np.abs(R).max(axis=0).astype(np.int16) + 2 k = keys(R, top) gs = {} for y in (0, 1): for y2 in (0, 1): Q = step(R[i], y - y2) kq = keys(np.clip(Q, -top, top), top) pos = np.minimum(np.searchsorted(k, kq), n - 1) ok = (k[pos] == kq) & np.all(np.abs(Q) < top, axis=1) r1 = np.where(y == 1, a + 1, 0) r2 = np.where(y2 == 1, a2 + 1, 0) ok &= (r1 < d) & (r2 < d) tf = (pos.astype(np.int64) * d + np.minimum(r1, d - 1)) * d + \ np.minimum(r2, d - 1) gs[y, y2] = np.where(ok, gid[tf], -1) NG = len(flat) start = (R[i, 0] == 0) & (a == 0) & (a2 == 0) def images(M): out = np.zeros_like(M) for ph in range(P): src = np.flatnonzero(M[:, ph]) for y2 in (0, 1): tg = gs[V[ph], y2][src] out[tg[tg >= 0], (ph + 1) % P] = True return out reach = np.zeros((NG, P), dtype=bool) reach[start, 0] = True front = reach.copy() while front.any(): front = images(front) & ~reach reach |= front live = reach.copy() while True: has = np.zeros_like(live) for ph in range(P): for y2 in (0, 1): tg = gs[V[ph], y2] has[:, ph] |= (tg >= 0) & live[np.maximum(tg, 0), (ph + 1) % P] new = live & has if (new == live).all(): break live = new zero = np.all(R[i] == 0, axis=1) off = np.flatnonzero(start & ~zero & live[:, 0]) if len(off): return ("start", tuple(int(x) for x in R[i[off[0]]])) stack = [(g, 0) for g in np.flatnonzero(start & zero & live[:, 0])] mark = set(stack) while stack: g, ph = stack.pop() for y2 in (0, 1): tg = gs[V[ph], y2][g] if tg < 0 or not live[tg, (ph + 1) % P]: continue if y2 != V[ph]: return ("branch", (tuple(int(x) for x in R[i[g]]), int(a[g]), int(a2[g]), ph)) if (tg, (ph + 1) % P) not in mark: mark.add((tg, (ph + 1) % P)) stack.append((tg, (ph + 1) % P)) return None def companion(d): return [[1] * d] + [[1 if j == i else 0 for j in range(d)] for i in range(d - 1)] def columns(d, P): """den and the integer columns den * (M^P - I)^(-1) M^k v_0, k < P.""" M = companion(d) A = [[Fraction(int(i == j)) for j in range(d)] for i in range(d)] for _ in range(P): A = [[sum(A[i][k] * M[k][j] for k in range(d)) for j in range(d)] for i in range(d)] cols, c = [], [1, 1] + [0] * (d - 2) for _ in range(P): cols.append(c) c = [sum(M[i][j] * c[j] for j in range(d)) for i in range(d)] A = [[A[i][j] - (i == j) for j in range(d)] + [Fraction(cols[k][i]) for k in range(P)] for i in range(d)] det = Fraction(1) for c in range(d): r = next(r for r in range(c, d) if A[r][c] != 0) if r != c: A[c], A[r], det = A[r], A[c], -det det *= A[c][c] A[c] = [x / A[c][c] for x in A[c]] for r in range(d): if r != c and A[r][c] != 0: A[r] = [x - A[r][c] * y for x, y in zip(A[r], A[c])] den = abs(int(det)) C = [[int(A[i][d + k] * den) for i in range(d)] for k in range(P)] return den, C def legal(word, d): """No d consecutive ones in the word read cyclically.""" if all(word): return False s = "".join(map(str, word)) * 2 return max(len(r) for r in s.split("0")) < d def hp(d, y, y2, ns=(200, 800, 1600)): """D_n alpha's distance to Z^(d-1), D_n = sum_(k Z^9", float(rows[-1]) < 1e-15 and float(ctl) > 0.1, f"at n = 200, 800, 1600: {rows}; control D_n + 1: {ctl}") def section_k(alive, starts): print("K G_d's live part against K_d's, read in w", flush=True) for d in sorted(alive): if d <= 5: bad = carry(d, flip=True) ^ {(w_of(E), a, a2) for E, a, a2 in alive[d]} check(f"control d = {d}: the flipped step misses", bad, f"symmetric difference {len(bad)}") for d in sorted(alive): G = {(w_of(E), a, a2) for E, a, a2 in alive[d]} K = carry(d) zs = {(0,) + w[:-1] for w in itertools.product((-1, 0, 1), repeat=d) if sum(w) == 0} check(f"K1 d = {d}: G_d live = K_d reachable live", G == K, f"G {len(G)}, K {len(K)}, symmetric difference {len(G ^ K)}") check(f"K2 d = {d}: live starts = the zero-sum windows", starts[d] == zs, f"{len(zs) - 1} nonzero, symmetric " f"difference {len(starts[d] ^ zs)}") def free_gap(B): """The largest excess, over the box's polygon sides, of the support of b_l sum_n delta_n lam_l^n with delta free in {-1, 0, 1}^N, sum_n |Re(u b_l lam_l^n)|, over the side's bound h (legal differences): at most 0, the sides hold every bounded relaxed path, whose contracting coordinates are such sums.""" gap = -math.inf for i, sides in B.poly.items(): lam, c = B.lam[i], B.b[i] T = math.ceil(math.log(1e-13 * (1 - abs(lam)) / abs(c)) / math.log(abs(lam))) pw = [c * lam ** n for n in range(T)] for u, h in sides: gap = max(gap, sum(abs((u * z).real) for z in pw) - h) return gap def section_g(top): """G_d's live part by BFS in its box: the relaxed set R_d (delta free in {-1, 0, 1}, the box the only filter), then legality on it.""" print("G_d the live part, by BFS in the box", flush=True) alive, starts, dense = {}, {}, {} for d in range(3, top + 1): N, t0 = Num(d), time.time() B = Box(N) gap = free_gap(B) check(f"F d = {d}: the box's sides hold every bounded relaxed " f"path, each side's free support within its legal one", gap <= 1e-9, f"largest excess {gap:.2e}") V = VBox(B, joint=300) P, ns, layers = reach(V, look=2) R, _ = relaxed_live(V, P) Rs, L = dense[d] = legal_dense(d, R) j, a, a2 = np.nonzero(L) G = {(tuple(int(x) for x in Rs[u]), int(p), int(q)) for u, p, q in zip(j, a, a2)} if d <= 8: alive[d] = G starts[d] = {E for E, p, q in G if E[0] == 0 and p == q == 0} i = np.unique(j) big = int(np.abs(windows(Rs[i])).max()) bs = int(np.abs(Rs[i, 0]).max()) Rw = {w_of(tuple(int(x) for x in E)) for E in R} rbig = int(np.abs(windows(R)).max()) print(f" d = {d}: starts {ns}, reached {len(P)} in {layers} layers;" f" R_d {len(R)}, max |carry| {rbig}; live nodes {len(G)} on " f"{len(i)} windows, max |carry| {big}, max |window sum| {bs}" f" [{time.time() - t0:.1f} s]", flush=True) c = max(1, (d - 1) // 4) check(f"B d = {d}: largest carry and window sum max(1, " f"floor((d - 1)/4)) = {c}", big == bs == rbig == c) if d <= 8: K = {w for w, a, a2 in carry(d)} e1 = tuple([1] + [0] * (d - 1)) check(f"Z d = {d}: R_d = the carry graph's windows and +-e_1", Rw ^ K == {e1, tuple(-x for x in e1)}, f"{len(Rw)} against {len(K)}") return alive, starts, dense def main(): shrink = 0.5 if "--plant" in sys.argv else 1.0 top = (int(sys.argv[sys.argv.index("--d") + 1]) if "--d" in sys.argv else 10 if shrink == 1.0 else 8) if shrink == 1.0: section_wave() alive, starts, dense = section_g(top) section_k(alive, starts) if "--rule" in sys.argv: print(f"{sum(CHECKS)}/{len(CHECKS)} checks passed") sys.exit(0 if all(CHECKS) else 1) print("C1 sigma's S and -sigma's A against G_d", flush=True) for d in range(4, top + 1): N, t0 = Num(d), time.time() S, A = words(d) if shrink == 1.0: rs = [paired_dense(d, *dense[d], S), "on G_d's live part"] ra = [paired_dense(d, *dense[d], A), ""] else: box, rs, ra = Box(N, shrink), [], [] for w, acc in ((S, rs), (A, ra)): nodes, succ, nst = build(N, box, ([], w)) live = live_part(succ) acc.append(paired(nodes, succ, nst, live)) acc.append((len(nodes), sum(live))) e3 = tuple(-1 if j == 2 else 0 for j in range(d)) ok = (rs[0] is None and ra[0] is not None and ra[0][0] == "start" and ra[0][1] in (e3, tuple(-x for x in e3))) check(f"d = {d}: S unpaired, A paired at E_0 = +-e_3", ok, f"S {rs[0]} {rs[1]}, A {ra[0]} {ra[1]} " f"[{time.time() - t0:.1f} s]") if shrink != 1.0: print(f"{sum(CHECKS)}/{len(CHECKS)} checks passed") sys.exit(0 if all(CHECKS) else 1) print("C2 census classes at d = 3, 4, m = 2..4", flush=True) for d in (3, 4): N = Num(d) box = Box(N) bad = [] for m in (2, 3, 4): out = census(N, m) for lab in ("single", "multiple"): for key, f in out[lab]: w = lasso(f["coding"]) nodes, succ, nst = build(N, box, w) pr = paired(nodes, succ, nst, live_part(succ)) if (pr is not None) != (lab == "multiple"): bad.append((m, key, lab, pr)) check(f"d = {d}: paired iff labelled multiple, m = 2..4", not bad, f"{bad[:4]}") print(f"{sum(CHECKS)}/{len(CHECKS)} checks passed") sys.exit(0 if all(CHECKS) else 1) if __name__ == "__main__": main()