fibre.py
Checks: READING: Scalings torn at one point rule.
Cited on: READING.
To run it: save fibre.py, with bonacci.py beside it; standard library only.
"""
fibre.py -- which torus points have one coding on the d-bonacci
completion, and the scalings n -> mn that one coding tears.
QUESTION. K_d, the completion of the trailing d-bonacci numeration
(places q_j = 2^j for j < d, q_k = q_(k-1) + ... + q_(k-d) after,
greedy digits = 0/1 strings with no d consecutive ones), maps onto the
torus T^(d-1) by phi(y) = sum y_k theta_k, theta_k = q_k alpha - nu_k,
alpha = (beta^-1, ..., beta^-(d-1)), so that phi(n) = n alpha
(bonacci.py). A coding of a torus point z is a y with phi(y) = z. If
n -> mn extends continuously to F, then phi F = m phi, and when z' has
one coding y while m z' has two, c and c', every neighbourhood of y
holds integers n with greedy(mn) near c and others with greedy(mn)
near c': n -> mn is torn. That needs one coding at 0 and integers
dense near z: the successor n -> n + 1 extends continuously, is onto
and injective off its fibre over 0, so z and z + alpha have equally
many codings for z != -alpha, and one coding at 0 gives one at every
n alpha. Which points are singly coded?
DESIGN.
- THE POINT FIBRE. Lift z = (g - r alpha)/m mod Z^(d-1). With M the
companion matrix, v_0 = (1, 1, 0, ..., 0) and r_beta the beta-
eigenvector, phi(y) = z iff the lattice point s_0 = T mod r_beta,
T = sum m y_j M^j v_0, has (s_0)_1 = -r and (s_0)_(i+1) = -g_i mod
m, and the sequence s_(k+1) = M^(-1)(s_k - m y_k v_0) stays
bounded. The bound: in eigen-coordinates each contracting
coordinate of s_k is m b_l times a legal tail sum_n y_(k+n) l^n,
which lies in a polygon Re(u w) <= h_l(u) for DIRS directions u
(support: an exact max over the first T terms by the run of ones,
plus a geometric tail below 1e-9); the beta-coordinate is bounded
by max(S0, m |b_beta|/(beta - 1)), S0 the bound (s_0)_1 = -r puts
on it at the start, which the recursion keeps. Two distinct starts with
infinite paths give distinct codings (their difference
M^(-k)(s_0 - s_0') stays bounded only at 0), and one start has
one path per label, so the codings are the infinite paths of this
graph, and z has one coding iff its live part (nodes with an
infinite continuation, found by peeling sinks) is one path.
- z' = (e_1 - alpha)/m is the torus limit of n_K alpha, n_K =
(q_K - 1)/m, along K = 0 mod the period of (q_K, ..., q_(K-d+1))
mod m, and m z' = -alpha, which has the d codings of Lemma V
(bonacci.py). So r = 1, g = e_1.
- Floating point: roots by Durand-Kerner, eigen-coordinates by
Gaussian elimination, every polygon side widened by 1e-7 and every
bound by a factor 1 + 1e-9 plus 1e-9. Widening only adds nodes, so
an error below the margins loses no coding; the margins are chosen,
not derived from an error bound.
PREDICTIONS AND CONTROLS (frozen before the run; N1 and N2 were frozen
at d = 2..10 with the kill 'any N1 or N2 failure at a d where P1
holds', and run at d <= 8 after a first --full spent minutes in N1
at d = 9).
P1 -alpha (m = 1, r = 1, g = 0) has exactly d live starts and no
branch node; 0 (r = 0) has one coding, the lemma's hypothesis.
N1 n alpha has one coding, greedy(n) on 60 digits, n = 1..30.
N2 -k alpha has exactly d codings, k = 1..5, read as the number of
live walks of length 200 from the live starts (the comparison
with length 400 added after the run, a strengthening). N1 and N2
at d <= 8.
N3 -alpha written at m >= 2, z = (0 - m alpha)/m, has d codings,
d = 2..6, m = 2..5. Its graph is m times the m = 1 graph, so it
reads only the code paths and margins at m >= 2, not a new point.
PLANT the box halved loses a coding of -alpha at some d = 2..6.
F3 z' has one coding at every pair read, so n -> mn is torn there.
P2 at a pair bonacci.py certifies, the coding found equals the
limit of A_K along K = 0 mod pi (y_j = W[(-j - 1) mod pi]) on
min(3 pi, 400) digits.
P3 the coding's low 40 digits equal greedy(n_K)'s at K = 0 mod the
state period, K the largest multiple of it up to 600, or the
period itself when longer, where the period is found below 5000.
The default run reads d = 2..6, --full d = 2..10, both m = 2..12 at
a node limit of 800000 per point.
FINDINGS (entered after the runs; prints copied from them).
- P1 holds at d = 2..10: -alpha d live starts, 0 branch nodes, 0
single (31658 + 28753 nodes at d = 10).
- N1, N2 hold at d = 2..8 (0 wrong; walks (d, d) at 200 and 400 at
every k, d = 7, 8 read by section_n alone after the --full), N3 at
d = 2..6, m = 2..5, as the successor's count says. The PLANT
fails as it must: the halved box loses a coding of -alpha at every
d = 2..6.
- F3: z' singly coded, so n -> mn TORN, at 76 of the 99 pairs
d = 2..10, m = 2..12; at no pair read two codings. 23 meet the
node limit: (6, 11), (7, 5), (7, 7), (7, 10), (7, 11), (8, 5),
(8, 7), (8, 10), (8, 11), (9, 5), (9, 6), (9, 7), (9, 10),
(9, 11), (9, 12), (10, 3), (10, 5), (10, 6), (10, 7), (10, 9),
(10, 10), (10, 11), (10, 12). The certificate tears 47, 46 of them
here; by the certificate alone (9, 6); torn by either 77 of 99.
P2 and P3 hold wherever run. The default run: 54 of 55 torn, (6,
11) unread; the certificate's 32 all among them.
- Tiers: the count and the lemma are proved (the argument above);
the tears are a finite check per pair given 0's single coding, itself P1's
finite check per d, in double precision with the chosen margins.
P2 ran at the 69 certified pairs read, P3 at the 75 torn pairs
whose state period was found ((5, 11) has none below 5000).
RUN RECORD: default 122/122, 16.3 s, peak commit 328.8 MB; --full
173/173, 471.0 s, peak commit 495.3 MB (before the 200/400 walk
check; d = 7, 8 rerun by section_n, 4/4), under a 512 MB memory guard.
"""
import cmath
import math
import sys
import time
from collections import deque
import bonacci
CHECKS = []
DIRS = 16
SHRINK = 1.0
def check(name, ok, detail=""):
CHECKS.append(bool(ok))
print(f" [{'PASS' if ok else 'FAIL'}] {name}" +
(f" -- {detail}" if detail else ""), flush=True)
# ------------------------------------------------------------ numeration
class Num:
def __init__(self, d):
self.d = self.dim = d
self.M = [[1] * d] + [[1 if j == i else 0 for j in range(d)]
for i in range(d - 1)]
self.v0 = (1, 1) + (0,) * (d - 2)
self.roots = roots(d)
self.beta = max(self.roots, key=abs).real
def places(self, n):
return bonacci.places(self.d, n)
def greedy(self, n, L):
s = bonacci.greedy(n, self.places(L))
return [1 if j in s else 0 for j in range(L)]
def roots(d):
"""Roots of x^d - x^(d-1) - ... - 1 by Durand-Kerner."""
p = [1] + [-1] * d # leading coefficient first
f = lambda z: sum(c * z ** (d - i) for i, c in enumerate(p))
z = [(0.4 + 0.9j) ** i for i in range(d)]
for _ in range(2000):
nz = []
for i, zi in enumerate(z):
den = 1
for j, zj in enumerate(z):
if i != j:
den *= zi - zj
nz.append(zi - f(zi) / den)
z = nz
return z
def solve(A, b):
n = len(A)
A = [list(map(complex, row)) + [complex(b[i])] for i, row in enumerate(A)]
for c in range(n):
p = max(range(c, n), key=lambda r: abs(A[r][c]))
A[c], A[p] = A[p], A[c]
for r in range(n):
if r != c:
f = A[r][c] / A[c][c]
A[r] = [x - f * y for x, y in zip(A[r], A[c])]
return [A[i][n] / A[i][i] for i in range(n)]
def eigen(N):
"""Roots l, the coordinates b_l of v_0 = sum b_l r_l, and the
inverse of R = (r_l) as rows: s = sum ell_l(s) r_l."""
d = N.dim
lam = N.roots
R = [[l ** (d - 1 - i) for l in lam] for i in range(d)]
b = solve(R, list(N.v0))
cols = [solve(R, [1 if j == i else 0 for j in range(d)])
for i in range(d)]
rows = [[cols[j][i] for j in range(d)] for i in range(d)]
return lam, b, rows
def Minv_col(N, u):
d = N.dim
v = list(u[1:]) + [0]
v[d - 1] = u[0] - sum(u[1:])
return tuple(v)
def support(l, c, d, u, eps=1e-9):
"""sup over legal 0/1 strings y of Re(u c sum_n y_n l^n): an exact
max over the first T terms (by the run of ones) plus a tail bound
sum_(n>=T) |c||l|^n below eps."""
a = abs(l)
T = max(1, math.ceil(math.log(eps * (1 - a) / abs(c)) / math.log(a)))
best = [0.0] + [-math.inf] * (d - 1) # by the current run
w = u * c
for n in range(T):
t = w.real
new = [max(best)] + [-math.inf] * (d - 1)
for r0 in range(d - 1):
if best[r0] > -math.inf:
new[r0 + 1] = best[r0] + t
best = new
w *= l
return max(best) + abs(c) * a ** T / (1 - a)
def fp_center(G, ctr, C):
"""Integer u with (u - ctr)^T G (u - ctr) <= C."""
n = len(G)
q = [row[:] for row in G]
for i in range(n):
for j in range(i + 1, n):
q[j][i] = q[i][j]
q[i][j] = q[i][j] / q[i][i]
for k in range(i + 1, n):
for l in range(k, n):
q[k][l] -= q[k][i] * q[i][l]
out, u = [], [0] * n
def rec(i, rem):
cen = ctr[i] - sum(q[i][j] * (u[j] - ctr[j])
for j in range(i + 1, n))
r = math.sqrt(max(rem, 0) / q[i][i])
for v in range(math.ceil(cen - r - 1e-9),
math.floor(cen + r + 1e-9) + 1):
u[i] = v
left = rem - q[i][i] * (v - cen) ** 2
if left < -1e-9:
continue
if i == 0:
out.append(tuple(u))
else:
rec(i - 1, left)
u[i] = 0
rec(n - 1, C)
return out
# ---------------------------------------------------- the point fibre
def point_fibre(N, m, r, g, limit=800_000, walk=0):
"""The codings y of z = (g - r alpha)/m mod Z^(d-1): the infinite
paths of s_(k+1) = M^(-1)(s_k - m y_k v_0) from the starts."""
d = N.dim
lam, b, rows = eigen(N)
ib = max(range(d), key=lambda i: abs(lam[i]))
beta = lam[ib].real
# each contracting coordinate m b_l sum y_(k+n) l^n of a legal tail
# lies in the polygon Re(u z) <= h_l(u), u over DIRS directions
dirs = [cmath.exp(2j * math.pi * t / DIRS) for t in range(DIRS)]
poly = {i: [(u, SHRINK * support(lam[i], m * 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 poly.items():
bnd[i] = max(h for _, h in hs) / math.cos(math.pi / DIRS)
S0 = (abs(r) + sum(bnd[i] * abs(lam[i]) ** (d - 1) for i in poly)) \
/ beta ** (d - 1)
bnd[ib] = SHRINK * max(S0, m * abs(b[ib]) / (beta - 1))
bnd = [x * (1 + 1e-9) + 1e-9 for x in bnd]
def inside(s):
for i in range(d):
z = sum(rows[i][j] * s[j] for j in range(d))
if abs(z) > bnd[i]:
return False
if i != ib and any((u * z).real > h for u, h in poly[i]):
return False
return True
# starts s_0 = c + m u, u_1 = 0, c_1 = -r, c_(i+1) = -g_i mod m
c = [-r] + [(-gi) % m for gi in g]
A, a = [], []
for i in range(d):
for part in (lambda z: z.real, lambda z: z.imag):
A.append([part(rows[i][j] * m) / bnd[i] for j in range(1, d)])
a.append(part(sum(rows[i][j] * c[j] for j in range(d))) / bnd[i])
n = d - 1
G = [[sum(A[k][i] * A[k][j] for k in range(len(A))) for j in range(n)]
for i in range(n)]
Ata = [sum(A[k][i] * a[k] for k in range(len(A))) for i in range(n)]
ctr = [-x.real for x in solve(G, Ata)]
const = sum(x * x for x in a) - sum(
ctr[i] * G[i][j] * ctr[j] for i in range(n) for j in range(n))
us = fp_center(G, ctr, d * (1 + 1e-9) - const)
starts = []
for u in us:
s0 = tuple([c[0]] + [c[j] + m * u[j - 1] for j in range(1, d)])
if inside(s0):
starts.append((s0, 0))
v0 = N.v0
idx = {nd: i for i, nd in enumerate(starts)}
nodes = list(starts)
succ = []
i = 0
while i < len(nodes):
s, run = nodes[i]
out = []
for y in (0, 1):
r1 = run + 1 if y else 0
if r1 >= d:
out.append(-1)
continue
s1 = Minv_col(N, tuple(si - m * y * vi for si, vi in zip(s, v0)))
if not inside(s1):
out.append(-1)
continue
nd = (s1, r1)
if nd not in idx:
idx[nd] = len(nodes)
nodes.append(nd)
if len(nodes) > limit:
return None
out.append(idx[nd])
succ.append(out)
i += 1
del idx, nodes
# live: nodes with an infinite continuation, by peeling sinks
nn = len(succ)
cnt = [sum(1 for t in o if t >= 0) for o in succ]
pred = [[] for _ in range(nn)]
for u, o in enumerate(succ):
for t in o:
if t >= 0:
pred[t].append(u)
live = [True] * nn
dq = deque(u for u in range(nn) if cnt[u] == 0)
while dq:
u = dq.popleft()
if not live[u]:
continue
live[u] = False
for p in pred[u]:
cnt[p] -= 1
if cnt[p] == 0 and live[p]:
dq.append(p)
del pred
ls = [i for i in range(len(starts)) if live[i]]
branch = sum(1 for u in range(nn) if live[u] and
sum(1 for t in succ[u] if t >= 0 and live[t]) > 1)
coding = None
if ls:
u, coding = ls[0], []
for _ in range(400):
t = [(y, t) for y, t in enumerate(succ[u]) if t >= 0 and live[t]]
coding.append(t[0][0])
u = t[0][1]
walks = None
if walk:
# live walks of length walk: every one extends, so this counts
# the codings' (start, prefix) pairs, the codings once it settles
w, half = {u: 1 for u in ls}, None
for _ in range(walk):
nw = {}
for u, k in w.items():
for t in succ[u]:
if t >= 0 and live[t]:
nw[t] = nw.get(t, 0) + k
w = nw
if _ == walk // 2 - 1:
half = sum(w.values())
walks = (half, sum(w.values()))
return dict(single=len(ls) == 1 and branch == 0, starts=len(starts),
live_starts=len(ls), branch=branch, nodes=nn,
live=sum(live), coding=coding, walks=walks)
# ------------------------------------------------------------ sections
def section_p1(ds):
print("P1 -alpha has d codings, 0 has one")
for d in ds:
N = Num(d)
t0 = time.time()
a = point_fibre(N, 1, 1, [0] * (d - 1))
z = point_fibre(N, 1, 0, [0] * (d - 1))
check(f"P1: d = {d}: -alpha {a['live_starts']} live starts, "
f"{a['branch']} branch nodes; 0 single",
a["live_starts"] == d and a["branch"] == 0 and z["single"],
f"{a['nodes']} + {z['nodes']} nodes, {time.time() - t0:.1f} s")
def section_n(ds):
print("N the successor's count along the orbit of 0 and of -alpha")
for d in ds:
N = Num(d)
bad1 = sum(1 for n in range(1, 31)
if not (f := point_fibre(N, 1, -n, [0] * (d - 1)))["single"]
or f["coding"][:60] != N.greedy(n, 60))
ws = [point_fibre(N, 1, k, [0] * (d - 1), walk=400)["walks"]
for k in range(1, 6)]
check(f"N1: d = {d}: n alpha has the one coding greedy(n), "
f"n = 1..30", bad1 == 0, f"{bad1} wrong")
check(f"N2: d = {d}: -k alpha has {d} codings, k = 1..5",
all(w == (d, d) for w in ws), f"walks at 200, 400: {ws}")
for d in [x for x in ds if x <= 6]:
N = Num(d)
ws = [point_fibre(N, m, m, [0] * (d - 1), walk=400)["walks"]
for m in range(2, 6)]
check(f"N3: d = {d}: -alpha at m = 2..5 has {d} codings",
all(w == (d, d) for w in ws), f"walks at 200, 400: {ws}")
def section_plant():
global SHRINK
print("PLANT the box halved")
SHRINK = 0.5
lost = [d for d in range(2, 7)
if point_fibre(Num(d), 1, 1, [0] * (d - 1))["live_starts"] < d]
SHRINK = 1.0
check("PLANT: the halved box loses a coding of -alpha", bool(lost),
f"at d = {lost}")
def section_f3(ds, ms):
print("F3 z' = (e_1 - alpha)/m", flush=True)
res, cert = {}, set()
for d in ds:
N = Num(d)
for m in ms:
t0 = time.time()
c = bonacci.certificate(d, m)
if c is not None and c["tear"]:
cert.add((d, m))
f = point_fibre(N, m, 1, [1] + [0] * (d - 2))
if f is None:
res[d, m] = "unread"
print(f" d = {d}, x{m}: UNREAD (node limit)", flush=True)
continue
res[d, m] = "TORN" if f["single"] else "open"
if c is not None and f["coding"]:
W, pi = c["W"], c["pi"]
want = [W[(-j - 1) % pi] for j in range(min(3 * pi, 400))]
check(f"P2: d = {d}, x{m}: the coding is A_K's class-0 limit",
f["coding"][:len(want)] == want)
pi = bonacci.state_period(d, m, bonacci.places(d, 6000), 100)
if pi and f["coding"]:
K = pi * max(1, 600 // pi)
q = N.places(K + 1)
nK = (q[K] - 1) // m
check(f"P3: d = {d}, x{m}: greedy(n_K) at K = {K} has the "
f"coding's low 40 digits", (q[K] - 1) % m == 0 and
N.greedy(nK, K)[:40] == f["coding"][:40])
print(f" d = {d}, x{m}: {res[d, m]} (starts "
f"{f['starts']}, live starts {f['live_starts']}, branch "
f"{f['branch']}, nodes {f['nodes']}, live {f['live']}, "
f"{time.time() - t0:.1f} s)", flush=True)
torn = sorted(k for k, v in res.items() if v == "TORN")
print(f" torn: {len(torn)} of {len(res)}; not torn: "
f"{sorted(k for k, v in res.items() if v != 'TORN')}")
both = cert & set(torn)
print(f" the certificate tears {len(cert)}, {len(both)} of them "
f"here; by the certificate alone {sorted(cert - both)}; torn by "
f"either {len(cert | set(torn))} of {len(res)}")
return res
def main():
full = "--full" in sys.argv
ds = range(2, 11) if full else range(2, 7)
section_p1(ds)
section_n([d for d in ds if d <= 8])
section_plant()
section_f3(ds, range(2, 13))
print(f"{sum(CHECKS)}/{len(CHECKS)} checks passed")
if __name__ == "__main__":
main()