pairs.py
Checks: READING: Pairs of codings theorem; READING: The triangle wave theorem; READING: The carry bound rule.
Cited on: READING.
To run it: save pairs.py, with bonacci.py, fibre.py beside it; needs mpmath, numpy, scipy.
"""
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<n) (y_k - y'_k) q_k,
the words read periodically; and D_n + 1's, the control."""
import mpmath
mpmath.mp.dps = 700
b = mpmath.findroot(lambda z: z ** d - sum(z ** k for k in range(d)),
mpmath.mpf(2))
q = [1]
while len(q) < max(ns):
q.append(sum(q[max(0, len(q) - d):]) + (1 if len(q) < d else 0))
P = len(y)
def dist(D):
return max(abs(D * b ** -j - mpmath.nint(D * b ** -j))
for j in range(1, d))
rows = []
for n in ns:
D = sum((int(y[k % P]) - int(y2[k % P])) * q[k] for k in range(n))
rows.append(mpmath.nstr(dist(D), 3))
return rows, mpmath.nstr(dist(D + 1), 3)
def wave(d):
"""The triangle wave of period P = d - 1, amplitude floor(P / 4), as
delta_t = x_t - x_(t+1)."""
P, a = d - 1, (d - 1) // 4
x = (list(range(a + 1)) + list(range(a - 1, -a - 1, -1)) +
list(range(-a + 1, 1)) + [0] * (P - 4 * a))
return [x[t] - x[t + 1] for t in range(P)], a
def section_wave(dmax=30):
"""W1: each wave's pair against the periodic-codings criterion, its
E-path closing, its largest carry a; W2: the d = 10 pair in high
precision, against the control D_n + 1."""
print("W the triangle wave of period d - 1", flush=True)
bad = []
for d in range(5, dmax + 1):
dl, a = wave(d)
P = len(dl)
den, C = columns(d, P)
num = [sum(dl[k] * C[k][j] for k in range(P)) for j in range(d)]
y = [int(v == 1) for v in dl]
y2 = [int(v == -1) for v in dl]
ok = (all(v % den == 0 for v in num) and num[0] == 0 and
legal(y, d) and legal(y2, d))
E, mc = [-v // den for v in num], 0
for k in range(P):
w = E[1:] + [E[0] - sum(E[1:])]
mc = max(mc, max(abs(v) for v in w), abs(E[0]))
u = [E[0] - dl[k], E[1] - dl[k]] + E[2:]
E = u[1:] + [u[0] - sum(u[1:])]
closes = E == [-v // den for v in num]
if not (ok and closes and mc == max(1, a)):
bad.append((d, ok, closes, mc, a))
if d in (5, 9, 13, 30):
print(f" d = {d}: y = ({''.join(map(str, y))})^inf, y' = "
f"({''.join(map(str, y2))})^inf, criterion {ok}, path "
f"closes {closes}, max |carry| {mc}, a = {a}", flush=True)
check(f"W1 d = 5..{dmax}: every wave a pair, its largest carry a",
not bad, f"failures {bad}")
dl, a = wave(10)
y = "".join(str(int(v == 1)) for v in dl)
y2 = "".join(str(int(v == -1)) for v in dl)
rows, ctl = hp(10, y, y2)
check(f"W2 d = 10, ({y}, {y2}): D_n alpha -> 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()