#!/usr/bin/env python3
"""check_proof.py: evidence for the DELTA proof of the twin 1RB---_0LC1RF_1LE0RD_0RB1RC_1RD0LE_0RC1RA.

Parent: kill #4, 1RB0LA_0RC1RD_0LD1RE_1LA0RB_0RD1RF_1RC--- (KILLED_NONHALT). The parent's check_proof.py is
loaded read-only for its abstract model only: the rule list of its proof.md 3.1 (words, claimed step counts),
its derivation runner (which applies exactly the chains of its proof.md 3.2-3.3 and sums the claimed step
counts), the families P/Q, the cycle map S/F and the hexagon. Everything that depends on the machine or on the
start is re-run here on the TWIN's own transition table.

What kind of check each section is:
  EXACT (proves a finite claim):  S (the renaming search and the start state), A (the 21 shift rules,
                                  renamed, executed on the twin's table), I (blank tape -> base, 16,590 steps),
                                  B (the invariant at the base), H (the polygon, parent's exact check re-run).
  SAMPLE/RANGE (evidence):        L (lockstep with the parent run from F>), M (macro rules on the twin's
                                  table), C (every canonical visit to 3e7), O (the twin's orbit vs the cycle
                                  map, to k = 14 by default).
Nothing here is a Rocq proof; only tools/check_kill.sh PASS counts.

  S  sigma = EDBCFA (parent A->E, B->D, C->B, D->C, E->F, F->A) is the only state renaming, with or without
     mirroring, that takes the parent's table to the twin's; the twin's start A = sigma(F)
  L  the twin from blank equals the relabelled parent run from 0^inf F> 0^inf, step by step
  A  the parent's 21 shift rules with head states renamed by sigma (<A -> <E, C> -> B>, E> -> F>), executed
     exactly on the twin's table: same step counts, the halting cell A1 = sigma(F1) never read
  M  R1, R2, R3, R4, R5, R5o (and R4o): the parent's derivation runner predicts the result and the step count;
     the twin's table is simulated directly from the renamed start for exactly that many steps
  I  blank -> Q(7,65,1,14,4) with head B = sigma(C) at exactly step 16,590; also the rule-governed prefix
     from Q(5,20,3,1,0) at step 1,711 (R5, R1^4, R2, R3^2, R5, R1^3, R2) with exact steps
  C  every canonical right-edge visit of the twin from step 1,711 to 3e7, vs the rule runner (config and step)
  B  the base S(6,14,4) = CfgF 5 0 14 4: the parent's Rocq Inv, the hexagon N*R (facet 4Y+2Z >= N tight),
     cycle-lemma hypotheses (a)-(d)
  H  the parent's exact rational hexagon check (its section 5), re-run unchanged
  O  the twin's orbit: S(k,Y,Z) -> S(k+1,F(Y,Z)) from the base; each predicted cycle state and its exact
     step (runner) is compared with direct simulation of the twin, chained, to S(kmax,...); all points in R
Usage: python3 -B check_proof.py [--quick] [--kmax K] [--lock N]
"""
import os, sys, json, random, argparse, itertools, importlib.util
sys.dont_write_bytecode = True
sys.path.insert(0, "[local working file]")
import sim  # noqa: E402

TM = '1RB---_0LC1RF_1LE0RD_0RB1RC_1RD0LE_0RC1RA'
PARENT = '1RB0LA_0RC1RD_0LD1RE_1LA0RB_0RD1RF_1RC---'
HERE = os.path.dirname(os.path.abspath(__file__))
PDIR = os.path.join(os.path.dirname(HERE), PARENT)

ap = argparse.ArgumentParser()
ap.add_argument('--quick', action='store_true')
ap.add_argument('--kmax', type=int, default=14, help='check the orbit up to S(kmax, ...) (14: step ~1.27e9)')
ap.add_argument('--lock', type=float, default=2e6, help='lockstep length')
ARGS = ap.parse_args()

