#!/usr/bin/env python3
"""check_proof.py: checks every claim of proof.md for W5-3 = 1RB1RA_1LC1LE_1RF1LD_1RE0LD_0LB0RA_---0RC.

Usage:  python3 -B check_proof.py [--seed S] [--quick | --full]

Self-contained (own simulators, no project imports).  Sections (proof.md section in brackets):
  T    transition table facts                                   [1, 6]
  L0   low-level rules: exact windowed runs, any context         [3.1]
  CQ   counter words C(n), Q_k(x,y): the string identities      [2, 3.3]
  L1   count, eat and Y rules: exact windowed runs, any context  [3.2, 3.3]
  X    the rule-level executor (L1 rules by pattern) = direct TM, G-configuration by G-configuration
  W    walk lemmas W, W2, W0 and the block stage                 [3.4]
  M    eat-to-eat macro rules: chain of proof.md = rule-level executor = direct TM; closed form;
       no eat configuration at any TM step strictly inside a macro rule                [3.5]
  INV  invariant J: coverage and closure                         [4]
  Z    initial segment (G form at 165, H form at 171), orbit from blank vs. stored and fresh direct runs [6, 8]
  NEG  the 6 mod 8 precondition is needed (the chain breaks below it)
  TEETH every single-transition mutation of the table is caught by the windowed L1 checks

A windowed run starts from a finite window of cells; every cell outside it is unknown and reading one is
a failure.  A rule checked this way therefore holds for every left and right context.
"""
import sys, os, re, random, time, argparse

TM = "1RB1RA_1LC1LE_1RF1LD_1RE0LD_0LB0RA_---0RC"
ST = "ABCDEF"
HERE = os.path.dirname(os.path.abspath(__file__))


def parse_tm(tm):
    T = {}
    for i, g in enumerate(tm.split('_')):
        for s in (0, 1):
            t = g[3 * s:3 * s + 3]
            T[(ST[i], s)] = None if t == '---' else (int(t[0]), 1 if t[1] == 'R' else -1, t[2])
    return T


T = parse_tm(TM)

# ---------------------------------------------------------------- bookkeeping
RES = {}
FAILS = []


def rec(sec, name, ok, info=''):
    k = (sec, name)
    r = RES.setdefault(k, [0, 0])
    r[0 if ok else 1] += 1
    if not ok and len(FAILS) < 60:
        FAILS.append(f"FAIL {sec}.{name}: {info}")


# ---------------------------------------------------------------- windowed runs
class OutOfWindow(Exception):
    pass


class Halted(Exception):
    pass


def win(spec):
    """spec: tokens separated by spaces. Cell tokens: strings over 01 (several cells per token allowed,
    empty tokens ignored). Head tokens: 'X>' (state X reads the next cell), '<X' (state X reads the
    previous cell). Returns (cells, head index, state); the head may sit just outside the window."""
    cells, head, st = [], None, None
    for tok in spec.split():
        m = re.fullmatch(r'([A-F])>', tok)
        if m:
            head, st = len(cells), m.group(1)
            continue
        m = re.fullmatch(r'<([A-F])', tok)
        if m:
            head, st = len(cells) - 1, m.group(1)
            continue
        assert re.fullmatch(r'[01]+', tok), tok
        cells.extend(int(c) for c in tok)
    assert head is not None, spec
    return tuple(cells), head, st


def run_window(cells, head, st, steps):
    c = list(cells)
    h = head
    tr = []
    for i in range(steps):
        if h < 0 or h >= len(c):
            raise OutOfWindow(f"reads outside the window at step {i}")
        r = c[h]
        t = T[(st, r)]
        if t is None:
            raise Halted(f"halts at step {i}")
        tr.append(st + str(r))
        c[h] = t[0]
        h += t[1]
        st = t[2]
    return tuple(c), h, st, tr


def check_rule(sec, name, lhs, rhs, steps, trans=None, info=''):
    L = win(lhs)
    R = win(rhs)
    if len(L[0]) != len(R[0]):
        rec(sec, name, False, f"window sizes differ {info}: {lhs} | {rhs}")
        return None
    try:
        c, h, s, tr = run_window(L[0], L[1], L[2], steps)
    except (OutOfWindow, Halted) as e:
        rec(sec, name, False, f"{e} {info}: {lhs}")
        return None
    ok = (c, h, s) == R
    if trans is not None:
        ok = ok and ' '.join(tr) == trans
    if 'F0' in tr:
        ok = False
    rec(sec, name, ok, f"{info}: got {c},{h},{s} want {R}; {' '.join(tr)}")
    return tr


def P(w, n):
    return w * n


# ---------------------------------------------------------------- counter words
def t1(n):
    """number of trailing 1 bits"""
    k = 0
    while n & 1:
        n >>= 1
        k += 1
    return k


def Cw(n):
    """C(n) as a finite string (followed by 0^inf): digits of n below its leading 1, lowest first,
    1 -> 11, 0 -> 00, then 11.  C(0) = '' (the blank half-tape)."""
    if n == 0:
        return ''
    d = n.bit_length() - 1
    return ''.join('11' if (n >> i) & 1 else '00' for i in range(d)) + '11'


def Dk(k, x):
    """k digits of x (0 <= x < 2^k), lowest first"""
    return ''.join('11' if (x >> i) & 1 else '00' for i in range(k))


def Qw(k, x, y):
    """Q_k(x, y) = D_k(x) 0 C(y)"""
    return Dk(k, x) + '0' + Cw(y)


def norm_R(R):
    return R.rstrip('0')


def norm_w(w):
    return w.lstrip('0')


def is_C(R):
    """R (normalized) == C(n) for some n >= 0 ?  returns n or None"""
    R = norm_R(R)
    if R == '':
        return 0
    if len(R) % 2 or not R.endswith('11'):
        return None
    n = 0
    pairs = [R[i:i + 2] for i in range(0, len(R), 2)]
    if any(q not in ('00', '11') for q in pairs):
        return None
    for i, q in enumerate(pairs[:-1]):
        if q == '11':
            n |= 1 << i
    return n | (1 << (len(pairs) - 1))


def lead11(R):
    """number of leading '11' pairs of R (R padded with zeros)"""
    k = 0
    while R[2 * k:2 * k + 2] == '11':
        k += 1
    return k


def incR(R):
    """R+ : R = (11)^k 00 R'  ->  (00)^k 11 R'.  Returns (R+, k) or None."""
    Rp = R + '0' * 8
    k = lead11(Rp)
    if Rp[2 * k:2 * k + 2] != '00':
        return None
    return norm_R('00' * k + '11' + Rp[2 * k + 2:]), k


