#!/usr/bin/env python3
"""check_proof.py -- checks the claims of proof.md for
1RB0LA_0RC1RD_0LD1RE_1LA0RB_0RD1RF_1RC---  (run: python3 check_proof.py [--quick]).

What kind of check each section is:
  * EXACT (a proof of a finite claim): [1] the 21 shift rules, the polygon vertex computation in [5],
    and the initial segment in [6].
  * RANGE/SAMPLE (evidence only): [2]-[4] and the orbit in [6]. The parametric claims (CTR for all k, n;
    QC for every G; R1-R5o; the cycle lemma) are PROVED by the hand derivations in proof.md section 3-4;
    this script re-runs those derivations over ranges and random samples.

[1] shift rules: exact finite execution on the rule's window (head never leaves it, no halt).
    For a fixed-word rule this IS a proof (it is what `es` does); it also lists the transitions used.
[2] repeated rules (AL^a, ..., CTR, OV): exact windowed execution for small parameters.
[3] macro rules R1..R5o (R4o too, though the proof never uses it): a derivation runner applies the rule
    chain written in proof.md to the concrete tape (each application checks the rule's pattern on the
    real tape), sums the step counts, and the result is compared with the predicted configuration AND
    with direct simulation for that many steps.
[4] cycle map F_N: runner + direct simulation from random invariant points; exhaustive integer check of
    invariance, and of which conversion rule each invariant point uses (R4o must never occur).
[5] invariant polygon: exact rational check (vertex images, all N >= 64), side conditions, including
    max(2u+v) = 27/32 (the even conversion never overflows).
[6] initial segment and full orbit: runner from step 1507 reproduces every observe.json history row
    (configuration and exact step number) up to 2.3e9 steps; rule-use tally; direct simulation spot checks.
"""
import sys, os, random, json, itertools
from fractions import Fraction as Fr

HERE = os.path.dirname(os.path.abspath(__file__))
sys.path.insert(0, "[local working file]")
import sim  # noqa: E402

TM = "1RB0LA_0RC1RD_0LD1RE_1LA0RB_0RD1RF_1RC---"
M = sim.parse(TM)
QUICK = "--quick" in sys.argv
FAILS = []


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


# own transition table (independent of sim.py)
DELTA = {}
for i, blk in enumerate(TM.split("_")):
    for s in (0, 1):
        t = blk[3 * s:3 * s + 3]
        DELTA["ABCDEF"[i], s] = None if t == "---" else (int(t[0]), 1 if t[1] == "R" else -1, t[2])

# ---------------------------------------------------------------------------------------------
# Rules.  Notation as in proof.md:  `L <S R` head on the last cell of L;  `L S> R` head on first of R.
# ---------------------------------------------------------------------------------------------
RULES = [  # name, lhs, rhs, steps (claimed in proof.md)
    ("ENDL", "00 C> 0^inf", "<A 001 0^inf", 7),
    ("AL", "00111 <A 001", "<A 001 00101", 25),
    ("BL", "10111 <A 001", "<A 001 00001", 9),
    ("D1L", "1010 <A 001", "<A 001 0101", 12),
    ("FLIP", "0010 <A 001", "1010 001 E>", 18),
    ("D0R", "001 E> 0101", "0010 001 E>", 4),
    ("AR", "001 E> 00101", "00111 001 E>", 7),
    ("BR", "001 E> 00001", "10111 001 E>", 23),
    ("ENDR", "001 E> 0^inf", "10111 00 C> 0^inf", 22),
    ("FL", "00111 111 <A 001", "<A 001 001 00001", 28),
    ("FR1", "001 E> 001 00001", "00111 <A 001 001", 13),
    ("FR2", "001 E> 001 00101", "00111 <A 001 101", 13),
    ("FR3", "001 E> 101 00101", "00111 111 001 E>", 10),
    ("FR4", "001 E> 101 00001", "<A 001 001 00101", 36),
    ("FE0", "001 E> 001 0^inf", "00111 00 C> 0^inf", 6),
    ("FE1", "001 E> 101 0^inf", "<A 001 001 001 0^inf", 36),
    ("FE2", "001 E> 001 001 0^inf", "00111 <A 001 1 0^inf", 13),
    ("FE3", "001 E> 1 0^inf", "<A 001 0011 0^inf", 36),
    ("FE4", "001 E> 0011 0^inf", "00111 <A 001 001 0^inf", 25),
    ("OVL", "0^inf <A 001 01", "0^inf 1010 001 E>", 12),
    ("OVR", "001 E> 01 00101", "0010 <A 001 101", 10),
]