_argv = sys.argv
sys.argv = [_argv[0]] + (['--quick'] if ARGS.quick else [])
spec = importlib.util.spec_from_file_location('parent_check_proof', os.path.join(PDIR, 'check_proof.py'))
cp = importlib.util.module_from_spec(spec)
spec.loader.exec_module(cp)                     # parent's checks, read-only (nothing runs at import)
sys.argv = _argv
assert cp.TM == PARENT

ST = 'ABCDEF'
SIGMA = dict(zip(ST, 'EDBCFA'))                 # parent state -> twin state
MT = sim.parse(TM)
FAILS = []


def ok(cond, msg):
    if not cond:
        FAILS.append(msg)
        print('   FAIL:', msg)
    return cond


def letters(s):
    t = {}
    for i, blk in enumerate(s.split('_')):
        for r in (0, 1):
            x = blk[3 * r:3 * r + 3]
            t[(ST[i], r)] = None if x == '---' else (int(x[0]), x[1], x[2])
    return t


TP, TT = letters(PARENT), letters(TM)


def image(tp, sig, mirror):
    out = {}
    for (q, r), v in tp.items():
        out[(sig[q], r)] = None if v is None else (v[0], {'L': 'R', 'R': 'L'}[v[1]] if mirror else v[1], sig[v[2]])
    return out


def tw(cfgstr):
    """sigma on a parent configuration string: rename the head state (tape words carry no states)."""
    out = []
    for t in cfgstr.split():
        if t.startswith('<') and t[1:] in ST:
            t = '<' + SIGMA[t[1:]]
        elif t.endswith('>') and t[:-1] in ST:
            t = SIGMA[t[:-1]] + '>'
        out.append(t)
    return ' '.join(out)


def tw_sim_norm(cfgstr, steps):
    r = sim.run(MT, sim.Config.parse(cfgstr), steps)
    if r.halted or r.step != steps:
        return ('HALT' if r.halted else 'SHORT', r.step)
    c = r.final
    return cp.norm_cells(list(c.tape), c.head, c.state)


# ------------------------------------------------------------------------------------------------ S
def section_S():
    print('[S] relation: renaming / mirror search, start state, Eqv_v2 TM7')
    sols = [(''.join(p), mir) for p in itertools.permutations(ST) for mir in (False, True)
            if image(TP, dict(zip(ST, p)), mir) == TT]
    print('   maps phi with twin = phi(parent) among 720 renamings x {plain, mirror}:',
          [(s, 'mirror' if m else 'plain') for s, m in sols])
    ok(sols == [('EDBCFA', False)], f'renaming not unique / not EDBCFA: {sols}')
    inv = {v: k for k, v in SIGMA.items()}
    print(f'   twin start A = sigma({inv["A"]}): the twin from blank is sigma(parent from 0^inf {inv["A"]}> 0^inf)')
    ok(inv['A'] == 'F', 'twin start is not sigma(F)')
    # halting cells
    ok([k for k, v in TT.items() if v is None] == [('A', 1)] and [k for k, v in TP.items() if v is None] == [('F', 1)],
       'halting cells')
    print('   only undefined cells: twin A1 = sigma(parent F1)')
    # parent's own start: A0 = 1RB, B0 = 0RC -> step 2 is 0^inf 1 0 C> 0^inf; from F: F0 = 1RC -> 0^inf 1 C> 0^inf
    ok(TP[('A', 0)] == (1, 'R', 'B') and TP[('B', 0)] == (0, 'R', 'C') and TP[('F', 0)] == (1, 'R', 'C'), 'starts')
    print('   parent from A: 1RB, 0RC -> 0^inf 1 0 C> 0^inf at step 2;  from F: 1RC -> 0^inf 1 C> 0^inf at step 1'
          ' (different tapes, so different orbits)')
    # Eqv_v2.v TM7 (vendor, read-only): tm' is our string, tm = 1LB0RC_1RC0LB_0RD1RA_0LA1RE_0RA1RF_1RD---
    e7 = '1LB0RC_1RC0LB_0RD1RA_0LA1RE_0RA1RF_1RD---'
    sols7 = [(''.join(p), mir) for p in itertools.permutations(ST) for mir in (False, True)
             if image(TP, dict(zip(ST, p)), mir) == letters(e7)]
    mir7 = '_'.join(b.translate(str.maketrans('LR', 'RL')) for b in e7.split('_'))
    print(f'   Eqv_v2 TM7 tm = {e7} = phi(parent) for {sols7} (started in the parent\'s D);')
    print(f'   its mirror {mir7} is the 1590-file member. Eqv_v2 proves only halts tm c0 <-> halts tm\' c0'
          ' (an equivalence, no nonhalt)')
    ok(sols7 == [('BCDAEF', False)] and mir7 == '1RB0LC_1LC0RB_0LD1LA_0RA1LE_0LA1LF_1LD---', 'Eqv_v2 TM7 relation')