# ---------------------------------------------------------------- the rule-level executor (proof.md 3.2-3.3)
def left_part(w):
    """left processing of a G-configuration G(w, R). Returns (kind, v, b, steps) where after the left part
    the machine is at  v 1 0 A> 1 1 R  (kinds INC1, INC0, INC011) or at  v 1 1 0 A> R  (kind EAT)."""
    W = '0' * 8 + w
    b = 0
    while W.endswith('10'):
        W = W[:-2]
        b += 1
    if W.endswith('1'):
        return 'INC1', norm_w(W[:-1] + '0' + '10' * b), b, 6 * b + 11
    if W.endswith('0100'):
        return 'INC0', norm_w(W[:-4] + '1010' + '10' * b), b, 6 * b + 19
    if W.endswith('101100'):
        return 'INC011', norm_w(W[:-6] + '010110' + '10' * b), b, 6 * b + 27
    if W.endswith('000'):
        return 'EAT', norm_w(W[:-3] + '10' * (b + 2)), b, 4 * b + 17
    return None


def gstep(w, R):
    """one G -> G transition by the L1 rules. Returns (rulename, w', R', steps) or None if no rule applies."""
    lp = left_part(w)
    if lp is None:
        return None
    kind, v, b, s0 = lp
    Rp = R + '0' * 12
    if kind != 'EAT':
        ir = incR(R)
        if ir is not None:
            Rn, k = ir
            return (kind, v, Rn, s0 + 4 * k + 13)
        # overflow: R = (11)^k 0 1 1 R''
        k = lead11(Rp)
        if Rp[2 * k:2 * k + 3] == '011':
            ir2 = incR(Rp[2 * k + 3:])
            if ir2 is None:
                return None
            Rn, j = ir2
            name = 'OVF' if kind == 'INC1' else kind + '+OVF'
            return (name, norm_w(v + '10' + '1' * (2 * k + 1)), Rn, s0 + 2 * k + 4 * j + 18)
        return None
    # Y processing on  v 1 1 0 A> R
    if Rp.startswith('11'):
        s = lead11(Rp)
        if Rp[2 * s:2 * s + 2] == '00':
            return ('EAT+Y1', norm_w(v + '1'), norm_R('00' * (s - 1) + '11' + Rp[2 * s + 2:]), s0 + 4 * s + 9)
        if s == 1 and Rp[2:5] == '011':
            ir2 = incR(Rp[5:])
            if ir2 is None:
                return None
            Rn, j = ir2
            return ('EAT+Y4', norm_w(v + '1101'), Rn, s0 + 4 * j + 18)
        return None
    if Rp.startswith('0000'):
        return ('EAT+Y2', norm_w(v + '111'), norm_R(Rp[4:]), s0 + 7)
    if Rp.startswith('0011'):
        s = lead11(Rp[4:])
        if Rp[4 + 2 * s:6 + 2 * s] == '00':
            return ('EAT+Y3', norm_w(v + '11'), norm_R('0' * (2 * s + 1) + '11' + Rp[6 + 2 * s:]), s0 + 4 * s + 20)
        return None
    return None