def parse_side(cfg):
    """-> (window bits, head offset in window (may be -1 or len), state, linf, rinf)."""
    toks = cfg.split()
    linf, rinf = toks[0] == "0^inf", toks[-1] == "0^inf"
    toks = toks[1 if linf else 0: len(toks) - (1 if rinf else 0)]
    h = [i for i, t in enumerate(toks) if t[0] == "<" or t[-1] == ">"]
    assert len(h) == 1, cfg
    h = h[0]
    L, R = "".join(toks[:h]), "".join(toks[h + 1:])
    st = toks[h].strip("<>")
    off = len(L) - 1 if toks[h][0] == "<" else len(L)
    return [int(c) for c in L + R], off, st, linf, rinf


def norm_rule(lhs, rhs):
    a, b = parse_side(lhs), parse_side(rhs)
    wa, oa, sa, la, ra = a
    wb, ob, sb, lb, rb = b
    assert (la, ra) == (lb, rb), (lhs, rhs)
    n = max(len(wa), len(wb))
    if ra:  # pad with zeros on the infinite right side
        wa, wb = wa + [0] * (n - len(wa)), wb + [0] * (n - len(wb))
    if la:  # pad on the infinite left side
        oa += n - len(wa); ob += n - len(wb)
        wa, wb = [0] * (n - len(wa)) + wa, [0] * (n - len(wb)) + wb
    assert len(wa) == len(wb), (lhs, rhs)
    return dict(wa=wa, oa=oa, sa=sa, wb=wb, ob=ob, sb=sb, linf=la, rinf=ra)


NR = {name: norm_rule(l, r) for name, l, r, _ in RULES}
STEPS = {name: s for name, _, _, s in RULES}


def window_exec(r, maxsteps=10 ** 5):
    """Run the machine on the rule's lhs window. Cells outside the window may be read only on an
    infinite (all-zero) side. Returns (steps, set of transitions used) when the rhs is reached."""
    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 = DELTA[st, sym]
        used.add(st + str(sym))
        if tr is None:
            return "HALT", used
        tape[pos] = tr[0]; pos += tr[1]; st = tr[2]
    return None, used


# ---------------------------------------------------------------------------------------------
# Derivation runner on a concrete tape
# ---------------------------------------------------------------------------------------------
class Tape:
    MARGIN = 64

    def __init__(self, cfg):
        w, off, st, linf, rinf = parse_side(cfg)
        assert linf and rinf
        self.t = bytearray(self.MARGIN) + bytearray(w) + bytearray(self.MARGIN)
        self.pos = self.MARGIN + off
        self.st = st
        self.steps = 0
        self.log = []

    def _room(self, lo, hi):
        if lo < 0:
            k = -lo + self.MARGIN
            self.t[0:0] = bytearray(k); self.pos += k; lo += k; hi += k
        if hi > len(self.t):
            self.t.extend(bytearray(hi - len(self.t) + self.MARGIN))
        return lo, hi

    def apply(self, name, count=1):
        if count < 0:
            raise ValueError(f"{name}^{count}")
        r = NR[name]
        for _ in range(count):
            n = len(r["wa"])
            start = self.pos - r["oa"]
            start, end = self._room(start, start + n)
            if self.st != r["sa"] or self.t[start:end] != bytes(r["wa"]) \
                    or (r["linf"] and any(self.t[:start])) or (r["rinf"] and any(self.t[end:])):
                raise ValueError(f"rule {name} does not match at step {self.steps}: {self.show()}")
            self.t[start:end] = bytes(r["wb"])
            self.pos, self.st = start + r["ob"], r["sb"]
            self.steps += STEPS[name]
        self.log.append((name, count))

    # repeated block rules (the L1 lemmas), applied in one go after checking the whole pattern
    def rep(self, name, count, blk_in, blk_out, side):
        """AL/BL (side 'L':  (in)^c <A 001 -> <A 001 (out)^c) or AR/BR (side 'R')."""
        if count == 0:
            return
        L = len(blk_in)
        if side == "L":
            # head A on last cell of the blocks, 001 to the right of the head
            s = self.pos + 1 - L * count
            s, _ = self._room(s, self.pos + 4)
            e = self.pos + 1
            if self.st != "A" or self.t[s:e] != bytes(blk_in) * count or self.t[e:e + 3] != b"\0\0\1":
                raise ValueError(f"{name}^{count} does not match: {self.show()}")
            self.t[s:e + 3] = b"\0\0\1" + bytes(blk_out) * count
            self.pos = s - 1
        else:
            s = self.pos - 3
            s, e = self._room(s, self.pos + L * count)
            if self.st != "E" or self.t[s:s + 3] != b"\0\0\1" or self.t[s + 3:e] != bytes(blk_in) * count:
                raise ValueError(f"{name}^{count} does not match: {self.show()}")
            self.t[s:e] = bytes(blk_out) * count + b"\0\0\1"
            self.pos = e
        self.steps += STEPS[name] * count
        self.log.append((name, count))

    def AL(self, a): self.rep("AL", a, [0, 0, 1, 1, 1], [0, 0, 1, 0, 1], "L")
    def BL(self, a): self.rep("BL", a, [1, 0, 1, 1, 1], [0, 0, 0, 0, 1], "L")
    def AR(self, a): self.rep("AR", a, [0, 0, 1, 0, 1], [0, 0, 1, 1, 1], "R")
    def BR(self, a): self.rep("BR", a, [0, 0, 0, 0, 1], [1, 0, 1, 1, 1], "R")

    def norm(self):
        return norm_cells(list(self.t), self.pos, self.st)

    def show(self):
        a, h, st = self.norm()
        s = "".join(map(str, a))
        return f"{s[:h]}[{st}]{s[h:]}"