# ------------------------------------------------------------------------------------------------ L
def section_L(n):
    print(f'[L] lockstep: twin from blank vs parent from 0^inf F> 0^inf, relabelled, {n:,} steps')
    sg = {ST.index(q): ST.index(SIGMA[q]) for q in ST}
    tabp = {(ST.index(q), r): (None if v is None else (v[0], 1 if v[1] == 'R' else -1, ST.index(v[2])))
            for (q, r), v in TP.items()}
    tabt = {(ST.index(q), r): (None if v is None else (v[0], 1 if v[1] == 'R' else -1, ST.index(v[2])))
            for (q, r), v in TT.items()}
    size = 2 * n + 16
    a = bytearray(size); b = bytearray(size); pa = pb = n + 8
    sa, sb = 0, ST.index('F')
    bad = None
    for step in range(n):
        if sa != sg[sb] or pa != pb or a[pa] != b[pb]:
            bad = step; break
        ta, tb = tabt[(sa, a[pa])], tabp[(sb, b[pb])]
        if (ta is None) != (tb is None):
            bad = step; break
        if ta is None:
            bad = ('halt', step); break
        a[pa] = ta[0]; pa += ta[1]; sa = ta[2]
        b[pb] = tb[0]; pb += tb[1]; sb = tb[2]
    ok(bad is None and a == b, f'lockstep differs at {bad}')
    print('   identical at every step (same tape and position, twin state = sigma(parent state)), no halt')


# ------------------------------------------------------------------------------------------------ A
def window_exec_twin(r, maxsteps=10 ** 5):
    tape = dict(enumerate(r['wa']))
    n = len(r['wa'])
    pos, st, used = r['oa'], r['sa'], set()
    tgt = dict(enumerate(r['wb']))
    for steps in range(maxsteps):
        if st == r['sb'] and pos == r['ob'] and all(tape.get(i, 0) == v for i, v in tgt.items()) \
                and all(v == 0 for i, v in tape.items() if not 0 <= i < n):
            return steps, used
        inside = 0 <= pos < n or (pos < 0 and r['linf']) or (pos >= n and r['rinf'])
        if not inside:
            return None, used
        sym = tape.get(pos, 0)
        tr = TT[(st, sym)]
        used.add(st + str(sym))
        if tr is None:
            return 'HALT', used
        tape[pos] = tr[0]; pos += 1 if tr[1] == 'R' else -1; st = tr[2]
    return None, used