def is_E(w, R):
    """E-configuration: w = (10)^p and R = C(N), N >= 1"""
    if not re.fullmatch(r'(10)*', w):
        return None
    n = is_C(R)
    if n is None or n == 0:
        return None
    return (len(w) // 2, n)


def exec_EE(p, N, maxg=10 ** 6):
    """rule-level execution from E(p,N) to the next E-configuration.
    Returns (p', N', [(rule, w, R, steps)...]) or (None, reason, chain)."""
    w, R = '10' * p, Cw(N)
    chain = []
    for i in range(maxg):
        g = gstep(w, R)
        if g is None:
            return None, f"no rule at G({w},{R})", chain
        name, w, R, st = g
        chain.append((name, w, R, st))
        e = is_E(w, R)
        if e is not None:
            return e[0], e[1], chain
    return None, 'too long', chain


# ---------------------------------------------------------------- the chain of proof.md 3.5 (closed forms + step formulas)
def walk(a, b, n, L):
    """W(a,b): G((10)^a 1 (10)^b, C(n)) -> E(a+b+[a odd], n+a+1) in a+1 counts"""
    while True:
        L.append(('INC1', 6 * b + 4 * t1(n) + 24))
        n += 1
        if a == 0:
            return b, n
        L.append(('INC0', 6 * b + 4 * t1(n) + 32))
        n += 1
        if a == 1:
            return b + 2, n
        a -= 2
        b += 2


def split6(N):
    """N = 2 mod 4  ->  (s, m) with N + 2 = (2m+1) 2^(s+2)"""
    v = N + 2
    s = 0
    while v % 2 == 0:
        v //= 2
        s += 1
    return s - 2, (v - 1) // 2


def chain_EE(p, N):
    """the eat-to-eat chain of proof.md 3.5 with the step formulas of 3.2-3.3.
    Returns (p', N', [(rule, steps), ...], case) or (None, reason, [], case)."""
    L = []
    if N % 2 == 1:
        L.append(('EAT+Y1', 4 * p + 17 + 4 * t1(N) + 9))
        pp, n = walk(p + 2, 0, (N + 1) // 2, L)
        return pp, n, L, 'odd'
    if N % 4 == 0:
        L.append(('EAT+Y2', 4 * p + 17 + 7))
        f = N // 4
        L.append(('INC1', 0 + 4 * t1(f) + 24))
        pp, n = walk(p + 2, 1, f + 1, L)
        return pp, n, L, '0mod4'
    s, m = split6(N)
    L.append(('EAT+Y3', 4 * p + 17 + 4 * s + 20))
    if s == 0:
        L.append(('OVF', 0 + 0 + 4 * t1(m) + 29))
        pp, n = walk(p + 4, 0, m + 1, L)
        return pp, n, L, '2mod8'
    L.append(('INC1', 24))  # Q_s(0) -> Q_s(1), b = 0, k = 0
    L.append(('EAT+' + ('Y4' if s == 1 else 'Y1'), 4 * (p + 3) + 17 + (4 * t1(m) + 18 if s == 1 else 13)))
    if s == 1:  # t = 0
        n = m + 1
        L.append(('INC1', 4 * t1(n) + 24)); n += 1
        L.append(('INC011', 4 * t1(n) + 40)); n += 1
        L.append(('INC1', 6 + 4 * t1(n) + 24)); n += 1
        L.append(('INC0', 6 + 4 * t1(n) + 32)); n += 1
        return p + 7, n, L, '6mod8,t=0'
    t = s - 1
    if p + 7 < 2 ** t:
        return None, 'precondition p+7 >= 2^t fails', [], f'6mod8,t={t}'
    a, b, x = p + 5, 0, 1
    for i in range(2 ** (t - 1) - 1):           # phase (d): W2 on Q_t
        L.append(('INC1', 6 * b + 4 * t1(x) + 24)); x += 1
        L.append(('INC0', 6 * b + 4 * t1(x) + 32)); x += 1
        a -= 2
        b += 2
    assert x == 2 ** t - 1 and a == p + 7 - 2 ** t and b == 2 ** t - 2
    L.append(('OVF', 6 * b + 2 * t + 4 * t1(m) + 29))   # phase (e)
    n = m + 1
    for j in range(t + 1):                       # phase (f): block stage
        L.append(('INC1', 6 * j + 4 * t1(n) + 24)); n += 1
    c, d = (2 ** t - 2) // 2, t                   # phase (g): walk stage W0
    for i in range(c):
        L.append(('INC0', 6 * (d + 2 * i) + 4 * t1(n) + 32)); n += 1
        L.append(('INC1', 6 * (d + 2 * i + 2) + 4 * t1(n) + 24)); n += 1
    L.append(('INC0', 6 * (d + 2 * c) + 4 * t1(n) + 32)); n += 1
    return p + 7 + t, n, L, f'6mod8,t={t}'


def model(p, N):
    """closed form of the eat-to-eat map (proof.md 1 and 3.5)"""
    po = p % 2
    if N % 2:
        return p + 2 + po, (N + 1) // 2 + p + 3
    if N % 4 == 0:
        return p + 3 + po, N // 4 + p + 4
    if N % 8 == 2:
        return p + 4 + po, (N + 6) // 8 + p + 5
    s, m = split6(N)
    if s == 1:
        return p + 7, m + 5
    t = s - 1
    return p + 7 + t, m + 2 ** t + t + 1


def macro_steps(p, N):
    r = chain_EE(p, N)
    return sum(s for _, s in r[2])


# ---------------------------------------------------------------- direct TM (fast)
TW = [0] * 12
TMV = [0] * 12
TN = [0] * 12
for (q, s), t in T.items():
    k = 2 * ST.index(q) + s
    if t is None:
        TN[k] = -1
    else:
        TW[k], TMV[k], TN[k] = t[0], t[1], 2 * ST.index(t[2])
HALTK = [k for k in range(12) if TN[k] < 0]
assert HALTK == [10]


def fast_run(tape, pos, st2, n):
    tw, tm, tn = TW, TMV, TN
    for _ in range(n):
        k = st2 + tape[pos]
        if k == 10:
            raise Halted('F0')
        tape[pos] = tw[k]
        pos += tm[k]
        st2 = tn[k]
    return pos, st2


TR01 = bytes.maketrans(b'01', b'\x00\x01')
TRback = bytes.maketrans(b'\x00\x01', b'01')


def gtape(w, R, extra):
    """tape for G(w,R) = 0^inf w 0 1 1 B> 0 R 0^inf; returns (bytearray, head pos, state*2)"""
    left = (w + '011').encode().translate(TR01)
    right = ('0' + R).encode().translate(TR01)
    pad = extra + 16
    tape = bytearray(pad) + bytearray(left) + bytearray(right) + bytearray(pad)
    return tape, pad + len(left), 2


def tm_config(tape, pos, st2):
    """config up to translation: (state, left string without leading zeros, right string without trailing zeros)"""
    lo = tape.find(1)
    hi = tape.rfind(1)
    if lo < 0:
        return (ST[st2 // 2], '', '')
    if pos < 8 or pos > len(tape) - 8:
        raise OutOfWindow('tape too small')
    left = bytes(tape[min(lo, pos):pos]).translate(TRback).decode().lstrip('0')
    right = bytes(tape[pos:max(hi + 1, pos)]).translate(TRback).decode().rstrip('0')
    return (ST[st2 // 2], left, right)


def gconfig(w, R):
    return ('B', (w + '011').lstrip('0'), ('0' + R).rstrip('0'))


def hconfig(w, R):
    """H(w, R) = 0^inf w <E 1 0 1 1 R  (G(w, R) after TS; w nonempty): E reads the last cell of w"""
    return ('E', w[:-1].lstrip('0'), (w[-1] + '1011' + R).rstrip('0'))


def e_form(tape, pos, st2):
    """(p, N) if the TM configuration is E(p, N) = 0^inf (10)^p 0 1 1 B> 0 C(N) with N >= 1, else None"""
    st, left, right = tm_config(tape, pos, st2)
    if st != 'B':
        return None
    if left == '11':
        p = 0
    elif re.fullmatch(r'(?:10)+011', left):
        p = (len(left) - 3) // 2
    else:
        return None
    if right[:1] == '1':
        return None
    n = is_C(right[1:])
    if not n:
        return None
    return (p, n)


def watch_E(tape, pos, st2, n):
    """run n steps one at a time; return (pos, st2, [(step, (p, N)) for every E-form at steps 1..n])"""
    tw, tm, tn = TW, TMV, TN
    found = []
    for i in range(1, n + 1):
        k = st2 + tape[pos]
        if k == 10:
            raise Halted('F0')
        tape[pos] = tw[k]
        pos += tm[k]
        st2 = tn[k]
        if st2 == 2 and tape[pos] == 0 and tape[pos - 1] == 1 and tape[pos - 2] == 1 and tape[pos - 3] == 0:
            e = e_form(tape, pos, st2)
            if e is not None:
                found.append((i, e))
    return pos, st2, found


# ================================================================ sections
def sec_T():
    und = [k for k, v in T.items() if v is None]
    rec('T', 'only_F0_undefined', und == [('F', 0)], str(und))
    intoF = [k for k, v in T.items() if v is not None and v[2] == 'F']
    rec('T', 'F_entered_only_by_C0', intoF == [('C', 0)], str(intoF))


def L0_rules():
    """(name, lhs, rhs, steps, transitions) families of proof.md 3.1"""
    R = []
    R.append(('TS', '0 1 1 B> 0', '<E 1 0 1 1', 6, 'B0 C1 D1 D0 E0 B1'))
    R.append(('X1', '1 <E', '0 A>', 1, 'E1'))
    R.append(('X0', '0 1 0 0 <E', '1 0 A> 1 0', 5, 'E0 B0 C1 D0 E1'))
    R.append(('A01', 'A> 0 1', '0 A> 1', 3, 'A0 B1 E1'))
    R.append(('A0', 'A> 0', '1 B>', 1, 'A0'))
    R.append(('CS', '1 0 C> 1', '1 1 0 A>', 3, 'C1 D0 E1'))
    R.append(('Z0', '0 A> 0 0', '1 0 1 A>', 6, 'A0 B0 C1 D0 E1 A1'))
    R.append(('F011', '1 0 1 1 0 0 <E', '0 1 0 1 1 A> 0', 14, 'E0 B0 C1 D1 D0 E0 B1 E1 A1 A0 B1 E1 A1 A1'))
    return R


def sec_L0(Nmax):
    for name, l, r, st, tr in L0_rules():
        check_rule('L0', name, l, r, st, tr)
    for n in range(0, Nmax + 1):
        check_rule('L0', 'EB(n)', P('10', n) + ' <E', '<E ' + P('10', n), 2 * n, ' '.join(['E0 B1'] * n), f'n={n}')
        check_rule('L0', 'AS(n)', 'A> ' + P('10', n) + ' 1', P('10', n) + ' A> 1', 4 * n, ' '.join(['A1 A0 B1 E1'] * n), f'n={n}')
        check_rule('L0', 'A1s(n)', 'A> ' + P('1', n), P('1', n) + ' A>', n, ' '.join(['A1'] * n), f'n={n}')
        check_rule('L0', 'D1s(n)', P('1', n) + ' <D', '<D ' + P('0', n), n, ' '.join(['D1'] * n), f'n={n}')
        check_rule('L0', 'CF(n)', 'C> ' + P('01', n), P('10', n) + ' C>', 2 * n, ' '.join(['C0 F1'] * n), f'n={n}')
        if n >= 1:
            check_rule('L0', 'DE(k)', '0 ' + P('1', n) + ' A> 0 0', '<E 1 ' + P('0', n) + ' 1 1', n + 6,
                       ' '.join(['A0 B0 C1'] + ['D1'] * n + ['D0 E0 B1']), f'k={n}')
        check_rule('L0', 'EAT0(n)', '0 0 0 <E ' + P('10', n) + ' 1', P('10', n + 2) + ' C>', 2 * n + 6,
                   ' '.join(['E0 B0 C0 F1'] + ['C0 F1'] * (n + 1)), f'n={n}')
        check_rule('L0', 'CARRY(k)', '1 0 1 1 ' + P('11', n) + ' A> 0 0', '0 1 1 B> 0 ' + P('00', n) + ' 1 1', 2 * n + 11,
                   ' '.join(['A0 B0 C1'] + ['D1'] * (2 * n + 2) + ['D0 E0 B1 E1 A1 A0']), f'k={n}')
    # compositions stated in 3.1: CARRY(k) = DE(2k+2), X1, A1, A0 ; DE(1) = A0 then TS
    for k in range(0, Nmax + 1):
        tr = check_rule('L0', 'CARRY=DE+X1+A1+A0', '1 0 1 1 ' + P('11', k) + ' A> 0 0', '0 1 1 B> 0 ' + P('00', k) + ' 1 1', 2 * k + 11)
        want = ['A0', 'B0', 'C1'] + ['D1'] * (2 * k + 2) + ['D0', 'E0', 'B1'] + ['E1', 'A1', 'A0']
        rec('L0', 'CARRY=DE+X1+A1+A0(trans)', tr == want, str(k))
    check_rule('L0', 'DE(1)=A0+TS', '0 1 A> 0 0', '<E 1 0 1 1', 7, 'A0 B0 C1 D1 D0 E0 B1')
    # C0 (the only way into F) occurs only in EAT0 and CF; every F step there reads 1
    users = set()
    for name, l, r, st, tr in L0_rules():
        if 'C0' in tr.split():
            users.add(name)
    for n in range(0, 6):
        fam = {'EB': ' '.join(['E0 B1'] * n), 'AS': ' '.join(['A1 A0 B1 E1'] * n), 'A1s': ' '.join(['A1'] * n),
               'D1s': ' '.join(['D1'] * n), 'CF': ' '.join(['C0 F1'] * n),
               'DE': ' '.join(['A0 B0 C1'] + ['D1'] * (n + 1) + ['D0 E0 B1']),
               'EAT0': ' '.join(['E0 B0 C0 F1'] + ['C0 F1'] * (n + 1)),
               'CARRY': ' '.join(['A0 B0 C1'] + ['D1'] * (2 * n + 2) + ['D0 E0 B1 E1 A1 A0'])}
        for name, tr in fam.items():
            if 'C0' in tr.split():
                users.add(name)
            rec('L0', 'every F step reads 1', 'F0' not in tr.split(), name)
    rec('L0', 'C0 only in EAT0 and CF', users == {'EAT0', 'CF'}, str(sorted(users)))


def sec_CQ(rng, nsmall, nrand):
    """C(2n) = 00 C(n), C(2n+1) = 11 C(n); C(n) = (11)^k 00 C(n') and C(n+1) = (00)^k 11 C(n') with k = t1(n);
    Q identities; the decompositions used by Y1-Y4 / OVF."""
    ns = list(range(0, nsmall)) + [rng.randrange(1, 2 ** rng.randint(2, 120)) for _ in range(nrand)]
    for n in ns:
        rec('CQ', 'C(2n)=00C(n)', norm_R('00' + Cw(n)) == Cw(2 * n), n)
        rec('CQ', 'C(2n+1)=11C(n)', norm_R('11' + Cw(n)) == Cw(2 * n + 1), n)
        k = t1(n)
        n2 = n >> (k + 1)
        rec('CQ', 'C(n)=(11)^k00C(n2)', norm_R('11' * k + '00' + Cw(n2)) == Cw(n), n)
        rec('CQ', 'C(n+1)=(00)^k11C(n2)', norm_R('00' * k + '11' + Cw(n2)) == Cw(n + 1), n)
        rec('CQ', 'incR(C(n))=C(n+1)', incR(Cw(n)) == (Cw(n + 1), k), n)
        rec('CQ', 'is_C(C(n))=n', is_C(Cw(n)) == n, n)
        if n % 2 == 1:   # Y1 on C(N), N odd: (00)^(s-1) 11 C(n2) = C((N+1)/2)
            rec('CQ', 'Y1:C(N)->C((N+1)/2)', norm_R('00' * (k - 1) + '11' + Cw(n2)) == Cw((n + 1) // 2), n)
        if n % 4 == 0 and n > 0:
            rec('CQ', 'Y2:C(N)=0000C(N/4)', norm_R('0000' + Cw(n // 4)) == Cw(n), n)
        if n % 4 == 2:
            s, m = split6(n)
            rec('CQ', 'Y3:C(N)=0011(11)^s00C(m)', norm_R('0011' + '11' * s + '00' + Cw(m)) == Cw(n), n)
            rec('CQ', 'Y3:result=Q_s(0,2m+1)', norm_R('0' * (2 * s + 1) + '11' + Cw(m)) == norm_R(Qw(s, 0, 2 * m + 1)), n)
    for k in range(0, 9):
        for x in range(2 ** k):
            for y in [1, 3, 5, 7, 9, 11, 2 * rng.randrange(1, 10 ** 6) + 1]:
                Q = Qw(k, x, y)
                if x + 1 < 2 ** k:
                    rec('CQ', 'incR(Q_k(x,y))=Q_k(x+1,y)', incR(Q) == (norm_R(Qw(k, x + 1, y)), t1(x)), (k, x, y))
                else:
                    m = (y - 1) // 2
                    rec('CQ', 'Q_k(2^k-1,2m+1)=(11)^k 011 C(m)', norm_R('11' * k + '011' + Cw(m)) == norm_R(Q), (k, y))
                if k >= 2 and x == 1:
                    rec('CQ', 'Y1:Q_k(1,y)->Q_(k-1)(1,y)', norm_R('11' + Q[4:]) == norm_R(Qw(k - 1, 1, y)), (k, y))
                if k == 1 and x == 1:
                    rec('CQ', 'Y4:Q_1(1,2m+1)=11011C(m)', norm_R('11011' + Cw((y - 1) // 2)) == norm_R(Q), y)
                # 2.3: for odd y, Q_k(x, y) is not a counter word (pair k is 0 1); for y = 0 it is C(x)
                rec('CQ', 'Q_k(x,y) odd y: not a counter word, pair k = 01', is_C(Q) is None and Q[2 * k:2 * k + 2] == '01', (k, x, y))
            rec('CQ', 'Q_k(x,0) = C(x) (the y = 0 exception)', is_C(Qw(k, x, 0)) == x, (k, x))


def L1_specs(b, k, j):
    """(name, lhs, rhs, steps) of the L1 rules (proof.md 3.2, 3.3) with parameters b, k, j"""
    S = []
    S.append(('INC1', f"1 {P('10', b)} 0 1 1 B> 0 {P('11', k)} 0 0", f"0 {P('10', b)} 0 1 1 B> 0 {P('00', k)} 1 1", 6 * b + 4 * k + 24))
    S.append(('INC0', f"0 1 0 0 {P('10', b)} 0 1 1 B> 0 {P('11', k)} 0 0", f"1 0 1 0 {P('10', b)} 0 1 1 B> 0 {P('00', k)} 1 1", 6 * b + 4 * k + 32))
    S.append(('INC011', f"1 0 1 1 0 0 {P('10', b)} 0 1 1 B> 0 {P('11', k)} 0 0", f"0 1 0 1 1 0 {P('10', b)} 0 1 1 B> 0 {P('00', k)} 1 1", 6 * b + 4 * k + 40))
    S.append(('OVF', f"1 {P('10', b)} 0 1 1 B> 0 {P('11', k)} 0 1 1 {P('11', j)} 0 0",
              f"0 {P('10', b + 1)} {P('1', 2 * k + 1)} 0 1 1 B> 0 {P('00', j)} 1 1", 6 * b + 2 * k + 4 * j + 29))
    S.append(('EAT', f"0 0 0 {P('10', b)} 0 1 1 B> 0", f"{P('10', b + 2)} 1 1 0 A>", 4 * b + 17))
    if k >= 1:
        S.append(('Y1', f"1 1 0 A> {P('11', k)} 0 0", f"1 0 1 1 B> 0 {P('00', k - 1)} 1 1", 4 * k + 9))
    S.append(('Y2', "1 1 0 A> 0 0 0 0", "1 1 1 0 1 1 B> 0", 7))
    S.append(('Y3', f"1 1 0 A> 0 0 1 1 {P('11', k)} 0 0", f"1 1 0 1 1 B> 0 {P('0', 2 * k + 1)} 1 1", 4 * k + 20))
    S.append(('Y4', f"1 1 0 A> 1 1 0 1 1 {P('11', j)} 0 0", f"1 1 0 1 0 1 1 B> 0 {P('00', j)} 1 1", 4 * j + 18))
    return S


def sec_L1(rng, small, nrand, big):
    for b in range(small + 1):
        for k in range(small + 1):
            for j in (0, 1, 2, k):
                for name, l, r, st in L1_specs(b, k, j):
                    if name in ('Y1', 'Y2', 'Y3', 'Y4', 'EAT') and (b > 0 and name != 'EAT'):
                        continue
                    check_rule('L1', name, l, r, st, info=f'b={b},k={k},j={j}')
    for _ in range(nrand):
        b, k, j = rng.randint(0, big), rng.randint(0, big // 4), rng.randint(0, big // 4)
        for name, l, r, st in L1_specs(b, k, j):
            check_rule('L1', name + '(random)', l, r, st, info=f'b={b},k={k},j={j}')
    # the L1 rules as compositions of the L0 rules: transition strings concatenate
    def cat(*parts):
        return ' '.join(x for x in parts if x)
    for b in range(0, 6):
        for k in range(0, 5):
            TSs = 'B0 C1 D1 D0 E0 B1'
            EBs = ' '.join(['E0 B1'] * b)
            def AS(n): return ' '.join(['A1 A0 B1 E1'] * n)
            def A1(n): return ' '.join(['A1'] * n)
            def CAR(kk): return ' '.join(['A0 B0 C1'] + ['D1'] * (2 * kk + 2) + ['D0 E0 B1 E1 A1 A0'])
            right = cat(A1(2 * k + 2), CAR(k))
            spec = {s[0]: s for s in L1_specs(b, k, 1)}
            for name, tr in [('INC1', cat(TSs, EBs, 'E1', AS(b + 1), right)),
                             ('INC0', cat(TSs, EBs, 'E0 B0 C1 D0 E1', AS(b + 2), right)),
                             ('INC011', cat(TSs, EBs, 'E0 B0 C1 D1 D0 E0 B1 E1 A1 A0 B1 E1 A1 A1', 'A0 B1 E1', AS(b + 1), right)),
                             ('OVF', cat(TSs, EBs, 'E1', AS(b + 1), A1(2 * k + 2), 'A0 B1 E1', A1(2), A1(2), CAR(1))),
                             ('EAT', cat(TSs, EBs, 'E0 B0 C0 F1', ' '.join(['C0 F1'] * (b + 2)), 'C1 D0 E1'))]:
                n_, l, r, st = spec[name]
                check_rule('L1', name + '=composition', l, r, st, tr, f'b={b},k={k}')
            if b == 0:
                for name, tr in [('Y1', cat(A1(2 * k), CAR(k - 1))) if k >= 1 else (None, None),
                                 ('Y2', 'A0 B0 C1 D0 E1 A1 A0'),
                                 ('Y3', cat('A0 B0 C1 D0 E1 A1', A1(2 * k + 2), 'A0 B0 C1', ' '.join(['D1'] * (2 * k + 3)), 'D0 E0 B1', 'E1 A1 A0')),
                                 ('Y4', cat(A1(2), 'A0 B1 E1', A1(2), A1(2), CAR(1)))]:
                    if name is None:
                        continue
                    n_, l, r, st = spec[name]
                    check_rule('L1', name + '=composition', l, r, st, tr, f'k={k}')


def sec_X(rng, nstates, pmax):
    """rule-level executor = direct TM, at every G-configuration (configuration and step count)"""
    cnt = 0
    for _ in range(nstates):
        p = rng.randint(0, pmax)
        N = rng.randint(1, 2 * p + 8)
        w, R = '10' * p, Cw(N)
        tape, pos, st2 = gtape(w, R, 8 * p + 400)
        ok = True
        for g in range(4 * p + 40):
            r = gstep(w, R)
            if r is None:
                break
            name, w, R, s = r
            try:
                pos, st2 = fast_run(tape, pos, st2, s)
            except Halted:
                ok = False
                break
            if tm_config(tape, pos, st2) != gconfig(w, R):
                ok = False
                break
            cnt += 1
        rec('X', 'executor=TM_at_every_G', ok, f'p={p},N={N},g={g}')
    rec('X', 'G_transitions_compared', cnt > 0, str(cnt))
    print(f"  X: {cnt} G transitions compared with the direct TM")


def sec_W(rng, amax, nrand):
    """W(a,b): G((10)^a 1 (10)^b, R) -> G((10)^(a+b+[a odd]), R^(+a+1)) in a+1 counts INC1,INC0,...;
    W2: G(u 1 0 1 0 1 (10)^b) -> G(u 1 (10)^(b+2)) by INC1, INC0 (any u);
    W0: G(u 0 (10)^(2c+1) 0 (10)^d) -> G(u (10)^(2c+2+d)) in 2c+1 counts INC0, INC1, ..., INC0 (any u);
    block stage: G(u 1^(2t+1)) -> G(u 0 (10)^t) in t+1 counts INC1 (any u)."""
    def run(w, R, n):
        names = []
        for _ in range(n):
            r = gstep(w, R)
            if r is None:
                return None, None, names
            names.append(r[0])
            w, R = r[1], r[2]
        return w, R, names
    for a in range(amax + 1):
        for b in range(0, 6):
            n0 = rng.randint(1, 10 ** 6)
            w, R, names = run('10' * a + '1' + '10' * b, Cw(n0), a + 1)
            exp = ['INC1', 'INC0'] * ((a + 1) // 2) + (['INC1'] if a % 2 == 0 else [])
            rec('W', 'W(a,b)', (w, R, names) == ('10' * (a + b + a % 2), Cw(n0 + a + 1), exp), f'a={a},b={b}')
    for _ in range(nrand):
        u = ''.join(rng.choice('01') for _ in range(rng.randint(0, 12)))
        b = rng.randint(0, 30)
        n0 = rng.randint(1, 10 ** 9)
        w, R, names = run(u + '10101' + '10' * b, Cw(n0), 2)
        rec('W', 'W2(any u)', (w, R, names) == (norm_w(u + '1' + '10' * (b + 2)), Cw(n0 + 2), ['INC1', 'INC0']), (u, b))
        c, d = rng.randint(0, 12), rng.randint(0, 20)
        w, R, names = run(u + '0' + '10' * (2 * c + 1) + '0' + '10' * d, Cw(n0), 2 * c + 1)
        rec('W', 'W0(any u)', (w, R, names) == (norm_w(u + '10' * (2 * c + 2 + d)), Cw(n0 + 2 * c + 1), ['INC0', 'INC1'] * c + ['INC0']), (u, c, d))
        t = rng.randint(0, 12)
        w, R, names = run(u + '1' * (2 * t + 1), Cw(n0), t + 1)
        rec('W', 'block(any u)', (w, R, names) == (norm_w(u + '0' + '10' * t), Cw(n0 + t + 1), ['INC1'] * (t + 1)), (u, t))
        # W2 on Q_t: x -> x+2 while x+2 <= 2^t - 1
        t = rng.randint(2, 10)
        x = 2 * rng.randrange(0, 2 ** (t - 1) - 1) + 1
        y = 2 * rng.randrange(0, 10 ** 6) + 1
        w, R, names = run(u + '10101' + '10' * b, Qw(t, x, y), 2)
        rec('W', 'W2 on Q_t', (w, R) == (norm_w(u + '1' + '10' * (b + 2)), norm_R(Qw(t, x + 2, y))), (u, t, x))
        # OVF on Q_t(2^t-1, 2m+1)
        m = rng.randrange(0, 10 ** 6)
        w, R, names = run(u + '1' + '10' * b, Qw(t, 2 ** t - 1, 2 * m + 1), 1)
        rec('W', 'OVF on Q_t', (w, R, names) == (norm_w(u + '0' + '10' * (b + 1) + '1' * (2 * t + 1)), Cw(m + 1), ['OVF']), (u, t, b, m))


def sec_M(rng, pmax, nrand, pbig, ntm, ptm, watch_cap):
    # (1) chain of proof.md = rule-level executor (rule names, steps, end state) = closed form, for every small state
    nst = 0
    for p in range(0, pmax + 1):
        for N in range(1, 2 * p + 1):
            pp, NN, L, case = chain_EE(p, N)
            q = exec_EE(p, N)
            names = [(c[0], c[3]) for c in q[2]]
            rec('M', 'chain=executor(names,steps)', pp is not None and names == L and (q[0], q[1]) == (pp, NN), f'p={p},N={N},{case}')
            rec('M', 'closed_form=chain', (pp, NN) == model(p, N), f'p={p},N={N}')
            nst += 1
    rec('M', 'states_checked', True, str(nst))
    # (2) random large states (executor on strings; N up to 2^60, p up to pbig)
    for _ in range(nrand):
        p = rng.randint(1, pbig)
        kind = rng.randrange(5)
        if kind == 4:   # 6 mod 8 with large t near the precondition
            t = rng.randint(1, max(1, (p + 7).bit_length() - 1))
            m = rng.randrange(0, 4)
            N = (2 * m + 1) * 2 ** (t + 3) - 2
        else:
            N = rng.randrange(1, 2 ** rng.randint(2, 60))
            N = N - (N % 8) + [1, 4, 2, 6][kind] if kind < 4 else N
            N = max(N, 1)
        pp, NN, L, case = chain_EE(p, N)
        if pp is None:
            continue
        q = exec_EE(p, N)
        names = [(c[0], c[3]) for c in q[2]]
        rec('M', 'chain=executor(random large)', names == L and (q[0], q[1]) == (pp, NN), f'p={p},N={N},{case}')
        rec('M', 'closed_form=chain(random large)', (pp, NN) == model(p, N), f'p={p},N={N}')
    # (3) direct TM: E(p,N) runs exactly sum(steps) to E(p',N'), and passes every G of the chain at the stated steps
    cases = []
    for _ in range(ntm):
        p = rng.randint(0, ptm)
        N = rng.randint(1, 2 * p + 2)
        cases.append((p, N))
    for t in range(1, 7):        # 6 mod 8 with each t, at and near the precondition boundary
        for m in (0, 1, 2):
            N = (2 * m + 1) * 2 ** (t + 3) - 2
            for p in sorted({max(0, 2 ** t - 7), max(0, 2 ** t - 6), 2 ** t, (N + 1) // 2}):
                cases.append((p, N))
    nwatch = [0, 0]
    for p, N in cases:
        pp, NN, L, case = chain_EE(p, N)
        if pp is None:
            continue
        q = exec_EE(p, N)
        w, R = '10' * p, Cw(N)
        tape, pos, st2 = gtape(w, R, 8 * p + 64 * (case.count('t=') + 1) + 400)
        ok = True
        try:
            for name, w2, R2, s in q[2]:
                pos, st2 = fast_run(tape, pos, st2, s)
                if tm_config(tape, pos, st2) != gconfig(w2, R2):
                    ok = False
                    break
        except Halted:
            ok = False
        ok = ok and tm_config(tape, pos, st2) == gconfig('10' * pp, Cw(NN))
        rec('M', f'direct_TM[{case.split(",")[0]}]', ok, f'p={p},N={N},{case}')
        # every TM step (not only rule boundaries): the first eat configuration after the start is E(p', N'),
        # reached exactly at the end of the chain
        S = sum(s for _, s in L)
        if S <= watch_cap:
            tape, pos, st2 = gtape('10' * p, Cw(N), 8 * p + 64 * (case.count('t=') + 1) + 400)
            try:
                pos, st2, found = watch_E(tape, pos, st2, S)
                okw = found == [(S, (pp, NN))]
            except Halted:
                okw, found = False, 'halt'
            rec('M', 'every step: first eat configuration = E(p\',N\') at the end', okw, f'p={p},N={N},{case}: {found[:3]}')
            nwatch[0] += 1
            nwatch[1] += S
    rec('M', 'per-step watched states', nwatch[0] > 0, str(nwatch))
    print(f"  M: per-step E-form watch on {nwatch[0]} direct runs, {nwatch[1]:,} TM steps")


def sec_INV(rng, pmax, nrand):
    """J(p,N): 1 <= N <= 2p.  Coverage: J -> the chain of 3.5 applies (6 mod 8: p+7 >= 2^t).  Closure: J(p',N')."""
    n = 0
    for p in range(1, pmax + 1):
        for N in range(1, 2 * p + 1):
            pp, NN = model(p, N)
            okc = True
            if N % 8 == 6:
                s, m = split6(N)
                if s >= 2:
                    t = s - 1
                    okc = p + 7 >= 2 ** t and 2 ** (t + 2) * (2 * m + 1) <= p + 1
            rec('INV', 'coverage(precondition)', okc, f'p={p},N={N}')
            rec('INV', 'closure', 1 <= NN <= 2 * pp, f'p={p},N={N}->({pp},{NN})')
            rec('INV', 'progress p\'>=p+2', pp >= p + 2, f'p={p},N={N}')
            n += 1
    rec('INV', 'exhaustive_states', True, str(n))
    # random large states, using the same bounds as proof.md 4
    for _ in range(nrand):
        p = rng.randrange(1, 2 ** rng.randint(1, 80))
        N = rng.randrange(1, 2 * p + 1)
        if rng.random() < 0.3:   # force 6 mod 8 with a large t
            t = rng.randint(0, max(0, (p + 1).bit_length() - 3))
            mmax = ((p + 1) // 2 ** (t + 2) - 1) // 2
            if mmax >= 0:
                m = rng.randint(0, mmax)
                N = (2 * m + 1) * 2 ** (t + 3) - 2
        if not (1 <= N <= 2 * p):
            continue
        pp, NN = model(p, N)
        rec('INV', 'closure(random large)', 1 <= NN <= 2 * pp, f'p={p},N={N}')
        if N % 8 == 6:
            s, m = split6(N)
            if s >= 2:
                t = s - 1
                rec('INV', '6mod8 bounds(random large)', 2 ** (t + 2) * (2 * m + 1) <= p + 1 and 2 ** t <= p + 7
                    and NN <= p + t + 2, f'p={p},N={N}')
        if N % 2 == 1:
            rec('INV', 'odd: (N+1)/2 <= p', (N + 1) // 2 <= p, f'p={p},N={N}')


def sec_Z(full, rng):
    # initial segment: blank -> E(0,1) at step 31 -> E(2,4) at step 165
    tape = bytearray(4096)
    pos, st2 = 2048, 0
    pos, st2 = fast_run(tape, pos, st2, 31)
    rec('Z', 'blank->E(0,1) at step 31', tm_config(tape, pos, st2) == gconfig('', Cw(1)), str(tm_config(tape, pos, st2)))
    s01 = macro_steps(0, 1)
    rec('Z', 'E(0,1)->E(2,4) in 134 steps (chain)', model(0, 1) == (2, 4) and s01 == 134, str(s01))
    pos, st2 = fast_run(tape, pos, st2, s01)
    rec('Z', 'blank->E(2,4) at step 165', tm_config(tape, pos, st2) == gconfig('1010', Cw(4)), str(tm_config(tape, pos, st2)))
    rec('Z', 'J(2,4)', 1 <= 4 <= 4)
    # H form (proof.md 8): S(p, N) := H((10)^p, C(N)) is G after TS, 6 steps later
    pos, st2 = fast_run(tape, pos, st2, 6)
    rec('Z', 'blank->H((10)^2,C(4)) at step 171 (G at 165, then TS)', tm_config(tape, pos, st2) == hconfig('1010', Cw(4)), str(tm_config(tape, pos, st2)))
    s24 = macro_steps(2, 4)
    pos, st2 = fast_run(tape, pos, st2, s24)
    rec('Z', 'H form to H form takes the macro steps (E(2,4)->E(5,7))', model(2, 4) == (5, 7) and tm_config(tape, pos, st2) == hconfig('10' * 5, Cw(7)), f'{s24}: {tm_config(tape, pos, st2)}')
    # the G-configurations before E(0,1): exactly one, G(1, C(0)) at step 7, and INC1 (b = 0, k = 0, 24 steps) takes it to E(0,1)
    tape = bytearray(4096)
    pos, st2 = 2048, 0
    gs = []
    for i in range(32):
        if st2 == 2 and tape[pos] == 0 and tape[pos - 1] == 1 and tape[pos - 2] == 1 and tape[pos - 3] == 0:
            gs.append((i, tm_config(tape, pos, st2)))
        if i < 31:
            pos, st2 = fast_run(tape, pos, st2, 1)
    rec('Z', 'G-configs to step 31: G(1,C(0)) at 7, E(0,1) at 31', gs == [(7, gconfig('1', Cw(0))), (31, gconfig('', Cw(1)))], str(gs))
    rec('Z', 'G(1,C(0)) -INC1-> E(0,1) by the executor in 24 steps', gstep('1', Cw(0)) == ('INC1', '', Cw(1), 24), str(gstep('1', Cw(0))))
    tape = bytearray(4096)
    pos, st2 = 2048, 0
    tr = []
    for i in range(7):
        tr.append(ST[st2 // 2] + str(tape[pos]))
        pos, st2 = fast_run(tape, pos, st2, 1)
    rec('Z', 'blank->G(1,C(0)): A0 B0 C1 D0 E1 A1 A0', ' '.join(tr) == 'A0 B0 C1 D0 E1 A1 A0' and tm_config(tape, pos, st2) == gconfig('1', ''), ' '.join(tr))
    # stored direct run (C, 1e11 steps): every eat E(p,N) with its step, versus the model and the step formula
    path = os.path.join(HERE, 'work', 'mev_1e11.txt')
    eats = []
    if os.path.exists(path):
        for line in open(path):
            f = line.split()
            if f[0] == 'E' and int(f[3]) % 2 == 0:
                step, n, d, v = int(f[1]), int(f[3]), int(f[4]), int(f[5])
                if d >= 62:
                    continue
                eats.append((step, n // 2, (1 << d) + v))
        p, N, step = 0, 1, 31
        ok = eats[0] == (31, 0, 1)
        bad = 0
        for (st_, pp_, NN_) in eats[1:]:
            step += macro_steps(p, N)
            p, N = model(p, N)
            if (step, p, N) != (st_, pp_, NN_):
                bad += 1
                if bad < 3:
                    FAILS.append(f"Z mismatch at stored eat {(st_, pp_, NN_)} vs model {(step, p, N)}")
        rec('Z', f'stored run: {len(eats)} eat configurations ({len(eats) - 1} transitions) = model incl. step numbers', ok and bad == 0 and len(eats) > 2000, f'bad={bad}')
    else:
        rec('Z', 'stored run file present', False, path)
    # fresh direct run (this script's simulator) from blank, compared at every eat
    budget = 4 * 10 ** 8 if full else 6 * 10 ** 7
    tape = bytearray(1 << 16)
    pos, st2 = 1 << 15, 0
    pos, st2 = fast_run(tape, pos, st2, 31)
    p, N, step, ne = 0, 1, 31, 0
    ok = True
    while True:
        s = macro_steps(p, N)
        if step + s > budget:
            break
        pos, st2 = fast_run(tape, pos, st2, s)
        step += s
        p, N = model(p, N)
        ne += 1
        if tm_config(tape, pos, st2) != gconfig('10' * p, Cw(N)):
            ok = False
            break
        if not (1 <= N <= 2 * p):
            ok = False
            break
    rec('Z', 'fresh direct run: every eat configuration = model, J along the orbit', ok, f'{ne + 1} eat configurations to step {step}, last E({p},{N})')
    print(f"  Z: fresh direct run: {ne + 1} eat configurations ({ne} transitions), last E({p},{N}) at step {step:,}")
    # every TM step from blank: the eat configurations met are exactly those of the model orbit
    wbudget = 2 * 10 ** 7 if full else 3 * 10 ** 6
    want = [(31, (0, 1))]
    p, N, step = 0, 1, 31
    while True:
        step += macro_steps(p, N)
        p, N = model(p, N)
        if step > wbudget:
            break
        want.append((step, (p, N)))
    tape = bytearray(1 << 14)
    pos, st2, found = watch_E(tape, 1 << 13, 0, wbudget)
    rec('Z', 'every step from blank: eat configurations = model orbit', found == want,
        f'{len(found)} found, {len(want)} expected to step {wbudget}; first diff ' +
        str(next(((a, b) for a, b in zip(found, want) if a != b), None)))
    print(f"  Z: every step to {wbudget:,}: {len(found)} eat configurations, all on the model orbit: {found == want}")
    # abstract orbit: J along the model orbit
    p, N = 2, 4
    okj = strict = True
    best = (0, 1, None)          # max N/p as (N, p, index)
    nmax = 2_000_000 if full else 300_000
    for i in range(nmax):
        p, N = model(p, N)
        if not (1 <= N <= 2 * p):
            okj = False
            break
        if N >= 2 * p:
            strict = False
        if N * best[1] > best[0] * p:
            best = (N, p, i + 1)
    rec('Z', f'J at the {nmax} model eat configurations after E(2,4)', okj, f'i={i},p={p},N={N}')
    rec('Z', 'N < 2p strictly after E(2,4) (proof.md 4.2 remark)', strict, f'max N/p = {best}')
    print(f"  Z: max N/p over {nmax:,} eat configurations after E(2,4): {best[0] / best[1]:.7f} at #{best[2]} E({best[1]},{best[0]})")


def sec_NEG(rng):
    """below the precondition (p+7 < 2^t) the chain of 3.5 does not describe the run"""
    for t in range(4, 8):
        for m in (0, 1):
            N = (2 * m + 1) * 2 ** (t + 3) - 2
            p = 2 ** t - 8 - rng.randint(0, 3)
            q = exec_EE(p, N)
            mo = (p + 7 + t, m + 2 ** t + t + 1)
            rec('NEG', 'p+7<2^t: run differs from the 6 mod 8 formula', q[0] is None or (q[0], q[1]) != mo, f'p={p},N={N}')


def sec_TEETH():
    """the windowed checks have teeth: each single-transition mutation that the rules use is caught"""
    global T
    T0 = dict(T)
    saved = (dict(RES), list(FAILS))
    used = ['B0', 'C1', 'D1', 'D0', 'E0', 'B1', 'E1', 'A1', 'A0', 'C0', 'F1']
    caught = {}
    for key in used:
        q, sym = key[0], int(key[1])
        w, d, nx = T0[(q, sym)]
        muts = [(1 - w, d, nx), (w, -d, nx)] + [(w, d, z) for z in ST if z != nx]
        for mt in muts:
            T = dict(T0)
            T[(q, sym)] = mt
            RES.clear()
            for b in range(0, 3):
                for k in range(0, 3):
                    for name, l, r, st in L1_specs(b, k, 1):
                        check_rule('tmp', name, l, r, st)
            bad = sum(v[1] for v in RES.values())
            caught[(key, mt)] = bad > 0
    T = T0
    RES.clear()
    RES.update(saved[0])
    FAILS[:] = saved[1]
    nc = sum(caught.values())
    rec('TEETH', 'every single-transition mutation fails some L1 check', nc == len(caught), f'{nc}/{len(caught)}: ' + str([k for k, v in caught.items() if not v][:5]))
    print(f"  TEETH: {nc}/{len(caught)} mutations caught")


# ================================================================ main
def main():
    ap = argparse.ArgumentParser()
    ap.add_argument('--seed', type=int, default=1)
    ap.add_argument('--quick', action='store_true')
    ap.add_argument('--full', action='store_true')
    a = ap.parse_args()
    rng = random.Random(a.seed)
    Q, F = a.quick, a.full
    t0 = time.time()
    sec_T()
    sec_L0(40 if not Q else 12)
    sec_CQ(rng, 4096 if not Q else 256, 3000 if not Q else 200)
    sec_L1(rng, 10 if not Q else 4, 600 if not Q else 50, 300 if not Q else 60)
    sec_X(rng, 200 if not Q else 20, 120 if not Q else 30)
    sec_W(rng, 40 if not Q else 10, 1500 if not Q else 100)
    sec_M(rng, (140 if F else 80) if not Q else 25, 400 if not Q else 30, 2000 if not Q else 200,
          (400 if F else 150) if not Q else 15, 400 if not Q else 60, (10 ** 6 if F else 3 * 10 ** 5) if not Q else 5 * 10 ** 4)
    sec_INV(rng, (3000 if F else 1500) if not Q else 200, 300000 if not Q else 10000)
    sec_Z(F, rng)
    sec_NEG(rng)
    sec_TEETH()
    npass = sum(1 for v in RES.values() if v[1] == 0)
    for (sec, name), (ok, bad) in sorted(RES.items(), key=lambda kv: kv[0]):
        print(f"  {'PASS' if bad == 0 else 'FAIL'} {sec}.{name}: {ok} ok, {bad} bad")
    for f in FAILS:
        print(f)
    print(f"{npass}/{len(RES)} PASS  (seed {a.seed}, {'quick' if Q else 'full' if F else 'default'}, {time.time() - t0:.1f} s)")
    sys.exit(0 if npass == len(RES) else 1)


if __name__ == '__main__':
    main()