def norm_cells(cells, pos, st):
    a, b = 0, len(cells)
    while a < pos and cells[a] == 0: a += 1
    while b - 1 > pos and cells[b - 1] == 0: b -= 1
    return tuple(cells[a:b]), pos - a, st


def cfg_norm(cfgstr):
    w, off, st, _, _ = parse_side(cfgstr)
    cells = [0] * 8 + w + [0] * 8
    return norm_cells(cells, off + 8, st)


def sim_norm(cfgstr, steps):
    """Direct simulation (tools/sim.py, C core for long runs) for exactly `steps` steps."""
    r = sim.run(M, 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 norm_cells(list(c.tape), c.head, c.state)


# ---------------------------------------------------------------------------------------------
# Families
# ---------------------------------------------------------------------------------------------
def C(k, n):
    assert 2 ** (k - 1) <= n < 2 ** k, (k, n)
    return " ".join("1010" if b == "1" else "0010" for b in format(n, "b"))


def P(k, n, X, z):
    return f"0^inf {C(k, n)} " + "00111 " * X + "10111 " * z + "00 C> 0^inf"


def Q(k, n, x, y, z):
    return f"0^inf {C(k, n)} " + "00111 " * x + "111 " + "00111 " * y + "10111 " * z + "00 C> 0^inf"


def tz1(n):  # trailing 1-bits
    t = 0
    while n & 1: n >>= 1; t += 1
    return t


# ---------------------------------------------------------------------------------------------
# Derivations (exactly the chains written in proof.md section 3)
# ---------------------------------------------------------------------------------------------
def CTR(T, k, n):  # 0^inf C_k(n) <A 001 -> 0^inf C_k(n+1) 001 E>,  n+1 < 2^k
    assert n + 1 < 2 ** k
    t = tz1(n)
    T.apply("D1L", t); T.apply("FLIP"); T.apply("D0R", t)


def OV(T, k):  # 0^inf C_k(2^k-1) <A 001 00101 00101 -> 0^inf C_{k+1}(2^k+1) 00111 111 001 E>
    T.apply("D1L", k); T.apply("OVL"); T.apply("D0R", k - 1); T.apply("OVR"); T.apply("FLIP"); T.apply("FR3")


def QC(T, k, n, x, y):
    """C_k(n) (00111)^x 111 (00111)^y <A 001 G  ->  C_k(n+3) (00111)^(x+2) 111 (00111)^(y-2) 001 E> G"""
    assert x >= 1 and y >= 2 and n + 3 <= 2 ** k - 1
    T.AL(y); T.apply("FL"); T.AL(x - 1); CTR(T, k, n); T.AR(x - 1); T.apply("FR1")
    T.AL(x); CTR(T, k, n + 1); T.AR(x); T.apply("FR2")
    T.AL(x + 1); CTR(T, k, n + 2); T.AR(x + 1); T.apply("FR3"); T.AR(y - 2)


def R1(T, k, n, X, z):
    T.apply("ENDL"); T.BL(z); T.AL(X); CTR(T, k, n); T.AR(X); T.BR(z); T.apply("ENDR")
    return ("P", k, n + 1, X, z + 1)


def R2(T, k, n, X, z):
    assert n == 2 ** k - 1 and X >= 2
    T.apply("ENDL"); T.BL(z); T.AL(X); OV(T, k); T.AR(X - 2); T.BR(z); T.apply("ENDR")
    return ("Q", k + 1, 2 ** k + 1, 1, X - 2, z + 1)


def R3(T, k, n, x, y, z):
    T.apply("ENDL"); T.BL(z); QC(T, k, n, x, y); T.BR(z); T.apply("ENDR")
    return ("Q", k, n + 3, x + 2, y - 2, z + 1)


DELTA_A = {"FR1": 1, "FR2": 1, "FR4": 0, "FE1": 0, "FE2": 1, "FE3": 0, "FE4": 1}


def actions(y, z):
    if y == 0:
        return ["FR1"] * (z + 1) + ["FE0"]
    return ["FR1", "FR2"] + ["FR4", "FR2"] * z + ["FE1", "FE2", "FE3", "FE4", "FE0"]


def CONV(T, k, n, x, y, z):
    """R4 / R5 / R4o / R5o: Q(k,n,x,y,z), y in {0,1} -> next right-edge configuration."""
    assert y in (0, 1) and x >= 1
    T.apply("ENDL"); T.BL(z); T.AL(y); T.apply("FL")
    acts = actions(y, z)
    a = x - 1          # A blocks between counter and old front (hL is at the counter after AL^a)
    T.AL(a)
    over = False       # after the overflow: new front present
    xn = yb = None     # new front: x' A blocks before it, yb A blocks between it and the old front
    for r, act in enumerate(acts, 1):
        if not over:
            if n < 2 ** k - 1:
                CTR(T, k, n); n += 1; T.AR(a)
            else:  # overflow round
                OV(T, k); n = 2 ** k + 1; k += 1; over = True; xn = 1
                T.AR(a - 2); yb = a - 2
        else:
            QC(T, k, n, xn, yb); n += 3; xn += 2; yb -= 2
        T.apply(act)
        if act == "FE0":
            break
        d = DELTA_A[act]
        if not over:
            a += d; T.AL(a)
        else:
            yb += d  # the hL is now right of the between region; QC starts with AL^yb
    if not over:
        return ("P", k, n, a + 1, 0)
    return ("Q", k, n, xn, yb + 1, 0)


def predict_conv(k, n, x, y, z):
    """The closed forms of proof.md: R4, R5, R4o, R5o."""
    T_ = z + 2 if y == 0 else 2 * z + 7
    X = x + z + 1 if y == 0 else x + z + 4
    m = n + T_ - 2 ** k
    if m < 0:
        return ("P", k, n + T_, X, 0)
    return ("Q", k + 1, 2 ** k + 1 + 3 * m, 1 + 2 * m, X - 2 - 2 * m, 0)


def build(c):
    return P(*c[1:]) if c[0] == "P" else Q(*c[1:])


# ---------------------------------------------------------------------------------------------
def section1():
    print("[1] shift rules: exact windowed execution")
    allused = set()
    for name, lhs, rhs, st in RULES:
        steps, used = window_exec(NR[name])
        allused |= used
        ok(steps == st, f"{name}: {lhs} -> {rhs}: got {steps}, claimed {st}")
        ok("F1" not in used, f"{name} reads F1")
        print(f"   {name:5s} {lhs:22s} -> {rhs:24s} {steps!s:>4} steps  F used: {sorted(u for u in used if u[0]=='F')}")
    print("   transitions used by the rule set:", sorted(allused), "(F1 never)")


def window_rule(lhs, rhs):
    return window_exec(norm_rule(lhs, rhs), 10 ** 7)[0]


def section2():
    print("[2] repeated rules (exact windowed execution)")
    A = 12 if QUICK else 25
    for a in range(A):
        ok(window_rule("00111 " * a + "<A 001", "<A 001 " + "00101 " * a) == 25 * a, f"AL^{a}")
        ok(window_rule("10111 " * a + "<A 001", "<A 001 " + "00001 " * a) == 9 * a, f"BL^{a}")
        ok(window_rule("001 E> " + "00101 " * a, "00111 " * a + "001 E>") == 7 * a, f"AR^{a}")
        ok(window_rule("001 E> " + "00001 " * a, "10111 " * a + "001 E>") == 23 * a, f"BR^{a}")
    print(f"   AL^a, BL^a, AR^a, BR^a: a = 0..{A-1}: steps 25a, 9a, 7a, 23a")
    K = 8 if QUICK else 10
    cnt = 0
    for k in range(1, K + 1):
        for n in range(2 ** (k - 1), 2 ** k - 1):
            s = window_rule(f"0^inf {C(k, n)} <A 001", f"0^inf {C(k, n + 1)} 001 E>")
            ok(s == 16 * tz1(n) + 18, f"CTR k={k} n={n}: {s}")
            cnt += 1
        s = window_rule(f"0^inf {C(k, 2 ** k - 1)} <A 001 00101 00101", f"0^inf {C(k + 1, 2 ** k + 1)} 00111 111 001 E>")
        ok(s == 16 * k + 46, f"OV k={k}: {s}")
    print(f"   CTR(k,n) for all {cnt} pairs with k<={K}, n+1<2^k: 16 t(n)+18 steps;  OV(k), k=1..{K}: 16k+46 steps")


def run_check(name, start, fn, args, expect, simulate=True):
    T = Tape(build(start))
    try:
        got = fn(T, *args)
    except (ValueError, AssertionError) as e:
        ok(False, f"{name} {start}: derivation failed: {e}")
        return None
    ok(got == expect, f"{name} {start}: derived {got}, expected {expect}")
    ok(T.norm() == cfg_norm(build(got)), f"{name} {start}: runner tape != {got}")
    if simulate:
        ok(sim_norm(build(start), T.steps) == cfg_norm(build(got)), f"{name} {start}: simulation disagrees after {T.steps}")
    return T.steps


def in_domain_conv(k, n, x, y, z):
    """Domain of R4 / R5 / R4o / R5o as stated in proof.md."""
    if y not in (0, 1) or x < 1 or not 2 ** (k - 1) <= n <= 2 ** k - 1:
        return False
    c = predict_conv(k, n, x, y, z)
    if c[0] == "P":
        return True
    m = (c[2] - 2 ** k - 1) // 3
    return c[4] >= 1 and 3 * m <= 2 ** k - 2


def section3():
    print("[3] macro rules: derivation runner (+ direct simulation)")
    stats = {}

    def tally(nm):
        stats[nm] = stats.get(nm, 0) + 1
    # (a) exhaustive over small parameters; every 3rd instance also by direct simulation
    i = 0
    K = 5 if QUICK else 6
    B = 5 if QUICK else 7
    for k in range(1, K + 1):
        for n in range(2 ** (k - 1), 2 ** k):
            for X in range(0, B):
                for z in range(0, B):
                    i += 1
                    if n < 2 ** k - 1:
                        s = run_check("R1", ("P", k, n, X, z), R1, (k, n, X, z), ("P", k, n + 1, X, z + 1), simulate=i % 3 == 0)
                        ok(s == 32 * X + 32 * z + 16 * tz1(n) + 47, f"R1 steps {s}"); tally("R1 small")
                    elif X >= 2:
                        s = run_check("R2", ("P", k, n, X, z), R2, (k, n, X, z), ("Q", k + 1, 2 ** k + 1, 1, X - 2, z + 1), simulate=i % 3 == 0)
                        ok(s == 32 * X + 32 * z + 16 * k + 61, f"R2 steps {s}"); tally("R2 small")
            for x in range(1, B):
                for y in range(0, B + 1):
                    for z in range(0, B):
                        i += 1
                        if y >= 2:
                            if n + 3 <= 2 ** k - 1:
                                s = run_check("R3", ("Q", k, n, x, y, z), R3, (k, n, x, y, z), ("Q", k, n + 3, x + 2, y - 2, z + 1), simulate=i % 3 == 0)
                                ok(s == 96 * x + 32 * (y - 2) + 32 * z + 16 * (tz1(n) + tz1(n + 1) + tz1(n + 2)) + 197, f"R3 steps {s}")
                                tally("R3 small")
                        elif in_domain_conv(k, n, x, y, z):
                            exp = predict_conv(k, n, x, y, z)
                            nm = ("R4" if y == 0 else "R5") + ("o" if exp[0] == "Q" else "")
                            run_check(nm, ("Q", k, n, x, y, z), CONV, (k, n, x, y, z), exp, simulate=i % 3 == 0)
                            tally(nm + " small")
    # (b) random instances with large parameters, all by direct simulation too
    rnd = random.Random(2026)
    N = 15 if QUICK else 60

    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, 11); n = rnd.randint(2 ** (k - 1), 2 ** k - 2); X = pick(0, 300); z = pick(0, 300)
        s = run_check("R1", ("P", k, n, X, z), R1, (k, n, X, z), ("P", k, n + 1, X, z + 1))
        ok(s == 32 * X + 32 * z + 16 * tz1(n) + 47, f"R1 steps {s}"); tally("R1 large")
        k = rnd.randint(1, 11); X = pick(2, 300); z = pick(0, 300)
        s = run_check("R2", ("P", k, 2 ** k - 1, X, z), R2, (k, 2 ** k - 1, X, z), ("Q", k + 1, 2 ** k + 1, 1, X - 2, z + 1))
        ok(s == 32 * X + 32 * z + 16 * k + 61, f"R2 steps {s}"); tally("R2 large")
        k = rnd.randint(3, 11); n = rnd.randint(2 ** (k - 1), 2 ** k - 4); x = pick(1, 300); y = pick(2, 300); z = pick(0, 300)
        s = run_check("R3", ("Q", k, n, x, y, z), R3, (k, n, x, y, z), ("Q", k, n + 3, x + 2, y - 2, z + 1))
        ok(s == 96 * x + 32 * (y - 2) + 32 * z + 16 * (tz1(n) + tz1(n + 1) + tz1(n + 2)) + 197, f"R3 steps {s}")
        tally("R3 large")
    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, 11); 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 in_domain_conv(k, n, x, y, z):
            continue
        exp = 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
        run_check(nm, ("Q", k, n, x, y, z), CONV, (k, n, x, y, z), exp)
        done[nm] += 1
        tally(nm + " large")
    print("   instances (all passed unless FAIL above):", stats)
    # (c) informational: just outside the domain the derivation breaks (the rules are not claimed there)
    out = tot = 0
    for k in range(3, 7):
        for n in range(2 ** (k - 1), 2 ** k):
            for x in range(1, 5):
                for y in (0, 1):
                    for z in range(0, 5):
                        if not in_domain_conv(k, n, x, y, z):
                            tot += 1
                            T = Tape(build(("Q", k, n, x, y, z)))
                            try:
                                CONV(T, k, n, x, y, z)
                            except (ValueError, AssertionError):
                                out += 1
    print(f"   (info) small y<=1 instances outside the R4o/R5o domain: the chain breaks on {out} of {tot} (domain is sufficient, not claimed necessary)")