def section_A():
    print('[A] the 21 shift rules, head states renamed by sigma, executed exactly on the twin\'s table')
    allused = set()
    for name, lhs, rhs, steps in cp.RULES:
        l2, r2 = tw(lhs), tw(rhs)
        got, used = window_exec_twin(cp.norm_rule(l2, r2))
        allused |= used
        ok(got == steps, f'{name}: {l2} -> {r2}: got {got}, parent claims {steps}')
        ok('A1' not in used, f'{name} reads A1')
        print(f'   {name:5s} {l2:22s} -> {r2:24s} {got!s:>4} steps')
    print('   transitions used:', sorted(allused), '(A1 never; all 11 others used)')
    ok(len(allused) == 11, 'not all defined transitions used')
    # the repeated rules, a few sizes, and CTR / OV on the twin's table
    K = 7 if ARGS.quick else 9
    cnt = 0
    for a in range(0, 13):
        for blk_in, blk_out, nm, s in (('00111', '00101', 'AL', 25), ('10111', '00001', 'BL', 9)):
            ok(window_exec_twin(cp.norm_rule((blk_in + ' ') * a + '<E 001', ('<E 001 ' + (blk_out + ' ') * a).strip()),
                                10 ** 6)[0] == s * a, f'{nm}^{a}')
        for blk_in, blk_out, nm, s in (('00101', '00111', 'AR', 7), ('00001', '10111', 'BR', 23)):
            ok(window_exec_twin(cp.norm_rule('001 F> ' + (blk_in + ' ') * a, (blk_out + ' ') * a + '001 F>'), 10 ** 6)[0] == s * a,
               f'{nm}^{a}')
    for k in range(1, K + 1):
        for n in range(2 ** (k - 1), 2 ** k - 1):
            s = window_exec_twin(cp.norm_rule(f'0^inf {cp.C(k, n)} <E 001', f'0^inf {cp.C(k, n + 1)} 001 F>'), 10 ** 6)[0]
            ok(s == 16 * cp.tz1(n) + 18, f'CTR k={k} n={n}: {s}')
            cnt += 1
        s = window_exec_twin(cp.norm_rule(f'0^inf {cp.C(k, 2 ** k - 1)} <E 001 00101 00101',
                                          f'0^inf {cp.C(k + 1, 2 ** k + 1)} 00111 111 001 F>'), 10 ** 6)[0]
        ok(s == 16 * k + 46, f'OV k={k}: {s}')
    print(f'   AL^a, BL^a, AR^a, BR^a for a = 0..12 (25a, 9a, 7a, 23a steps); CTR for all {cnt} (k,n) with k <= {K}'
          f' (16 t(n)+18); OV for k = 1..{K} (16k+46)')


# ------------------------------------------------------------------------------------------------ M
def check_macro(name, start, fn, args, expect):
    T = cp.Tape(cp.build(start))
    try:
        got = fn(T, *args)
    except (ValueError, AssertionError) as e:
        ok(False, f'{name} {start}: parent derivation failed: {e}')
        return
    ok(got == expect, f'{name} {start}: derived {got}, expected {expect}')
    ok(tw_sim_norm(tw(cp.build(start)), T.steps) == cp.cfg_norm(tw(cp.build(got))),
       f'{name} {start}: twin simulation disagrees after {T.steps} steps')


def section_M():
    print('[M] macro rules on the twin\'s table: parent runner (prediction + steps) vs direct twin simulation')
    rnd = random.Random(4417)
    N = 12 if ARGS.quick else 40
    cnt = {}

    def pick(lo, hi):
        return rnd.choice([rnd.randint(lo, min(hi, lo + 3)), rnd.randint(lo, hi)])
    for _ in range(N):
        k = rnd.randint(2, 10); n = rnd.randint(2 ** (k - 1), 2 ** k - 2); X = pick(0, 200); z = pick(0, 200)
        check_macro('R1', ('P', k, n, X, z), cp.R1, (k, n, X, z), ('P', k, n + 1, X, z + 1))
        k = rnd.randint(1, 10); X = pick(2, 200); z = pick(0, 200)
        check_macro('R2', ('P', k, 2 ** k - 1, X, z), cp.R2, (k, 2 ** k - 1, X, z), ('Q', k + 1, 2 ** k + 1, 1, X - 2, z + 1))
        k = rnd.randint(3, 10); n = rnd.randint(2 ** (k - 1), 2 ** k - 4); x = pick(1, 200); y = pick(2, 200); z = pick(0, 200)
        check_macro('R3', ('Q', k, n, x, y, z), cp.R3, (k, n, x, y, z), ('Q', k, n + 3, x + 2, y - 2, z + 1))
    cnt.update({'R1': N, 'R2': N, 'R3': N})
    done = {'R4': 0, 'R5': 0, 'R4o': 0, 'R5o': 0}
    tries = 0
    while min(done.values()) < N and tries < 500 * N:
        tries += 1
        y = rnd.randint(0, 1); k = rnd.randint(4, 10); x = pick(1, 2 ** (k - 1)); z = pick(0, 2 ** (k - 1))
        T_ = z + 2 if y == 0 else 2 * z + 7
        if rnd.random() < 0.5:
            lo, hi = max(2 ** (k - 1), 2 ** k - T_), 2 ** k - 1
        else:
            lo, hi = 2 ** (k - 1), 2 ** k - 1 - T_
        if lo > hi:
            continue
        n = rnd.randint(lo, hi)
        if not cp.in_domain_conv(k, n, x, y, z):
            continue
        exp = cp.predict_conv(k, n, x, y, z)
        nm = ('R4' if y == 0 else 'R5') + ('o' if exp[0] == 'Q' else '')
        if done[nm] >= N:
            continue
        check_macro(nm, ('Q', k, n, x, y, z), cp.CONV, (k, n, x, y, z), exp)
        done[nm] += 1
    cnt.update(done)
    print('   random instances per rule (k <= 10, blocks up to 200 or 2^(k-1)), all simulated on the twin:', cnt)


