pairs.py
Checks: READING: Pairs of codings theorem; READING: The carry bound rule.
Cited on: READING.
To run it: save pairs.py, with bonacci.py, fibre.py beside it; standard library only.
"""
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. The
start nodes: the lattice points (0, e) inside, enumerated over the
bounding box of the polygon bounds. G_d by BFS, its live part by
peeling sinks (fibre.peel). 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.
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.
SK: the products of K_d with a lasso word u V^inf (y forced, y'
free), in exact integers, a depth-first search for the longest
non-diagonal path, one arm at a time with its own memo: the START arm
from the nonzero starts, the BRANCH arm off the diagonal at E = 0 at
every time up to len(u) + 2 len(V); a node met again on the stack is
an infinite path. Words: S^inf and A^inf; R = (1^(d-1) 0)^inf, a coding
of -alpha, which has d codings, its partners found on the start
arm; Q, read off K_d's live set (branch_word: an edge
y != y' from the zero start, then live edges to a repeat), whose
partner branches at time 0.
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.
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.
SK's control: the A-product does not die, -sigma being multiply
coded. Added after the S-product's read, before the run that
ships: R lives on the start arm, Q on the branch arm.
FINDINGS (entered after the runs; prints copied from them).
- C1 holds at d = 4..8: the S-product has no non-diagonal live path
(live 10, 12, 14, 16, 18 nodes, the cycle alone), 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}. G_d: nodes 33, 232, 1225, 13428, 92777, 1144068; live
33, 130, 463, 1566, 5117, 16402; live E's 13, 35, 93, 253, 697,
1945; live nonzero starts 6, 18, 50, 140, 392, 1106. P2 holds at
d = 3 alone (6) and fails from d = 4 (18 against 14).
- 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.
- SK: inside K_d the S-product has no infinite non-diagonal path at
d = 4..11 on either arm, its longest from a nonzero start 7, 13,
17, 26, 33, 40, 56, 63 steps and off the diagonal 15, 23, 29, 36,
44, 55, 65, 82; the A-product has one on the start arm at each d
(C1's start at d = 4..8). The controls hold at d = 4..11: R lives
on the start arm, Q on the branch arm. Inside K_d is all of G_d's live part where K1 holds
(d <= 8); past it the read covers the second codings whose carries
and window sums lie in {-1, 0, 1}.
RUN RECORD: the bare run, d = 3..8, 204.1 s, peak commit 326.5 MB
(G_8 184.9 s; the products 0.3 s seeded from G_d's live starts and kept
in its live set, the box check skipped as that set lies in the box,
against 316.1 s seeded from the box). SK at d = 12 passes 512 MB in its
memo, so it stops at 11 (9.1 s with its controls). --plant d = 4..8
4.5 s, 15.3 MB.
"""
import cmath
import itertools
from array import array
import math
import sys
import time
from fibre import DIRS, Minv_col, Num, census, eigen, support
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, allow=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'.
With allow, only nodes whose (E, a, a') is in it: a live product path
projects to a live G_d path. 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 = {}
seeds = (box.starts() if allow is None else
sorted(E for E, a, a2 in allow if (a, a2, E[0]) == (0, 0, 0)))
for E in seeds:
if allow is None or (E, 0, 0) in allow:
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 allow is not None:
if (E1, r1, r2) not in allow:
continue
elif 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 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]}
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 s_life(d, word, arm):
"""The product of K_d with the lasso word (u, V) (y = u V^inf forced,
y' free): the longest non-diagonal path of one arm, "start" (from a
nonzero start) or "branch" (off the diagonal at E = 0, at every time
up to len(u) + 2 len(V), past which time repeats phase and run), None
if one is infinite (a node met again on the stack). Each arm has its
own memo, so a cycle one arm met cannot hide from the other."""
u, V = word
seq, P, memo, inf = u + V, len(u) + len(V), {}, []
def nxt(ph):
return ph + 1 if ph + 1 < P else len(u)
def life(w, a2, ph):
key = (sum((x + 1) * 3 ** i for i, x in enumerate(w)) * d
+ a2) * P + ph
if key in memo:
if memo[key] is None:
inf.append(key)
return 0
return memo[key]
memo[key] = None
s, out = sum(w), 0
for y2 in (0, 1):
dl, r2 = seq[ph] - y2, (a2 + 1 if y2 else 0)
c = 2 * w[0] - s - dl
if r2 < d and abs(c) <= 1 and abs(w[0] - dl) <= 1:
out = max(out, 1 + life(w[1:] + (c,), r2, nxt(ph)))
memo[key] = out
return out
sys.setrecursionlimit(10000)
worst = 0
if arm == "start":
worst = max(life(w, 0, 0) for w in itertools.product((-1, 0, 1),
repeat=d)
if sum(w) == 0 and any(w))
else:
ph, run = 0, 0
for _ in range(len(u) + 2 * len(V)):
y2 = 1 - seq[ph]
if not y2 or run + 1 < d:
worst = max(worst, 1 + life((0,) * (d - 1) + (y2 - seq[ph],),
run + 1 if y2 else 0, nxt(ph)))
run = run + 1 if seq[ph] else 0
ph = nxt(ph)
return None if inf else worst
def branch_word(d):
"""A word whose second coding branches off the diagonal, read off K_d's
live set: from the zero start an edge with y != y' into a live node,
then the first live edge at each node until a node repeats, the labels
y read as (u, V)."""
K = carry(d)
nd, seen, path = ((0,) * d, 0, 0), {}, []
while nd not in seen:
seen[nd] = len(path)
w, a, a2 = nd
nexts = [(y, (w[1:] + (2 * w[0] - sum(w) - (y - y2),),
a + 1 if y else 0, a2 + 1 if y2 else 0))
for y in (0, 1) for y2 in (0, 1) if path or y != y2]
y, nd = next((y, t) for y, t in nexts if t in K)
path.append(y)
return path[:seen[nd]], path[seen[nd]:]
def section_sk(top):
print("SK sigma's S-product inside K_d, exact integers, by arm",
flush=True)
for d in range(4, top + 1):
t0 = time.time()
S, A = words(d)
R = [1] * (d - 1) + [0]
got = {(nm, arm): s_life(d, V, arm) for nm, V in
(("S", ([], S)), ("A", ([], A)), ("R", ([], R)),
("Q", branch_word(d))) for arm in ("start", "branch")}
txt = ", ".join(f"{nm} {arm} {'inf' if v is None else v}"
for (nm, arm), v in got.items())
check(f"d = {d}: S dies on both arms; A and R live on the start "
f"arm, Q on the branch arm", got[("S", "start")] is not None and
got[("S", "branch")] is not None and
got[("A", "start")] is None and
got[("R", "start")] is None and got[("Q", "branch")] is None,
f"{txt} [{time.time() - t0:.1f} s]")
def section_g(top):
print("G_d the live part", flush=True)
p1, alive, starts = [], {}, {}
for d in range(3, top + 1):
N, t0 = Num(d), time.time()
box = Box(N)
nodes, succ, nst = build(N, box)
live = live_part(succ)
alive[d] = {nodes[i][:3] for i in range(len(nodes)) if live[i]}
Es = sorted({nd[0] for nd in alive[d]})
starts[d] = {nodes[i][0] for i in range(nst) if live[i]}
st = sorted(starts[d] - {tuple([0] * d)})
big = max(max(abs(x) for x in E) for E in Es)
if big > 1:
p1.append((d, big))
print(f" d = {d}: nodes {len(nodes)}, live {sum(live)}, live E's "
f"{len(Es)}, live nonzero starts {len(st)} (P2 "
f"{2 * (2 ** (d - 1) - 1)}), max |entry| {big} "
f"[{time.time() - t0:.1f} s]", flush=True)
if d <= 5:
print(f" starts {st}", flush=True)
print(f" E's {Es}", flush=True)
del nodes, succ, live
print(f"P1: {'holds' if not p1 else 'fails ' + str(p1)}")
section_k(alive, starts)
return alive
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 8)
alive = section_g(top) if shrink == 1.0 else {}
if shrink == 1.0:
section_sk(11)
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()
box = Box(N, shrink)
S, A = words(d)
rs, ra = [], []
for w, acc in ((S, rs), (A, ra)):
nodes, succ, nst = build(N, box, ([], w), alive.get(d))
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()