# ---------------------------------------------------------------------------------------------
# The cycle map
# ---------------------------------------------------------------------------------------------
def S(k, Y, Z):
    """Start of the level-(k+1) cycle.  N = 2^k."""
    N = 2 ** k
    m = max(0, -Z)
    return ("Q", k + 1, N + 1 + 3 * m, 1 + 2 * m, Y - 2 * m, Z + m)


def F(k, Y, Z):
    N = 2 ** k
    if Y % 2 == 0:
        j = Y // 2
        return 3 * j + Z, N - 3 - 4 * j - Z
    j = (Y - 1) // 2
    return 3 * j + 3 + Z, N - 8 - 5 * j - 2 * Z


def conv_kind(k, Y, Z):
    """Which conversion rule the cycle lemma uses from S(k,Y,Z): R4, R4o, R5 or R5o."""
    N = 2 ** k
    j = Y // 2
    n, z = N + 1 + 3 * j, Z + j
    s = n + (z + 2 if Y % 2 == 0 else 2 * z + 7)
    return ("R4" if Y % 2 == 0 else "R5") + ("o" if s >= 2 * N else "")


def cycle(T, k, Y, Z, visit=None, use=None):
    """S(k,Y,Z) -> S(k+1, F(k,Y,Z)) with the rules; returns the end state. visit(cfg, steps) at every
    right-edge configuration; use (a dict) tallies the macro rules applied."""
    def tal(nm, c=1):
        if use is not None: use[nm] = use.get(nm, 0) + c
    c = S(k, Y, Z)
    if visit: visit(c, T.steps)
    while c[0] == "Q" and c[4] >= 2:
        c = R3(T, *c[1:]); tal("R3")
        if visit: visit(c, T.steps)
    y = c[4]
    c = CONV(T, *c[1:])
    tal(("R4" if y == 0 else "R5") + ("o" if c[0] == "Q" else ""))
    if visit: visit(c, T.steps)
    if c[0] == "P":
        while c[2] < 2 ** c[1] - 1:
            c = R1(T, *c[1:]); tal("R1")
            if visit: visit(c, T.steps)
        c = R2(T, *c[1:]); tal("R2")
        if visit: visit(c, T.steps)
    return c