# ------------------------------------------------------------------------------------------------ I
BASE_STEP = 16590
BASE = (6, 14, 4)                      # S(6,14,4) = Q(7,65,1,14,4) = CfgF 5 0 14 4
PREFIX_STEP, PREFIX = 1711, ('Q', 5, 20, 3, 1, 0)


def section_I():
    print('[I] initial segment (exact): blank -> base')
    r = sim.run(MT, sim.Config.blank(), BASE_STEP)
    ok(not r.halted and r.step == BASE_STEP, 'run short / halted')
    got = cp.norm_cells(list(r.final.tape), r.final.head, r.final.state)
    want = cp.cfg_norm(tw(cp.build(cp.S(*BASE))))
    ok(got == want, f'blank -> S{BASE} at {BASE_STEP}')
    ok(cp.S(*BASE) == ('Q', 7, 65, 1, 14, 4), 'S(6,14,4) form')
    print(f'   c0 --{BASE_STEP:,}--> {tw(cp.build(cp.S(*BASE)))[:60]}...  = sigma Q(7,65,1,14,4): agrees')
    r = sim.run(MT, sim.Config.blank(), PREFIX_STEP)
    got = cp.norm_cells(list(r.final.tape), r.final.head, r.final.state)
    ok(got == cp.cfg_norm(tw(cp.build(PREFIX))), 'prefix config')
    T = cp.Tape(cp.build(PREFIX))
    c = cp.CONV(T, *PREFIX[1:]); seq = ['R5']
    while c[0] == 'P' and c[2] < 2 ** c[1] - 1:
        c = cp.R1(T, *c[1:]); seq.append('R1')
    c = cp.R2(T, *c[1:]); seq.append('R2')
    while c[0] == 'Q' and c[4] >= 2:
        c = cp.R3(T, *c[1:]); seq.append('R3')
    c = cp.CONV(T, *c[1:]); seq.append('R5' if c[0] == 'P' else 'R5o')
    while c[0] == 'P' and c[2] < 2 ** c[1] - 1:
        c = cp.R1(T, *c[1:]); seq.append('R1')
    c = cp.R2(T, *c[1:]); seq.append('R2')
    ok(c == cp.S(*BASE) and PREFIX_STEP + T.steps == BASE_STEP, f'prefix rules: {c} at {PREFIX_STEP + T.steps}')
    comp = []
    for s in seq:
        if comp and comp[-1][0] == s:
            comp[-1][1] += 1
        else:
            comp.append([s, 1])
    print(f'   also by the rules: c0 --1,711--> Q(5,20,3,1,0) (exact), then '
          + ', '.join(f'{s}^{n}' if n > 1 else s for s, n in comp)
          + f' --> S(6,14,4) at step {PREFIX_STEP + T.steps:,} (runner step sum)')


