periodic.py
Checks: READING: Periodic codings theorem.
Cited on: READING.
To run it: save periodic.py, with bonacci.py, fibre.py beside it; needs sympy.
"""
periodic.py -- the codings that are periodic words, read exactly in
Q(beta): the class -sigma = (2, 1, 2, 0, ..., 0) of P_4 has the codings
A^inf != B^inf at every d >= 4, and sigma = (2, 3, 2, 0, ..., 0) the
coding S^inf.
QUESTION. fibre.py codes a point z = (g - r alpha)/m of P_m by the
paths s_(k+1) = M^(-1)(s_k - m y_k v_0) from a lattice start s_0 with
(s_0)_1 = -r, (s_0)_(i+1) = -g_i mod m: y codes z iff its states stay
bounded (each contracting eigen-coordinate is expanded by M^(-1) unless
it equals m b_l times y's tail sum, and the beta-coordinate is
contracted). Its float graphs read -sigma multiply coded and sigma
singly coded at d = 4..10, each by one d-shape. Do the shapes hold at
every d, read exactly?
THE THEOREM. Let V be a word of length P with no d consecutive ones
read cyclically, psi(V) = sum v_k beta^k, X_V = beta psi(V)/(beta^P -
1). If m X_V lies in Z[beta], V^inf codes the point of P_m with (r, g) =
iota(m X_V) EXACTLY, iota(xi) = xi(M) e_1, M with first row ones and
subdiagonal ones (Z[beta] -> Z^d, beta acting as M). Proof: M^k v_0 =
iota(beta^(k+1)), so s_P = s_0 reads (M^P - I) s_0 = -iota(m beta
psi(V)), whose solution s_0 = -iota(m X_V) is a lattice point; from it
the states repeat with period P, so they stay bounded, and V^inf codes
the point its start names. Two such words code one point iff r is equal
and g agrees mod m, that is iff X_A - X_B lies in Z[beta] with
lambda(X_A - X_B) = 0, lambda the first iota coordinate; and words
differing as sequences are two codings.
The words, low digit first, one period each:
S = 1 0^(d+1) 1^(d-3) 0 1 0, period 2d + 2;
A = U W, B = W U, U = 00 1^(d-1) 0^(d-2) 1 00, W = 0 1^(d-2) 0 11
0 1^(d-4) 00 1, period 4d + 4 (B is A shifted by 2d + 2).
EVERY d (FORMAL). In Q(beta), beta^d = 1/(2 - beta) at every d,
x^(d+1) - 2x^d + 1 being (x - 1) times the d-bonacci polynomial. The
words are runs of lengths c0 + c1 d (RUNS), so t psi(V) - X (t^P - 1),
psi summed run by run, becomes at t^d = u = 1/(2 - t) a rational
function of t alone, its poles at t = 1, 2 only; where it is 0 the
period relation holds at t = beta for every d >= 4, and with it X_V =
X. So X_A = (2 beta^2 - beta - 1)/4, X_B = (-2 beta^2 + 3 beta + 3)/4,
X_S = (2 beta^2 - 3 beta + 1)/4 at every d >= 4, and iota(beta) = e_1 +
e_2, iota(beta^2) = 2e_1 + e_2 + e_3 give iota(4 X_A) = (2, 1, 2, 0,
..., 0), iota(4 X_B) = (2, 1, -2, 0, ..., 0), iota(4 X_S) = (2, -1, 2,
0, ..., 0): A and B start at one r = 2 with g equal mod 4, and X_A -
X_B = beta^2 - beta - 1 = iota^(-1)(e_3).
LEGALITY at every d >= 4: A's runs of ones are 1^(d-1), 1, 1^(d-2),
11, 1^(d-4), 1 and S's 1, 1^(d-3), 1, each run of zeros between them
nonempty, and A's wrap (W ends 001, U starts 00) and B's joins (U ends
00, W ends 1 and U starts 00) merge no runs: no d consecutive ones. A
!= B: U and W differ at index 1. So -sigma = (2, 1, 2, 0, ..., 0) is
multiply coded at every d >= 4, by A^inf and B^inf, and S^inf codes a
point of sigma. That S^inf is sigma's only coding is the float
graphs' read alone, d = 4..10.
DESIGN. Elements of Q(beta) as coefficient vectors on 1, beta, ...,
beta^(d-1) in Fractions, beta^d = beta^(d-1) + ... + 1; division by
beta^P - 1 by solving its multiplication matrix exactly. CONTROLS
read fibre.read_class's codings (one per live start, each start's
least live path on 400 digits) and labels at d = 4..10. FORMAL in
sympy.
PREDICTIONS AND KILL (frozen before the exact run).
KILL: any d = 4..40 where A or B is illegal, X_A - X_B is not in
Z[beta] or has lambda != 0, or a key iota(4 X_V) mod 4 differs
from (2, 1, 2, 0, ..., 0) for A, B and (2, 3, 2, 0, ..., 0) for S.
CONTROLS, read first: at d = 4..8 (4..10 since, the d = 9, 10 rows
added when the script was carried, predicted to hold) the engine's
codings equal the words built here, and iota(4 X_V) of each engine
coding is its class's key; the PLANT, A reversed against A, reads
not one point.
FORMAL (frozen after the exact run): the formal difference is 0 for
A, B and S with the X above. CONTROL d = 4..10: the runs are the
words, and the claimed X meets the period relation at t = beta.
PLANT: A's 1^(d-2) cut to 1^(d-3) leaves a nonzero difference.
FINDINGS.
- CONTROLS 15/15: at every d = 4..10 sigma single with the one
coding S, -sigma multiple with the codings A and B, the d = 9, 10
rows as predicted; no branch node in -sigma's graph, so A^inf and
B^inf are its only codings there; iota(4 X_V) of each engine
coding its key.
- KILL misses: A, B, S legal, A and B one point, keys (2, 1, 2, 0,
..., 0), (2, 1, 2, 0, ..., 0) and S's (2, 3, 2, 0, ..., 0) at every
d = 4..40, exact; the plant reads two points.
- FORMAL: the formal difference is 0 for A, B and S; X_A - X_B =
t^2 - t - 1; the plant leaves t(-2t^3 + t^2 - 3)/(4(t^2 - 4t + 4));
the controls' relative error at most 7e-16 (8e-16 in the lab run,
beta then found by bisection).
- So the theorem above holds as stated: -sigma has the codings
A^inf != B^inf at every d >= 4, so at least two; exactly two is
the float graphs' read, d = 4..10.
RUN RECORD: in the lab, the controls at d = 4..8 with d = 4..40 6.0 s,
peak commit 13.2 MB, and FORMAL 0.6 s, 52.5 MB, the codings compared
on their first period; carried, the
bare run 27/27, 30.8 s, peak commit 218.4 MB (d = 10's reads 21.6 s);
with the 400-digit comparison 27/27, 29.2 s, 214.6 MB; with -sigma's
branch nodes read, 27/27, 30.3 s, 214.5 MB.
"""
import sys
import time
from fractions import Fraction
from fibre import Num, label, read_class
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)
def words(d):
"""S (sigma), A and B (-sigma), low digit first, one period each."""
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])
h = 2 * d + 2
return S, A, A[h:] + A[:h]
def legal(V, d):
run = 0
for y in V + V:
run = run + 1 if y else 0
if run >= d:
return False
return True
def mul_beta(x):
"""x beta, beta^d = beta^(d-1) + ... + 1."""
top = x[-1]
return [top] + [x[i - 1] + top for i in range(1, len(x))]
def mul(x, y):
out, p = [Fraction(0)] * len(x), list(y)
for c in x:
if c:
out = [o + c * q for o, q in zip(out, p)]
p = mul_beta(p)
return out
def power(d, P):
x = [Fraction(1)] + [Fraction(0)] * (d - 1)
for _ in range(P):
x = mul_beta(x)
return x
def divide(x, y):
"""x / y in Q(beta): solve the matrix of multiplication by y."""
d = len(x)
cols, e = [], [Fraction(1)] + [Fraction(0)] * (d - 1)
for _ in range(d):
cols.append(mul(e, y))
e = mul_beta(e)
A = [[cols[j][i] for j in range(d)] + [x[i]] for i in range(d)]
for c in range(d):
k = next(i for i in range(c, d) if A[i][c] != 0)
A[c], A[k] = A[k], A[c]
pv = A[c][c]
A[c] = [v / pv for v in A[c]]
for i in range(d):
if i != c and A[i][c] != 0:
f = A[i][c]
A[i] = [a - f * b for a, b in zip(A[i], A[c])]
return [A[i][d] for i in range(d)]
def X(V, d):
"""beta psi(V)/(beta^P - 1)."""
psi = [Fraction(0)] * d
p = [Fraction(1)] + [Fraction(0)] * (d - 1)
for y in V:
p = mul_beta(p)
if y:
psi = [a + b for a, b in zip(psi, p)]
den = power(d, len(V))
den[0] -= 1
return divide(psi, den)
def iota(x):
"""x(M) e_1, M with first row ones and subdiagonal ones."""
d = len(x)
v, out = [Fraction(1)] + [Fraction(0)] * (d - 1), [Fraction(0)] * d
for c in x:
out = [o + c * a for o, a in zip(out, v)]
v = [sum(v)] + v[:-1]
return out
def key(V, d):
"""iota(4 X_V) mod 4, or None when 4 X_V is not in Z[beta]."""
x = [4 * c for c in X(V, d)]
if any(c.denominator != 1 for c in x):
return None
return tuple(int(c) % 4 for c in iota(x))
def one_point(V, W, d):
diff = [a - b for a, b in zip(X(V, d), X(W, d))]
return all(c.denominator == 1 for c in diff) and iota(diff)[0] == 0
# run-length words, each run (bit, (c0, c1)) of length c0 + c1 d
U = [(0, (2, 0)), (1, (-1, 1)), (0, (-2, 1)), (1, (1, 0)), (0, (2, 0))]
W = [(0, (1, 0)), (1, (-2, 1)), (0, (1, 0)), (1, (2, 0)), (0, (1, 0)),
(1, (-4, 1)), (0, (2, 0)), (1, (1, 0))]
RUNS = {"A": U + W, "B": W + U,
"S": [(1, (1, 0)), (0, (1, 1)), (1, (-3, 1)), (0, (1, 0)),
(1, (1, 0)), (0, (1, 0))]}
CLAIM = {"A": "(2*t**2 - t - 1)/4", "B": "(-2*t**2 + 3*t + 3)/4",
"S": "(2*t**2 - 3*t + 1)/4"}
def expand(runs, d):
return [b for b, (c0, c1) in runs for _ in range(c0 + c1 * d)]
def formal(runs, X, t, u):
"""t psi(V) - X (t^P - 1) with t^d = u: psi summed block by block."""
psi, a = 0, (0, 0)
for b, (c0, c1) in runs:
if b:
psi += t ** a[0] * u ** a[1] * (t ** c0 * u ** c1 - 1) / (t - 1)
a = (a[0] + c0, a[1] + c1)
return t * psi - X * (t ** a[0] * u ** a[1] - 1)
def inf(V, n=400):
"""V^inf on n digits, the length of an engine coding."""
return (V * (n // len(V) + 1))[:n]
def section_controls():
print("CONTROLS the engine's codings against the words", flush=True)
for d in range(4, 11):
N, t0 = Num(d), time.time()
S, A, B = words(d)
sp, ns = (2, 3, 2) + (0,) * (d - 3), (2, 1, 2) + (0,) * (d - 3)
fs, fn = read_class(N, 4, sp), read_class(N, 4, ns)
es = [y[:len(S)] for y in fs["codings"]]
en = [y[:len(A)] for y in fn["codings"]]
check(f"d = {d}: sigma {label(fs)}, coding S^inf; -sigma "
f"{label(fn)}, codings A^inf, B^inf",
label(fs) == "single" and fs["codings"] == [inf(S)] and
label(fn) == "multiple" and fn["branch"] == 0 and
sorted(fn["codings"]) == sorted([inf(A), inf(B)]),
f"[{time.time() - t0:.1f} s]")
check(f"d = {d}: iota(4 X_V) of the engine codings is the key",
[key(V, d) for V in es] == [sp] and
{key(V, d) for V in en} == {ns})
S, A, B = words(6)
check("PLANT d = 6: A reversed against A, not one point",
not one_point(A, A[::-1], 6))
def section_exact():
print("EXACT d = 4..40", flush=True)
kill = []
for d in range(4, 41):
t0 = time.time()
S, A, B = words(d)
sp, ns = (2, 3, 2) + (0,) * (d - 3), (2, 1, 2) + (0,) * (d - 3)
ok = (legal(A, d) and legal(B, d) and legal(S, d) and A != B and
one_point(A, B, d) and key(A, d) == ns and key(B, d) == ns
and key(S, d) == sp)
if not ok:
kill.append(d)
if d <= 6 or d % 6 == 4 or not ok:
print(f" d = {d}: {'holds' if ok else 'FAILS'}, legal "
f"{legal(A, d)}, {legal(B, d)}, {legal(S, d)}; one point "
f"{one_point(A, B, d)}; keys {key(A, d)}, {key(B, d)}, "
f"S {key(S, d)} [{time.time() - t0:.1f} s]", flush=True)
check("KILL misses: A, B one point of -sigma, S in sigma, d = 4..40",
not kill, f"fires at {kill}")
def section_formal():
import sympy
print("FORMAL every d at once, t^d = 1/(2 - t)", flush=True)
t = sympy.Symbol("t")
for d in range(4, 11):
b = Num(d).beta
same, rel = True, 0.0
for name, runs in RUNS.items():
same &= expand(runs, d) == dict(zip("SAB", words(d)))[name]
x = float(sympy.sympify(CLAIM[name]).subs(t, b))
P = len(expand(runs, d))
rel = max(rel, abs(formal(runs, x, b, b ** d)) / (x * b ** P))
check(f"CONTROL d = {d}: the runs are the words; the claimed X "
f"meets the period relation at t = beta, relative {rel:.0e}",
same and rel < 1e-12)
bad = [r if i != 1 else (1, (-3, 1)) for i, r in enumerate(W)]
pl = sympy.simplify(formal(U + bad, sympy.sympify(CLAIM["A"]), t,
1 / (2 - t)))
check("PLANT: A with 1^(d-2) cut to 1^(d-3) leaves a nonzero difference",
pl != 0, f"{pl}")
for name, runs in RUNS.items():
r = sympy.simplify(formal(runs, sympy.sympify(CLAIM[name]), t,
1 / (2 - t)))
print(f" {name}: X = {CLAIM[name]}, formal difference {r}",
flush=True)
check(f"{name}: identity in t", r == 0)
diff = sympy.expand(sympy.sympify(CLAIM["A"]) - sympy.sympify(CLAIM["B"]))
print(f" X_A - X_B = {diff}", flush=True)
def main():
section_controls()
section_exact()
section_formal()
print(f"{sum(CHECKS)}/{len(CHECKS)} checks passed")
sys.exit(0 if all(CHECKS) else 1)
if __name__ == "__main__":
main()