# invariant polygon  R = { (u,v) : a u + b v <= h }  with u = Y/N, v = Z/N
POLY = [(-8, 0, -1), (8, 0, 3), (0, 3, 1), (0, -8, 1), (-4, -2, -1), (20, 16, 9)]  # a*Y + b*Z <= h*N


def inR(k, Y, Z):
    N = 2 ** k
    return all(a * Y + b * Z <= h * N for a, b, h in POLY)


def section4():
    print("[4] cycle map: S(k,Y,Z) -> S(k+1,F_N(Y,Z)) by the rules; direct simulation for k <= 8")
    rnd = random.Random(7)
    cnt = {}
    rk = {}
    for k in ([6, 7, 8] if QUICK else [6, 7, 8, 9, 10, 11]):
        N = 2 ** k
        done = 0
        while done < (6 if QUICK else 12):
            Y = rnd.randint(N // 8, 3 * N // 8); Z = rnd.randint(-N // 8, N // 3)
            if not inR(k, Y, Z):
                continue
            T = Tape(build(S(k, Y, Z)))
            u = {}
            try:
                got = cycle(T, k, Y, Z, use=u)
            except (ValueError, AssertionError) as e:
                ok(False, f"cycle k={k} Y={Y} Z={Z}: {e}")
                continue
            Y2, Z2 = F(k, Y, Z)
            kind = [nm for nm in u if nm[:2] in ("R4", "R5")]
            ok(kind == [conv_kind(k, Y, Z)] and "R4o" not in u, f"cycle k={k} Y={Y} Z={Z}: conversion {u}")
            for nm in u: rk[nm] = rk.get(nm, 0) + u[nm]
            ok(got == S(k + 1, Y2, Z2), f"cycle k={k} Y={Y} Z={Z}: got {got}, F gives {S(k+1, Y2, Z2)}")
            ok(inR(k + 1, Y2, Z2), f"invariant not preserved at k={k} Y={Y} Z={Z}")
            if k <= 8:
                ok(sim_norm(build(S(k, Y, Z)), T.steps) == cfg_norm(build(got)), f"cycle sim k={k} Y={Y} Z={Z}")
            done += 1
        cnt[k] = done
    print("   random invariant points per level k (N=2^k):", cnt, "-- all mapped by F_N, into R, sim agrees for k<=8")
    print("   macro rules used by those cycles:", dict(sorted(rk.items())))
    # exhaustive integer check of invariance (complements the exact rational proof of section 5), of the
    # cycle-lemma hypotheses (a)-(d), and of which conversion rule each point uses
    tot = 0
    kinds = {}
    for k in range(6, 10 if QUICK else 12):
        N = 2 ** k
        for Y in range(N // 8, 3 * N // 8 + 1):
            for Z in range(-N // 8, N // 3 + 1):
                if inR(k, Y, Z):
                    tot += 1
                    Y2, Z2 = F(k, Y, Z)
                    if not inR(k + 1, Y2, Z2):
                        ok(False, f"F leaves R: k={k} Y={Y} Z={Z}")
                    hyp = (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 (2 * Y + Z <= N - 4 and Z2 >= 1)))
                    if not hyp:
                        ok(False, f"cycle-lemma hypotheses fail: k={k} Y={Y} Z={Z}")
                    kd = conv_kind(k, Y, Z)
                    kinds[kd] = kinds.get(kd, 0) + 1
    ok("R4o" not in kinds, f"an invariant point uses R4o: {kinds}")
    print(f"   exhaustive: all {tot} integer points of N*R for N=64..{2**(9 if QUICK else 11)} map into 2N*R"
          f" and satisfy cycle-lemma hypotheses (a)-(d)")
    print(f"   conversion rule used, over those points: {dict(sorted(kinds.items()))} (R4o never)")


def section5():
    print("[5] invariant polygon R (exact rationals)")
    C_ = [(Fr(a), Fr(b), Fr(h)) for a, b, h in POLY]
    vs = []
    for (a1, b1, h1), (a2, b2, h2) in itertools.combinations(C_, 2):
        d = a1 * b2 - a2 * b1
        if d == 0: continue
        u = (h1 * b2 - h2 * b1) / d; v = (a1 * h2 - a2 * h1) / d
        if all(a * u + b * v <= h for a, b, h in C_): vs.append((u, v))
    vs = sorted(set(vs))
    print("   vertices:", ", ".join(f"({u},{v})" for u, v in vs))
    slacks = []
    for (u, v) in vs:
        for par in (0, 1):
            for t in (Fr(0), Fr(1, 64)):  # t = 1/N, affine in t: endpoints suffice for all N >= 64
                if par == 0:
                    u2, v2 = Fr(3, 4) * u + Fr(1, 2) * v, Fr(1, 2) - u - Fr(1, 2) * v - Fr(3, 2) * t
                else:
                    u2, v2 = Fr(3, 4) * u + Fr(1, 2) * v + Fr(3, 4) * t, Fr(1, 2) - Fr(5, 4) * u - v - Fr(11, 4) * t
                for a, b, h in C_:
                    slacks.append((h - (a * u2 + b * v2), (u, v), par, t, (int(a), int(b), int(h))))
    m = min(s[0] for s in slacks)
    ok(m > 0, f"polygon not invariant: {m}")
    ties = [s for s in slacks if s[0] == m]
    print(f"   min slack of image vertices over both parities and all N>=64: {m} (>0), attained "
          + "; ".join(f"at vertex ({s[1][0]},{s[1][1]}) parity {'odd' if s[2] else 'even'} 1/N={s[3]} facet {s[4][0]}Y{s[4][1]:+d}Z<={s[4][2]}N"
                      for s in ties))
    mins = min(u + 2 * v for u, v in vs); maxu = max(u for u, v in vs); minu = min(u for u, v in vs)
    ok(mins > 0 and minu > 0 and Fr(3, 2) * maxu < 1, "side conditions")
    print(f"   side conditions on R: min u = {minu} > 0;  min (u+2v) = {mins} (so Y+2Z >= N/16 >= 4);  max 3u/2 = {Fr(3,2)*maxu} < 1")
    # (d): the even conversion never overflows: 2Y+Z <= 27N/32 <= N-4 for N >= 64
    max2uv = max(2 * u + v for u, v in vs)
    arg = ", ".join(f"({u},{v})" for u, v in vs if 2 * u + v == max2uv)
    ok(max2uv == Fr(27, 32) and (1 - max2uv) * 64 >= 4, f"max 2u+v = {max2uv}")
    print(f"   max (2u+v) on R = {max2uv} at {arg}: 2Y+Z <= 27N/32 <= N-4 for N >= 64, so the even branch takes R4 (Z' >= 5N/32-3 > 0)")


def section6():
    print("[6] initial segment and full orbit")
    # blank tape -> Q(5,17,1,2,2) at 1507 and -> S(6,16,0) = Q(7,65,1,16,0) at 17215 by direct simulation
    for steps, c in [(1507, ("Q", 5, 17, 1, 2, 2)), (17215, S(6, 16, 0))]:
        r = sim.run(M, sim.Config.blank(), steps)
        got = norm_cells(list(r.final.tape), r.final.head, r.final.state)
        ok(got == cfg_norm(build(c)), f"c0 -> {c} at {steps}")
        print(f"   c0 --{steps}--> {c}: simulation agrees")
    ok(inR(6, 16, 0), "S(6,16,0) in R")
    print(f"   S(6,16,0) = Q(7,65,1,16,0) is in N*R for N=64:", inR(6, 16, 0))
    # rules from Q(5,17,1,2,2) (step 1507) reproduce the history rows of observe.json (config + step)
    hist = json.load(open(os.path.join(HERE, "observe.json")))["history"]["rows"]
    want = {}
    for row in hist:
        ph, k, n, x, y, z = row["values"]
        c = ("P", k, n, x, z) if ph == "P" else ("Q", k, n, x, y, z)
        want[c] = row["step"]
    seen = {}
    T = Tape(Q(5, 17, 1, 2, 2))
    base = 1507

    def visit(c, s):
        if c in want:
            seen[c] = base + s
    kfull = 11 if QUICK else 14   # full cycles k = 4 .. kfull-1, i.e. counter lengths 5 .. kfull
    k, Y, Z = 4, 2, 2
    path = []
    use = {}
    while k < kfull:
        c = cycle(T, k, Y, Z, visit, use)
        Y, Z = F(k, Y, Z); k += 1
        path.append((k, Y, Z, inR(k, Y, Z)))
        ok(c == S(k, Y, Z), f"orbit cycle {k}")
    if not QUICK:  # the last history row (P start at 2.3e9) lies inside the next cycle: run its Q phase + conversion
        c = S(k, Y, Z)
        while c[0] == "Q" and c[4] >= 2:
            c = R3(T, *c[1:]); visit(c, T.steps); use["R3"] = use.get("R3", 0) + 1
        y = c[4]
        c = CONV(T, *c[1:]); visit(c, T.steps)
        nm = ("R4" if y == 0 else "R5") + ("o" if c[0] == "Q" else "")
        use[nm] = use.get(nm, 0) + 1
    ok("R4o" not in use, f"orbit uses R4o: {use}")
    matched = [c for c in want if c in seen and seen[c] == want[c]]
    early = [c for c in want if want[c] < base]
    reach = base + T.steps
    missing = [c for c in want if c not in early and c not in matched and want[c] <= reach]
    ok(not missing, f"history rows not reproduced: {missing}")
    print(f"   orbit levels (k, Y, Z, in R): {path}")
    print(f"   macro rules applied by the runner from step 1507: {dict(sorted(use.items()))}")
    print(f"   history rows reproduced with exact step numbers: {len(matched)} of {len(want)}"
          f" ({len(early)} rows precede step 1507; rows beyond step {reach:,} not reached in this mode: {len(want)-len(matched)-len(early)})")
    print(f"   runner reached step {base + T.steps:,} applying only proof.md rules")


if __name__ == "__main__":
    section1(); section2(); section3(); section4(); section5(); section6()
    print()
    print("ALL CHECKS PASSED" if not FAILS else f"{len(FAILS)} FAILURES")
    sys.exit(1 if FAILS else 0)