# ------------------------------------------------------------------------------------------------ C
def section_C(limit):
    """Every canonical right-edge visit of the twin (state B = sigma(C), parsed as P/Q) from step 1,711 to `limit`
    against the rule runner started at Q(5,20,3,1,0): same configuration, same step, same order."""
    print(f'[C] every canonical right-edge visit of the twin from step {PREFIX_STEP:,} to {limit:,}: rules vs simulation')
    sys.path.insert(0, os.path.join(HERE, 'work', 'scripts'))
    import canon
    real, other = [], 0
    for ev in sim.events(MT, max_steps=limit, on_record=True, snapshot=True):
        if ev.side != 'R' or ev.state != 'B':
            continue
        c = canon.canon(ev.config)
        if c is None:
            other += 1
        elif ev.step >= PREFIX_STEP:
            real.append((ev.step, c))
    pred = [(PREFIX_STEP, PREFIX)]
    T = cp.Tape(cp.build(PREFIX))
    v = lambda c, s: pred.append((PREFIX_STEP + s, c))
    c = cp.CONV(T, *PREFIX[1:]); v(c, T.steps)
    while c[2] < 2 ** c[1] - 1:
        c = cp.R1(T, *c[1:]); v(c, T.steps)
    c = cp.R2(T, *c[1:]); v(c, T.steps)
    k, Y, Z = 5, 5, 5
    ok(c == cp.S(k, Y, Z), 'prefix to S(5,5,5)')
    while PREFIX_STEP + T.steps <= limit:
        cp.cycle(T, k, Y, Z, visit=v)
        Y, Z = cp.F(k, Y, Z); k += 1
    ded = []
    for q in pred:
        if not ded or ded[-1] != q:
            ded.append(q)
    ded = [q for q in ded if q[0] <= limit]
    ok(real == ded, f'canonical visits differ: {len(real)} real vs {len(ded)} predicted')
    print(f'   {len(real)} of {len(ded)} canonical visits identical (configuration and step);'
          f' {other} other right-edge B visits (the start and mid-conversion edge touches)')


# ------------------------------------------------------------------------------------------------ B
def inv_rocq(k, i, y, z):
    """The parent's .v `Inv k i y z` (index coordinates), verbatim."""
    N = 2 ** (k + 1)
    return (64 <= N and N <= 8 * (y + 2 * i) and 8 * (y + 2 * i) <= 3 * N and 3 * z <= N + 3 * i
            and 8 * i <= 8 * z + N and N <= 4 * y + 6 * i + 2 * z and 20 * y + 24 * i + 16 * z <= 9 * N)


