fibre.py
Checks: READING: The fibre lemma theorem; READING: Scalings torn at one point rule; READING: Multiply coded classes 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? 0 is, at every d, by the
cycles of the graph below; is z'? And which points of denominator m,
z' among them, tear x n?
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, whose codings are the d limits of Lemma V
(bonacci.py), and only those (THE CYCLES). So r = 1, g = e_1.
- THE CYCLES. At m = 1, with iota(xi) = xi(M) e_1 (Z[beta] -> Z^d,
M^k e_1, k < d, unitriangular), M^k v_0 = iota(beta^(k+1)), so a
closed walk through s reading V = v_0 ... v_(p-1) gives
(M^p - I) s = -sum v_k M^k v_0 and Z = psi(V)/(beta^p - 1) =
-iota^(-1)(s)/beta in Z[beta], psi(V) = sum v_k beta^k. Its digit
string has no d consecutive ones, so it is Z's greedy expansion
unless a tail is (1^(d-1)0)^inf (Parry 1960), and property (F)
(Frougny and Solomyak 1992, a_1 >= ... >= a_d >= 1) makes that
expansion finite: V is 0^p at s = 0 (M^p - I invertible, beta
Pisot) or a power of a rotation of 1^(d-1)0, off s = 0 (Z > 0). On
a pattern the node's run fixes the next label, and reading u^j, u
the rotation from that node, the
state is -iota(beta psi(u)/(beta^d - 1)) whatever j, so the walk is
back after d. The graph is deterministic, so each live SCC holding a
cycle is that one cycle; cycle_words reads it once around from its
least node (Tarjan), kind() names it zero, pattern, other or tangle.
A coding y != 0 of 0 repeats a node past its lowest 1, and the
periodic coding it gives is an integer n >= 1 (phi = n alpha) or
lim (q_K - r), r >= 1 (phi = -r alpha): neither is 0. So 0 has one
coding at every d, and by the successor every n alpha one and every
-k alpha exactly d: the counts of P1, N1 and N2 at every d (P1's d
distinct live starts and N1's digits stay checks). At m >= 2 the
state s_k/m runs the m = 1 recursion on (1/m) Z^d, its cycles giving
psi(V)/(beta^p - 1) in (1/m) Z[beta], out of (F)'s reach: C3's
word is such a cycle.
- THE CLASSES. z = (g - r alpha)/m depends on (r mod m, g mod m) up
to alpha: g = g' mod m gives the same torus point and z + alpha =
(g - (r - m) alpha)/m. So P_m, the points of denominator m, falls
into m^d classes, each one alpha-orbit, -alpha in the zero class,
and by the successor every point of a nonzero class has its class's
count, read at r = r_0, or m when r_0 = 0 (read_class). If c is
singly coded and n c = (n r, n g) mod m is multiply coded or zero,
x n is torn: at r >= 1, n z lies in n c, and in the zero class it is
-(n r/m) alpha, which has d codings. So (E1) a singly coded nonzero
class of P_m tears x km for every k >= 1, and (E4) for n prime to m
the class n^(-1) w of a multiply coded w depends on n only mod m:
one read tears a residue class of n.
- CENSUS labels each nonzero class single (one live start, no branch
node), multiple, or unread at the node limit; uncoded would be a
fault of the box, every torus point having a coding. O1 reads
Mult_d(m), the multiply coded classes, against M mod m, which
carries (r, g) to (r + sum g, r, g_1, ...): the real-linear B with
theta_(k+1) = B theta_k (B alpha = alpha - e_1, B e_j = e_(j+1) -
alpha) carries class kappa to M kappa and phi(0y) = B phi(y), so
y -> 0y carries the codings of a point of kappa to codings of
points of M kappa, onto one point only when the lifts' integer parts
v share their coordinate sum (B v = -(sum v) alpha mod Z^(d-1)):
invariance predicted, not proved. Counts are the live walks at 100
and 200 the census already reads. T_d(m), the residues n with n
kappa multiply coded or zero for some singly coded kappa: by the
fibre lemma's (3) every x n with n mod m in T_d(m) is torn.
- SHORT. s_(k+1) = M^(-1) s_k mod m, so a class's graph lies over its
orbit under M mod m (E2, a heuristic for its size). At a pair F3
leaves unread, the classes of the first three L whose kernel of
M^L - I mod p, p prime, leaves the earlier span (a multiple of L*
repeats L*'s
kernel), up to 60 per L, read to the first single (--short d:p).
- PREIMAGE. z' = n^(-1) w over the multiply coded short-orbit classes
w of P_m', m' = 2, 3, 5, 7, 11 with n mod m' neither 0 nor 1 (u = 1
gives z' = w), to the first single (--preimage d:n); n z' = w -
j alpha, j = (n r' - r_w)/m', is checked on the torus below 1e-9 at
each witness. At m' = 4, n = 3, z' = -w, every n = 3 mod 4 torn by
E4: at d >= 4 w = -sigma = (2, 1, 2, 0, ..., 0), multiply coded at
every d >= 4 (periodic.py), at d = 3 the first multiply coded w of
P_4 in itertools order whose -w is single (--neg d,...).
The bare and --full runs read the witnesses these searches found
(SHORT, PREIMAGE), and CLOSURE adds E1 and E4 to the tears.
- 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.
C2 0's live part is one node, with its loop 0, and -alpha's cycles
are zero or pattern, at P1's d.
C1 every live cycle of N1's and N2's graphs is zero or a pattern.
C3 POSITIVE CONTROL: at m = 2, F3's graph of z' holds a live cycle
that is neither, a rotation of 0 1^(d-1) 0, d = 2..6 (--full
2..10).
VOID: any C1 or C2 failure, printed as the cycle's kind and word, is
a fault of the code or of the proof, and voids both until found; C3
or the cycle PLANT missed voids the run.
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; the
cycle reader on hand-built graphs reads 1100 'other', 1101
'pattern', a loop at 0 'zero', and two tangles 'tangle', one at
state 0 whose first edges read 0 (with the tangle test removed in
a scratch copy it reads that one 'zero', and the plant FAILS).
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.
H1 every nonzero class of P_m is singly coded; H0 some has two.
Predicted H0 somewhere, no confidence: the Rauzy boundary is an
uncountable fractal (a TRANSPLANT). KILL of H1: any nonzero class
printed multiply coded. Arms: d = 2, m = 2..12; d = 3, m = 2..12;
d = 4, m = 2..7; d = 5, m = 2..4 (the default d = 2, 3, m <= 8).
CONTROLS, read first: -alpha at (m, 0) reads d codings and is
labelled multiple; the class (1, e_1), z', single; none uncoded.
SHORT's KILL: no unread pair torn. Its CONTROLS: the L = 1 classes
multiple at (3, 2) and (4, 3); at (6, 5), a pair F3 tears, the
L = 1 graphs below z''s.
PREIMAGE's KILL: no unread pair torn. Its CONTROL: (3, 5), torn by
F3, torn again here.
O1 Mult_d(m) is a union of M-orbits with one count per orbit at
every census arm. KILL: a class of Mult whose M-image is not in
Mult, or an orbit with two counts. CONTROLS: every class labelled,
its PLANT the census at (3, 4) rerun at a node limit of 8;
5 in T_3(11), E4's tear; 3 in T_3(4). T_d(m) printed, no
prediction.
NEG at m' = 4 (frozen after the census printed 3 in T_d(4) at d =
3, 4, 5): -1 in T_d(4) at d = 6, 9, 10, read by --neg. KILL at d:
none among the first 40 multiply coded w, or none in the first
1500 classes. CONTROL first: d = 3, 4, 5.
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). C2 holds at d = 2..10: 0's
live part the node 0 with its loop, -alpha's cycles one pattern.
- C1 holds at d = 2..8: the 35 graphs' live cycles are 30 zero loops
and 5 patterns, none other, no tangle. C3 sees one cycle 'other' at
every d = 2..10, 001, 0011, ..., rotations of 0 1^(d-1) 0, as
predicted; the cycle PLANT reads ['other', 'pattern', 'tangle',
'zero'] and ['tangle'].
- 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.
- CENSUS: KILL of H1 FIRES. 371 of 12675 nonzero classes multiply
coded, counted per (d, m), so a point of P_m' recurs in every P_m,
m' | m; by d {2: 0, 3: 28, 4: 54, 5: 289}. Each has 2, 3, 4, 5, 6
or 10 codings (live walks settled at 100 and 200; 292, 71, 2, 1, 2,
3 classes), no branch node, one live cycle but at 5 classes with
two. Nothing unread or uncoded; the controls hold at every pair.
- SHORT: the controls as predicted, the L = 1 graphs at (6, 5) of 72,
66, 52 and 65 nodes against z''s 6306. TORN, so SHORT's KILL
misses: (7, 5), (7, 11), (8, 7), (8, 11), (10, 3), (10, 5),
(10, 7), orbits 12, 10, 240, 120, 8, 24, 38, graphs up to 433690
nodes. --short at (6, 11): no class of orbit L <= 2000; at (9, 7)
only L = 6's six classes, multiple.
- PREIMAGE: TORN, so its KILL misses: (7, 7) by m' = 11, (8, 5),
(9, 5), (9, 11) by m' = 3, with the control (3, 5) and (7, 5), n z'
in w's class below 1e-15 at each. --preimage at (6, 11) (m' = 3, 7)
and (9, 7) (m' = 5, 11) reads no multiply coded short-orbit w: those
pairs are untested by it, not refuted.
- O1: KILL misses. At every census arm Mult_d(m) is a union of
M-orbits, one count per orbit, every class labelled (the PLANT, 48
classes unread at a limit of 8, caught): orbits of 1 to 12, at (3, 11)
one orbit of 10 with 2 codings. T_d(m) is the non-units of Z/m
alone at every pair but (3, 4), (3, 8), (3, 11), (3, 12), (4, 4),
(4, 6), (5, 4), whose units in T are [3], [3, 7], every residue
but 1, [7, 11], [3], [5], [3]: so at d = 3 every x n with n != 1
mod 11 is torn. Its controls 3 in T_3(4) and 5 in T_3(11) hold.
- NEG: -w single and w multiple at m' = 4, w = (1, 3, 3) at d = 3
(the 31st class read by --neg) and (0, ..., 0, 3, 0, 3) at d = 4,
5, 6, 9, 10 (the 51st), n z' in w's class below 1e-15: every n = 3
mod 4 is torn at those d. The table reads w = -sigma = (2, 1, 2,
0, ..., 0) at d = 4..10, the class periodic.py proves multiply
coded at every d >= 4: multiple, -w = sigma single, n z' in w's
class below 1e-15 at each.
- CLOSURE: torn 99 of 99; by E1 alone (7, 10), (8, 10), (9, 10),
(9, 12), (10, 6), (10, 9), (10, 10), (10, 12); by E4 alone (6, 11),
(7, 7), (8, 5), (9, 5), (9, 7), (9, 11), (10, 11); E4's families
every n = 3 mod 4 at d = 3..10, n = 5 mod 11 at d = 3
and 7, n = 7 mod 11 at d = 7, n = 2 mod 3 at d = 8 and 9.
- Tiers: the count, the lemma and 0's one coding are proved (the
arguments above, (F) and Parry cited); P1, N and C1-C2 are finite
checks of them and kill nothing; the tears are a finite check per
pair, in double precision with the chosen margins, and so are the
classes' labels; E1 and E4 are proved (THE CLASSES).
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 178/178, 18.0 s, peak commit 330.4 MB (the census
plant and O1's kill as a check added after the --full); --full
287/287, 597.9 s, peak commit 498.3 MB, under a 600 MB memory guard;
--short and --preimage at (6, 11), (9, 7) 3.8 s and 27.9 s; --neg
3, 4, 5 0.3 s; at d = 6, 9, 10 the same search ran in the lab before
its port here (584.5 s, d = 10's 51 reads 532 s), not rerun. With
the -sigma rows: --full 290/290, 640.8 s, peak commit 500.0 MB under
the 600 MB guard.
"""
import cmath
import itertools
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, cycles=False):
"""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; with
cycles, the live part's cycles as (kind, word)."""
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
nn = len(succ)
if not cycles:
del nodes
live = peel(succ)
cyc = [(kind(c, d), c[1]) for c in cycle_words(nodes, succ, live)] \
if cycles else None
if cycles:
del nodes
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)
# from each live start its least live path on 400 digits: with no
# branch node, that start's one coding
codings = []
for u in ls:
y = []
for _ in range(400):
t = [(x, t) for x, t in enumerate(succ[u]) if t >= 0 and live[t]]
y.append(t[0][0])
u = t[0][1]
codings.append(y)
coding = codings[0] if codings else None
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, codings=codings, walks=walks,
cycles=cyc)
def peel(succ):
"""live: the 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)
return live
def sccs(succ, live):
"""Tarjan, iterative, on the live subgraph."""
nn = len(succ)
index, low, on, st, out = [None] * nn, [0] * nn, [False] * nn, [], []
cnt = 0
for v0 in range(nn):
if not live[v0] or index[v0] is not None:
continue
work = [(v0, 0)]
while work:
v, pi = work.pop()
if pi == 0:
index[v] = low[v] = cnt
cnt += 1
st.append(v)
on[v] = True
nxt = [t for t in succ[v] if t >= 0 and live[t]]
if pi < len(nxt):
work.append((v, pi + 1))
t = nxt[pi]
if index[t] is None:
work.append((t, 0))
elif on[t]:
low[v] = min(low[v], index[t])
continue
for t in nxt:
if on[t] and index[t] > index[v]:
low[v] = min(low[v], low[t])
if low[v] == index[v]:
comp = []
while True:
w = st.pop()
on[w] = False
comp.append(w)
if w == v:
break
out.append(comp)
return out
def cycle_words(nodes, succ, live):
"""Each live SCC holding a cycle as ('cycle', word, state), read once
around from its least node, or ('tangle', size, None)."""
res = []
for comp in sccs(succ, live):
cs = set(comp)
inner = {u: [(y, t) for y, t in enumerate(succ[u]) if t in cs]
for u in comp}
if len(comp) == 1 and not inner[comp[0]]:
continue # no cycle
if any(len(e) != 1 for e in inner.values()):
res.append(("tangle", len(comp), None))
continue
u0 = min(comp, key=lambda u: nodes[u])
word, u = [], u0
while True:
y, u = inner[u][0]
word.append(y)
if u == u0:
break
res.append(("cycle", "".join(map(str, word)), nodes[u0][0]))
return res
def kind(c, d):
"""zero (the loop 0 at state 0), pattern (a rotation of 1^(d-1)0),
other, or tangle."""
if c[0] != "cycle":
return "tangle"
w = c[1]
if w == "0" and not any(c[2]):
return "zero"
pat = "1" * (d - 1) + "0"
if len(w) == d and w in pat + pat:
return "pattern"
return "other"
# ------------------------------------------------------- the classes of P_m
def label(f):
if f is None:
return "unread"
if f["live_starts"] == 0:
return "uncoded"
return "single" if f["single"] else "multiple"
def read_class(N, m, key, **kw):
"""point_fibre at the class key = (r_0,) + g, representative r = r_0
or m."""
return point_fibre(N, m, key[0] or m, list(key[1:]), **kw)
def census(N, m, limit=800_000):
"""Every nonzero class of P_m by its label."""
d = N.dim
out = {"single": [], "multiple": [], "unread": [], "uncoded": []}
for key in itertools.product(range(m), repeat=d):
if any(key):
f = read_class(N, m, key, limit=limit)
out[label(f)].append((key, f))
return out
def rref(A, p):
"""Row-reduce A over F_p, p prime: (rows, pivot columns)."""
A, piv, r = [list(row) for row in A], [], 0
for c in range(len(A[0]) if A else 0):
k = next((i for i in range(r, len(A)) if A[i][c] % p), None)
if k is None:
continue
A[r], A[k] = A[k], A[r]
inv = pow(A[r][c], p - 2, p)
A[r] = [x * inv % p for x in A[r]]
for i in range(len(A)):
if i != r and A[i][c]:
f = A[i][c]
A[i] = [(x - f * y) % p for x, y in zip(A[i], A[r])]
piv.append(c)
r += 1
return A, piv
def kernel(A, p):
"""A basis of the kernel of the square A over F_p."""
n = len(A)
R, piv = rref(A, p)
basis = []
for f in (c for c in range(n) if c not in piv):
v = [0] * n
v[f] = 1
for i, c in enumerate(piv):
v[c] = -R[i][f] % p
basis.append(v)
return basis
def mat_mod(A, B, p):
n = len(A)
return [[sum(A[i][k] * B[k][j] for k in range(n)) % p for j in range(n)]
for i in range(n)]
def orbit(N, p, key):
"""The length of key's orbit under M mod p, on s_0 mod p = -key."""
c = s = tuple(-x % p for x in key)
for L in range(1, p ** N.dim + 1):
s = tuple(sum(N.M[i][j] * s[j] for j in range(N.dim)) % p
for i in range(N.dim))
if s == c:
return L
def short_orbits(N, p, count=3, top=2000, cap=60):
"""The first count L <= top whose kernel of M^L - I over F_p leaves
the span of the earlier ones (a multiple of L* repeats L*'s kernel),
each with up to cap of its kernel vectors outside that span."""
d = N.dim
rank = lambda V: len(rref(V, p)[1])
P, out, span = [row[:] for row in N.M], [], []
for L in range(1, top + 1):
B = kernel([[(P[i][j] - (i == j)) % p for j in range(d)]
for i in range(d)], p)
if B and rank(span + B) > rank(span):
vs = []
for co in itertools.product(range(p), repeat=len(B)):
v = tuple(sum(c * b[i] for c, b in zip(co, B)) % p
for i in range(d))
if rank(span + [list(v)]) > rank(span):
vs.append(v)
if len(vs) == cap:
break
out.append((L, len(B), vs))
span = span + B
if len(out) == count:
break
P = mat_mod(P, N.M, p)
return out
def short_keys(N, p):
"""The classes of short_orbits as keys, s_0 mod p = c giving r_0 =
-c_1, g = -c_rest, with their L."""
return [(L, tuple(-x % p for x in c)) for L, _, vs in short_orbits(N, p)
for c in vs]
def torus(N, m, key):
"""The point (g - r alpha)/m mod Z^(d-1) of key at r = r_0 or m."""
r = key[0] or m
al = [N.beta ** -(i + 1) for i in range(N.dim - 1)]
return [((gi - r * a) / m) % 1 for gi, a in zip(key[1:], al)]
def scaled_off(N, mp, n, z, w):
"""|n z' - (w - j alpha)| mod Z^(d-1), j = (n r' - r_w)/m': n z' lies
in w's class."""
j = (n * (z[0] or mp) - (w[0] or mp)) // mp
al = [N.beta ** -(i + 1) for i in range(N.dim - 1)]
nz = [(n * x) % 1 for x in torus(N, mp, z)]
wj = [(x - j * a) % 1 for x, a in zip(torus(N, mp, w), al)]
return max(min(abs(u - v), 1 - abs(u - v)) for u, v in zip(nz, wj))
def act(N, m, key):
"""M key mod m: the class (r, g) to (r + sum g, r, g_1, ...)."""
return tuple(sum(N.M[i][j] * key[j] for j in range(N.dim)) % m
for i in range(N.dim))
def orbits(N, m, S):
"""S cut into M-orbits, or None when some class of S leaves S."""
if any(act(N, m, k) not in S for k in S):
return None
left, out = set(S), []
while left:
orb = [min(left)]
while act(N, m, orb[-1]) != orb[0]:
orb.append(act(N, m, orb[-1]))
left -= set(orb)
out.append(orb)
return out
def tears(m, single, mult):
"""T_d(m): the residues n with n kappa multiply coded or zero for some
singly coded kappa."""
return [n for n in range(m) if any(
tuple(n * x % m for x in k) in mult or not any(n * x % m for x in k)
for k in single)]
def o1(N, m, out, cnt):
"""O1 at one census: every class labelled, Mult a union of M-orbits,
one count per orbit; and T_d(m)."""
single = [k for k, _ in out["single"]]
mult = {k for k, _ in out["multiple"]}
orb = orbits(N, m, mult)
return dict(labelled=len(single) + len(mult) == m ** N.dim - 1,
orbits=orb, mixed=orb is not None and
any(len({cnt[k] for k in o}) > 1 for o in orb),
T=tears(m, single, mult))
# ------------------------------------------------------------ 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), cycles=True)
z = point_fibre(N, 1, 0, [0] * (d - 1), cycles=True)
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")
check(f"C2: d = {d}: 0's live part is one node with the loop 0; "
f"-alpha's cycles {sorted(k for k, _ in a['cycles'])}",
z["live"] == 1 and z["cycles"] == [("zero", "0")] and
set(k for k, _ in a["cycles"]) <= {"zero", "pattern"},
f"0's cycles {z['cycles']}")
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)
fs = [point_fibre(N, 1, -n, [0] * (d - 1), cycles=True)
for n in range(1, 31)]
bad1 = sum(1 for n, f in enumerate(fs, 1) if not f["single"]
or f["coding"][:60] != N.greedy(n, 60))
gs = [point_fibre(N, 1, k, [0] * (d - 1), walk=400, cycles=True)
for k in range(1, 6)]
ws = [f["walks"] for f in gs]
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}")
seen = {}
for f in fs + gs:
for k, _ in f["cycles"]:
seen[k] = seen.get(k, 0) + 1
check(f"C1: d = {d}: every live cycle of the 35 graphs is zero or a "
f"rotation of 1^({d - 1})0", set(seen) <= {"zero", "pattern"},
f"cycles by kind {seen}")
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}")
# the cycle reader at d = 4: a loop 0 at state 0; 1 -> 2 -> 3 -> 4 -> 1
# reading 1100; 5 -> .. -> 5 reading 1101; 9 with a loop 0 and a loop 1
# through 10; and, apart, a tangle at state 0 whose first edges read 0
nodes = [((i, 0), 0) for i in range(11)]
succ = [[0, -1], [-1, 2], [-1, 3], [4, -1], [1, -1],
[-1, 6], [-1, 7], [8, -1], [-1, 5], [9, 10], [9, -1]]
ks = sorted(kind(c, 4) for c in cycle_words(nodes, succ, peel(succ)))
s2 = [[0, 1], [0, -1]]
k2 = [kind(c, 4) for c in cycle_words([((0, 0), 0), ((1, 0), 1)], s2,
peel(s2))]
check("PLANT: 1100 reads 'other', 1101 'pattern', the loop at 0 "
"'zero', the two loops at 9 and at 0 'tangle'",
ks == ["other", "pattern", "tangle", "zero"] and k2 == ["tangle"],
f"{ks}, {k2}")
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), cycles=m == 2)
if m == 2:
want = "0" + "1" * (d - 1) + "0"
check(f"C3: d = {d}, x2: a live cycle neither zero nor a "
f"pattern, a rotation of {want}", f is not None and
any(k == "other" and len(w) == d + 1 and w in want * 2
for k, w in f["cycles"]),
f"cycles {f and f['cycles']}")
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, cert
ARMS = [(2, range(2, 13)), (3, range(2, 13)), (4, range(2, 8)),
(5, range(2, 5))]
# singly coded short-orbit classes, found by --short
SHORT = {(7, 5): (1, 4, 3, 3, 2, 4, 4), (7, 11): (9, 6, 9, 9, 2, 0, 10),
(8, 7): (2, 0, 5, 0, 0, 0, 0, 6), (8, 11): (5, 4, 6, 7, 1, 0, 0, 10),
(10, 3): (0, 2, 0, 1, 1, 2, 1, 0, 0, 2),
(10, 5): (3, 2, 4, 0, 1, 4, 1, 0, 0, 1),
(10, 7): (6, 1, 6, 2, 1, 5, 3, 0, 0, 4)}
# (d, n, m', w): w multiply coded in P_m', z' = n^(-1) w, found by --preimage
PREIMAGE = [(3, 5, 11, (8, 0, 10)), (7, 5, 11, (5, 7, 5, 5, 6, 0, 8)),
(7, 7, 11, (5, 7, 5, 5, 6, 0, 8)),
(8, 5, 3, (0, 0, 0, 0, 2, 0, 2, 1)),
(9, 5, 3, (0, 0, 0, 0, 0, 2, 2, 1, 2)),
(9, 11, 3, (0, 0, 0, 0, 0, 2, 2, 1, 2))]
# (d, 3, 4, w), 3^(-1) w = -w singly coded: at d = 3 the first such w of
# P_4 in itertools order, found by --neg; at d >= 4 w = -sigma
PREIMAGE += [(3, 3, 4, (1, 3, 3))] + [
(d, 3, 4, (2, 1, 2) + (0,) * (d - 3)) for d in range(4, 11)]
def section_census(arms):
print("CENSUS every nonzero class of P_m", flush=True)
total = {"single": 0, "multiple": 0, "unread": 0, "uncoded": 0}
mult, kill = {}, []
for d, ms in arms:
N = Num(d)
for m in ms:
t0 = time.time()
z = point_fibre(N, m, m, [0] * (d - 1), walk=400)
out = census(N, m)
for k in total:
total[k] += len(out[k])
mult[d] = mult.get(d, 0) + len(out["multiple"])
print(f" d = {d}, x{m}: {m ** d - 1} classes, single "
f"{len(out['single'])}, multiple {len(out['multiple'])}, "
f"unread {len(out['unread'])}, uncoded "
f"{len(out['uncoded'])} [{time.time() - t0:.1f} s]",
flush=True)
cnt = {}
for c, f in out["multiple"]:
f2 = read_class(N, m, c, walk=200, cycles=True)
cnt[c] = f2["walks"]
print(f" MULTIPLE {c}: branch {f['branch']}, walks at "
f"100, 200 {f2['walks']}, cycles {f2['cycles']}",
flush=True)
zp = [k for k, lst in out.items()
if any(c == (1, 1) + (0,) * (d - 2) for c, _ in lst)]
check(f"CONTROL d = {d}, x{m}: -alpha at (m, 0) {d} codings, "
f"labelled multiple; z' single; none uncoded",
z is not None and z["walks"] == (d, d) and
label(z) == "multiple" and zp == ["single"] and
not out["uncoded"], f"walks {z and z['walks']}, z' {zp}")
o = o1(N, m, out, cnt)
check(f"CONTROL d = {d}, x{m}: every class labelled",
o["labelled"])
if o["orbits"] is None or o["mixed"]:
kill.append((d, m))
units = [n for n in o["T"] if math.gcd(n, m) == 1]
sizes = o["orbits"] and sorted((len(b), cnt[b[0]][1])
for b in o["orbits"])
print(f" O1: orbits (size, count) {sizes}, mixed "
f"{o['mixed']}; T = {o['T']}, units {units}", flush=True)
if (d, m) == (3, 4):
check("CONTROL: 3 in T_3(4)", 3 in o["T"])
pl = census(N, m, limit=8)
check(f"PLANT d = 3, x4: the census at a node limit of 8 "
f"leaves {len(pl['unread'])} unread, not all labelled",
pl["unread"] and not o1(N, m, pl, cnt)["labelled"])
if (d, m) == (3, 11):
check("CONTROL: 5 in T_3(11), E4's tear", 5 in o["T"])
print(f" census: {total}; multiple by d {mult}", flush=True)
check("KILL of O1 misses: every Mult a union of M-orbits, one count "
"per orbit", not kill, f"fires at {kill}")
print(f" KILL of O1: {'FIRES at ' + str(kill) if kill else 'misses'}",
flush=True)
return mult
def section_short(ds):
print("SHORT singly coded classes of short orbit under M mod p",
flush=True)
for d, p in [(3, 2), (4, 3)]:
N = Num(d)
fix = [c for L, c in short_keys(N, p) if L == 1]
labs = [label(read_class(N, p, c)) for c in fix]
check(f"CONTROL d = {d}, x{p}: the L = 1 classes {fix} multiple",
fix and set(labs) == {"multiple"}, f"{labs}")
if 6 in ds:
N = Num(6)
zn = point_fibre(N, 5, 1, [1] + [0] * 4)["nodes"]
ns = [read_class(N, 5, c)["nodes"] for L, c in short_keys(N, 5)
if L == 1]
check("CONTROL d = 6, x5: the L = 1 graphs below z''s",
ns and max(ns) < zn, f"{ns} against {zn} nodes")
torn = []
for (d, p), key in SHORT.items():
if d not in ds:
continue
t0 = time.time()
f = read_class(Num(d), p, key)
if label(f) == "single":
torn.append((d, p))
check(f"d = {d}, x{p}: class {key}, orbit {orbit(Num(d), p, key)}, "
f"single", label(f) == "single", f"{label(f)}, nodes "
f"{f and f['nodes']} [{time.time() - t0:.1f} s]")
return torn
def section_preimage(ds):
print("PREIMAGE z' = n^(-1) w for a multiply coded class w of P_m'",
flush=True)
fam = []
for d, n, mp, w in PREIMAGE:
if d not in ds:
continue
N, t0 = Num(d), time.time()
z = tuple(pow(n, -1, mp) * x % mp for x in w)
lw, lz = label(read_class(N, mp, w)), label(read_class(N, mp, z))
off = scaled_off(N, mp, n, z, w)
check(f"d = {d}, x{n}, m' = {mp}: w {w} multiple, z' {z} single, "
f"n z' in w's class", lw == "multiple" and lz == "single" and
off < 1e-9, f"{lw}, {lz}, off {off:.1e} "
f"[{time.time() - t0:.1f} s]")
if lw == "multiple" and lz == "single":
fam.append((d, n % mp, mp))
return fam
def section_closure(ds, ms, f3, cert, short, fam):
"""E1: a point tear z, mz = -r alpha, tears x km; E4: a preimage tears
every n = n_0 mod m'."""
print("CLOSURE the pairs torn by a point, the certificate or E1, E4")
base = {(d, m) for (d, m), v in f3.items() if v == "TORN"} | set(short)
how = {}
for d in ds:
for m in ms:
h = [w for w, ok in [
("F3", f3.get((d, m)) == "TORN"), ("cert", (d, m) in cert),
("short", (d, m) in short),
("E1", any((d, b) in base for b in ms if b < m and
m % b == 0)),
("E4", any(e == d and m % mp == u for e, u, mp in fam))] if ok]
how[d, m] = h
torn = [k for k, h in how.items() if h]
print(f" torn: {len(torn)} of {len(how)}; open: "
f"{sorted(k for k, h in how.items() if not h)}")
for key in ("cert", "short", "E1", "E4"):
print(f" {key} alone: "
f"{sorted(k for k, h in how.items() if h == [key])}")
for e, u, mp in sorted(set(fam)):
print(f" E4: every n = {u} mod {mp} torn at d = {e}")
def short_search(d, p):
"""Arm 2: the shortest orbits' classes, read to the first single."""
N, t0 = Num(d), time.time()
seen = {}
for L, c in short_keys(N, p):
if c not in seen:
seen[c] = label(read_class(N, p, c))
if seen[c] == "single":
print(f" d = {d}, x{p}: TORN at class {c}, L = {L}, read "
f"{len(seen)} [{time.time() - t0:.1f} s]", flush=True)
return c
print(f" d = {d}, x{p}: open, read {len(seen)}: "
f"{sorted(set(seen.values()))} [{time.time() - t0:.1f} s]")
def preimage_search(d, n, mps=(2, 3, 5, 7, 11)):
"""Arm 3: z' = n^(-1) w over the multiply coded short-orbit w of P_m',
to the first single."""
N = Num(d)
for mp in mps:
if n % mp in (0, 1):
continue
t0 = time.time()
inv, ws, unread = pow(n, -1, mp), 0, 0
for _, w in dict.fromkeys(short_keys(N, mp)):
lw = label(read_class(N, mp, w))
unread += lw == "unread"
if lw != "multiple":
continue
ws += 1
z = tuple(inv * x % mp for x in w)
if label(read_class(N, mp, z)) == "single":
print(f" d = {d}, x{n}: TORN, m' = {mp}, w {w}, z' {z}, "
f"off {scaled_off(N, mp, n, z, w):.1e} "
f"[{time.time() - t0:.1f} s]", flush=True)
return mp, w
print(f" d = {d}, x{n}, m' = {mp}: none of {ws} multiply coded w, "
f"{unread} w unread [{time.time() - t0:.1f} s]", flush=True)
def neg_search(d, m=4, top_w=40, top=1500):
"""The first multiply coded w of P_m in itertools order whose -w is
singly coded."""
N, t0 = Num(d), time.time()
seen = ws = unread = 0
for key in itertools.product(range(m), repeat=d):
if not any(key):
continue
seen += 1
lab = label(read_class(N, m, key))
unread += lab == "unread"
if lab == "multiple":
ws += 1
f = read_class(N, m, tuple(-x % m for x in key))
if label(f) == "single":
print(f" d = {d}, x{m}: TORN by -1, w {key}, read {seen}, "
f"unread {unread} [{time.time() - t0:.1f} s]",
flush=True)
return key
if ws == top_w or seen == top:
break
print(f" d = {d}, x{m}: open, read {seen}, multiple {ws}, unread "
f"{unread} [{time.time() - t0:.1f} s]", flush=True)
def pairs(flag):
return [tuple(map(int, x.split(":"))) for x in
sys.argv[sys.argv.index(flag) + 1].split(",")]
def main():
if "--short" in sys.argv:
for d, p in pairs("--short"):
short_search(d, p)
return
if "--preimage" in sys.argv:
for d, n in pairs("--preimage"):
preimage_search(d, n)
return
if "--neg" in sys.argv:
for d in sys.argv[sys.argv.index("--neg") + 1].split(","):
neg_search(int(d))
return
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()
f3, cert = section_f3(ds, range(2, 13))
mult = section_census(ARMS if full else [(2, range(2, 9)),
(3, range(2, 9))])
check("CENSUS: nonzero classes multiply coded at d = 3, none at d = 2",
mult[2] == 0 and mult[3] > 0, f"{mult}")
short = section_short(ds)
fam = section_preimage(ds)
section_closure(ds, range(2, 13), f3, cert, short, fam)
print(f"{sum(CHECKS)}/{len(CHECKS)} checks passed")
sys.exit(0 if all(CHECKS) else 1)
if __name__ == "__main__":
main()