def hyps(k, Y, Z):
    N = 2 ** k
    Y2, Z2 = cp.F(k, Y, Z)
    return (Y + 2 * Z >= 1 and 3 * (Y // 2) <= N - 2 and Y2 >= 0 and Y2 + 2 * Z2 >= 1 and 3 * (Y2 // 2) <= 2 * N - 2
            and (Y % 2 == 1 or Z2 >= 1))


def section_B():
    print('[B] the base S(6,14,4) = CfgF 5 0 14 4')
    k, Y, Z = BASE
    N = 2 ** k
    ok(inv_rocq(5, 0, 14, 4), 'Rocq Inv 5 0 14 4')
    ok(cp.inR(k, Y, Z), 'base in N*R')
    ok(hyps(k, Y, Z), 'cycle-lemma hypotheses at the base')
    facets = [(a * Y + b * Z, h * N, f'{a}Y{b:+d}Z <= {h}N') for a, b, h in cp.POLY]
    tight = [f for f in facets if f[0] == f[1]]
    print('   facets (lhs, rhs):', [(f[2], f[0], f[1]) for f in facets])
    print('   tight:', [f[2] for f in tight], '(i.e. 4Y+2Z >= N with 4*14+2*4 = 64 = N)')
    ok([f[2] for f in tight] == ['-4Y-2Z <= -1N'], 'tight facet')
    print('   parent .v Inv 5 0 14 4 holds; in N*R; cycle-lemma (a)-(d) hold; even branch (Y = 14)')
    ok(not inv_rocq(5, 0, 13, 4) and not inv_rocq(5, 0, 14, 3), 'tightness controls')
    print('   negative controls: Inv 5 0 13 4 and Inv 5 0 14 3 fail (the base is on the facet)')


# ------------------------------------------------------------------------------------------------ O
def section_O(kmax):
    print(f'[O] twin orbit from the base: cycle map vs direct simulation of the twin, to S({kmax},...)')
    k, Y, Z = BASE
    cfg = sim.run(MT, sim.Config.blank(), BASE_STEP).final
    step = BASE_STEP
    use = {}
    path = [(k, Y, Z, step)]
    bounds = []                                    # phase boundaries (first/last config of each phase run)
    nb = 0
    while k < kmax:
        T = cp.Tape(cp.build(cp.S(k, Y, Z)))
        seen = []
        c = cp.cycle(T, k, Y, Z, visit=lambda c_, s_: seen.append((c_, s_)), use=use)
        Y2, Z2 = cp.F(k, Y, Z)
        ok(c == cp.S(k + 1, Y2, Z2), f'runner cycle {k}')
        ok(hyps(k, Y, Z), f'hypotheses at {(k, Y, Z)}')
        ph = lambda c_: (c_[0], c_[1])
        marks = [i for i in range(len(seen) - 1)
                 if i == 0 or ph(seen[i - 1][0]) != ph(seen[i][0]) or ph(seen[i + 1][0]) != ph(seen[i][0])]
        base_step = step
        for i in marks + [len(seen) - 1]:          # chained direct simulation through every boundary
            c_, s_ = seen[i]
            r = sim.run(MT, cfg, base_step + s_ - step)
            ok(not r.halted and r.step == base_step + s_, f'twin sim halted/short in cycle {k}: {r.step}')
            got = cp.norm_cells(list(r.final.tape), r.final.head, r.final.state)
            ok(got == cp.cfg_norm(tw(cp.build(c_))), f'twin sim != {c_} at {base_step + s_}')
            cfg, step = r.final, base_step + s_
            nb += 1
            if i < len(seen) - 1:
                bounds.append((step, c_))
        k, Y, Z = k + 1, Y2, Z2
        ok(cp.inR(k, Y, Z), f'{(k, Y, Z)} not in R')
        path.append((k, Y, Z, step))
    bounds.append((step, cp.S(k, Y, Z)))
    ok('R4o' not in use, f'R4o used: {use}')
    for k_, Y_, Z_, s_ in path:
        print(f'   S({k_},{Y_},{Z_}) at step {s_:,}  (in R: {cp.inR(k_, Y_, Z_)}, parity {"odd" if Y_ % 2 else "even"})')
    print(f'   {len(path) - 1} of {len(path) - 1} predicted cycle states, and all {nb} phase-boundary configurations'
          ' (phase starts/ends), hit at the exact predicted step; macro rules used:', dict(sorted(use.items())))
    # the parent's orbit for comparison (abstract only)
    k, Y, Z = 6, 16, 0
    par = [(k, Y, Z)]
    while k < kmax:
        Y, Z = cp.F(k, Y, Z); k += 1; par.append((k, Y, Z))
    print('   parent orbit (for comparison, from S(6,16,0)):', par)
    if not ARGS.quick:
        rows = [{'step': s_, 'values': [c_[0], c_[1], c_[2], c_[3], c_[4] if c_[0] == 'Q' else '-', c_[-1]]}
                for s_, c_ in bounds]
        with open(os.path.join(HERE, 'work', 'history_rows.json'), 'w') as f:
            json.dump(rows, f, indent=0)
        print(f'   wrote work/history_rows.json ({len(rows)} phase-boundary rows from step {BASE_STEP:,})')
    return path


if __name__ == '__main__':
    section_S()
    section_L(int(ARGS.lock if not ARGS.quick else 2e5))
    section_A()
    section_M()
    section_I()
    section_C(int(2e6 if ARGS.quick else 3e7))
    section_B()
    print('[H] parent section 5 (exact rational hexagon check), re-run unchanged')
    nf = len(cp.FAILS)
    cp.section5()
    ok(len(cp.FAILS) == nf, 'parent section 5')
    section_O(min(ARGS.kmax, 11) if ARGS.quick else ARGS.kmax)
    print()
    print('ALL CHECKS PASSED' if not FAILS else f'{len(FAILS)} FAILURES')
    sys.exit(1 if FAILS else 0)