#!/usr/bin/env python3
"""Checks for proof.md of 1RB0RA_1LC0RE_0LC0LD_1RD1LB_0RF1RA_0RB--- (twin #2 of pilot #2).

The parent 1RB0RF_1LC0RE_0LC0LD_1RD1LB_0RF1RA_0RB--- (KILLED) differs only in A1 (0RF there, 0RA here).
Every check prints `name: passed/total` (or a count). The run ends with ALL PASS or a list of failures.
Simulation is evidence, not proof. The `--rocq` part is a proof check (Coq 8.20.1, unmodified busycoq).

Notation (proof.md). Tokens: L=110, 0=000, 1=100 (3 cells) and the stray S=1 (1 cell, NEW). <w> = w's cells.
S(w) := 0^inf 11 D> <w> 0^inf. A tape is written canonically with the fewest strays (`parse`).
Passes over a right part starting at cell o (nothing left of the exit cell is read):
  h  : `B>` at cell o, exit `<C` on cell o-1.     h' : `1 D>` with the 1 at o-1, exit `<C` on cell o-2.
R(tp,k,w): hx = passes [h'] + h^k over <w>; h1 = passes h^k over `1 <w>` (exactly TwinSpec.v below).
TwinSpec.v also has RIncs_eq (re-parse an input up to tape equality); here that is `canon`, applied where proof.md
says (Lemma C: once, after the eater's leading L-run).

Run:  python3 -B check_proof.py          (default, ~25 s: from blank to 1e8 steps, small samples; checks 8a-8f = proof.md 3.3,
                                         9a-9j = pass 3: review gaps and the V* reformulation of 3.3.5)
      python3 -B check_proof.py --full   (~10 min: from blank to 1e9, bigger samples, abstract orbit 3000 sweeps, Rocq)
      python3 -B check_proof.py --rocq   (only the Rocq check, ~4 s)
"""
import itertools
import os
import random
import re
import subprocess
import sys
import tempfile
from collections import Counter

sys.path.insert(0, "[local working file]")
import sim  # noqa: E402

sys.setrecursionlimit(1_000_000)
TM = "1RB0RA_1LC0RE_0LC0LD_1RD1LB_0RF1RA_0RB---"
PARENT = "1RB0RF_1LC0RE_0LC0LD_1RD1LB_0RF1RA_0RB---"
TM7 = "1RB0RC_1LC0RE_0LC0LD_1RD1LB_0RF1RA_0RB---"   # SOMDCv6.v Module TM7
INV = {"L": "110", "0": "000", "1": "100", "S": "1"}
FULL = "--full" in sys.argv
ROCQ_ONLY = "--rocq" in sys.argv
ROOT = "[local working file]"
FAILS = []


def report(name, ok, tot, extra=""):
    print(f"{name}: {ok}/{tot}{('  ' + extra) if extra else ''}", flush=True)
    if ok != tot:
        FAILS.append(name)


def cells(w):
    return "".join(INV[t] for t in w)


def norm(c):
    return c.rstrip("0")


def parse(x):
    """Canonical token word of the cell string x (trailing zeros ignored): fewest strays S."""
    x = norm(x)
    n, INF = len(x), 10 ** 9
    cost, choice = [INF] * (n + 1), [None] * (n + 1)
    cost[n] = 0
    for i in range(n - 1, -1, -1):
        for t in "L01":
            w, seg = INV[t], x[i:i + 3]
            if len(seg) < 3:
                if w.startswith(seg) and set(w[len(seg):]) <= {"0"} and 0 < cost[i]:
                    cost[i], choice[i] = 0, (t, n)
            elif seg == w and cost[i + 3] < cost[i]:
                cost[i], choice[i] = cost[i + 3], (t, i + 3)
        if x[i] == "1" and cost[i + 1] + 1 < cost[i]:
            cost[i], choice[i] = cost[i + 1] + 1, ("S", i + 1)
    out, i = [], 0
    while i < n:
        t, i = choice[i]
        out.append(t)
    return "".join(out)


# ---------------------------------------------------------------- the relation of proof.md section 2 (= TwinSpec.v)
class Fail(Exception):
    pass


def inc(w, cnt):
    """Inc(w, w') as a function (I1-I4, IS). Returns w'."""
    w = list(w)
    i = 0
    while i < len(w):
        t = w[i]
        if t == "L":
            i += 1                                   # I3 (skipped L's are counted in bulk)
        elif t == "0":
            cnt["I2"] += 1; w[i] = "1"; return "".join(w)
        elif t == "1":
            cnt["I1"] += 1; w[i] = "0"; i += 1
        else:                                        # IS (NEW): Inc(S::r, r') <- R(h1,1,r,r')
            cnt["IS"] += 1
            return "".join(w[:i]) + R("h1", 1, "".join(w[i + 1:]), cnt)
    cnt["I4"] += 1
    return "".join(w) + "1"


def R(tp, k, w, cnt, parent=False):
    """R(tp,k,w,w') as a partial function (raises Fail where no clause applies).
    parent=True gives the parent's rules: X1/C1 add one pass (k+1) and there is no stray (Z0/Z1 only)."""
    out = []
    i = 0
    while True:
        if tp == "x":
            if i == len(w):
                raise Fail("hx at the right end")
            t = w[i]
            if t == "L":                             # X1: Inc the rest, eater keeps k (parent: k+1)
                j = i
                while j < len(w) and w[j] == "L":
                    j += 1
                rest = w[j:]
                for _ in range(j - i):
                    rest = inc(rest, cnt)
                    k += parent
                cnt["X1"] += j - i
                out.append("L" * (j - i))
                w, i = rest, 0
            elif t == "0":
                cnt["X2"] += 1; out.append("L"); k += 1; i += 1
            elif t == "1":
                cnt["X3"] += 1; out.append("L"); tp = "1"; i += 1
            else:
                raise Fail("hx on S")
        else:
            if k == 0:
                if i == len(w):
                    cnt["R0"] += 1; out.append("1")
                elif parent:
                    if w[i:] in ("1", "0"):
                        cnt["Z1" if w[i] == "1" else "Z0"] += 1
                        out.append("L" if w[i] == "1" else "1")
                    else:
                        raise Fail("parent: h1 with 0 passes")
                else:
                    cnt["Z"] += 1; out.append("S" + w[i:])   # Z (NEW): R(h1,0,r,S::r)
                return "".join(out)
            if i == len(w):                           # R1/R2
                if k % 2:
                    cnt["R1"] += 1; out.append("0"); k = (k - 1) // 2
                else:
                    cnt["R2"] += 1; out.append("1"); k = (k - 2) // 2
                continue
            t = w[i]
            i += 1
            if t == "L":                              # C1: Inc the rest, eater with k-1 (parent: k)
                cnt["C1"] += 1; out.append("L")
                w, i, tp = inc(w[i:], cnt), 0, "x"
                k = k if parent else k - 1
            elif t == "1":
                cnt["C2"] += 1; out.append("L"); k -= 1
            elif t == "0":
                if k % 2:
                    cnt["C3"] += 1; out.append("0"); k = (k - 1) // 2
                else:
                    cnt["C4"] += 1; out.append("1"); k = (k - 2) // 2
            else:
                raise Fail("h1 on S (k >= 1)")


def sweep(w, cnt=None, parent=False):
    return R("x", 0, w, Counter() if cnt is None else cnt, parent=parent)


def csweep(w, cnt=None, parent=False):
    """sweep followed by the canonical re-parse (a tape identity)."""
    return parse(cells(sweep(w, cnt, parent)))


# ---------------------------------------------------------------- real machine
def table(tm):
    T = {}
    for qi, part in enumerate(tm.split("_")):
        for s in range(2):
            t = part[3 * s:3 * s + 3]
            T[("ABCDEF"[qi], s)] = None if t[2] == "-" else (int(t[0]), 1 if t[1] == "R" else -1, t[2])
    return T


TAB = {m: table(m) for m in (TM, PARENT, TM7)}


class Run:
    def __init__(self, tm):
        self.T, self.steps, self.a1, self.fread = TAB[tm], 0, 0, Counter()

    def go(self, tape, pos, q, stop, cap=5 * 10 ** 7, until_state=None):
        for _ in range(cap):
            if pos == stop and (until_state is None or q == until_state):
                return q
            s = tape.get(pos, 0)
            if q == "A" and s == 1:
                self.a1 += 1
            if q == "F":
                self.fread[s] += 1
            tr = self.T[(q, s)]
            if tr is None:
                raise RuntimeError("HALT")
            w, mv, q = tr
            tape[pos] = w
            pos += mv
            self.steps += 1
        raise RuntimeError("TIMEOUT")


def put(tape, o, s):
    for j, ch in enumerate(s):
        tape[o + j] = int(ch)


def read(tape, o):
    hi = max([c for c, v in tape.items() if v] + [o - 1])
    return norm("".join(str(tape.get(c, 0)) for c in range(o, hi + 1)))


def passes(tm, tp, k, w):
    """The pass list of R(tp,k,w) on the real machine; returns (cells read back, Run)."""
    tape, run, o = {}, Run(tm), 0
    if tp == "x":
        put(tape, 0, cells(w))
        tape[-1] = 1
        if run.go(tape, 0, "D", -2) != "C":
            raise RuntimeError("h' exit state")
        o = -1
    else:
        tape[0] = 1
        put(tape, 1, cells(w))
    for _ in range(k):
        if run.go(tape, o, "B", o - 1) != "C":
            raise RuntimeError("h exit state")
    return read(tape, o), run


def sweep_tm(tm, w):
    """Real machine from S(w) until the head first reaches cell -3, then 2 steps; S(w') one cell left?"""
    tape, run = {-2: 1, -1: 1}, Run(tm)
    put(tape, 0, cells(w))
    q = run.go(tape, 0, "D", -3)
    pos = -3
    for _ in range(2):
        wr, mv, q = run.T[(q, tape.get(pos, 0))]
        tape[pos] = wr
        pos += mv
        run.steps += 1
    ok = q == "D" and pos == -1 and tape.get(-3) == 1 and tape.get(-2) == 1
    return (read(tape, -1) if ok else None), run


def seg_table(tm):
    """Per-token pass transduction, assuming every outgoing pass comes back as `<C` (proof.md section 2)."""
    T = TAB[tm]
    res = {}
    for name, w, it in [("L", "110", "h"), ("L", "110", "h'"), ("0", "000", "h"), ("0", "000", "h'"),
                        ("1", "100", "h"), ("1", "100", "h'"), ("S0", "1000", "h"), ("S1", "1100", "h"),
                        ("SL", "1110", "h")]:
        reg = [int(c) for c in w]
        pos, st = 0, "B"
        if it == "h'":
            reg, pos, st = [1] + reg, 1, "D"
        outs = []
        for _ in range(1000):
            if pos < 0:
                break
            if pos >= len(reg):
                if st == "B":
                    outs.append("h"); pos = len(reg) - 1; st = "C"
                elif st == "D" and reg[-1] == 1:
                    outs.append("h'"); reg = reg[:-1]; pos = len(reg) - 1; st = "C"
                else:
                    outs.append("?" + st); break
                continue
            wr, mv, st = T[(st, reg[pos])]
            reg[pos] = wr
            pos += mv
        res[(name, it)] = (st, "".join(map(str, reg)), ",".join(outs))
    return res


# ---------------------------------------------------------------- words
def rword(rng, strays=True):
    kind = rng.random()
    if kind < 0.35:
        w = "".join(rng.choice("L01") for _ in range(rng.randint(0, 14)))
    elif kind < 0.7:
        w = "L" * rng.randint(0, 40) + "".join(rng.choice("L01") for _ in range(rng.randint(0, 10)))
    else:
        w = "".join(rng.choice(["LLL01", "L01", "01", "0001", "LLLLLL1", "11L", "L"]) for _ in range(rng.randint(0, 6)))
    if strays and rng.random() < 0.6:
        j = rng.randint(0, len(w))
        w = w[:j] + "S" + w[j:]
    return w + rng.choice(["1", "01", "L1", "11"])


def main():
    rng = random.Random(20261010)

    # 1. per-token transduction table: twin vs parent vs TM7
    tabs = {m: seg_table(m) for m in (TM, PARENT, TM7)}
    exp = {("L", "h"): "h", ("L", "h'"): "h,h'", ("0", "h"): "", ("0", "h'"): "h',h", ("1", "h"): "h",
           ("1", "h'"): "", ("S0", "h"): "", ("S1", "h"): "", ("SL", "h"): "h,h'"}
    ok = sum(tabs[TM][key][2] == v and tabs[TM][key][0] == "C" for key, v in exp.items())
    report("token table (twin): L h'->[h,h'], SL h->[h,h'], all others as parent", ok, len(exp))
    diff = [key for key in exp if tabs[TM][key] != tabs[PARENT][key]]
    print("   twin vs parent differ only at:", diff,
          {key: (tabs[PARENT][key][2], tabs[TM][key][2]) for key in diff},
          " TM7 at L h':", tabs[TM7][("L", "h'")][2], flush=True)
    report("twin differs from parent exactly at (L,h') and (SL,h)", int(diff == [("L", "h'"), ("SL", "h")]), 1)

    # 2. machine-level L entry: eater and drain, every X <= 5 digits + random, twin and parent
    xs = ["".join(p) for n in range(6) for p in itertools.product("L01", repeat=n)]
    xs += ["".join(rng.choice("L01") for _ in range(rng.randint(6, 40))) for _ in range(150)]
    for tm, pre, ce, cd in [(TM, "0111", 9, 8), (PARENT, "1111", 13, 12)]:
        for ctx in ("eater", "drain"):
            good = 0
            for X in xs:
                # j = number of digits of X before its first 0 digit (or |X|)
                j = next((idx for idx, d in enumerate(X) if d == "0"), len(X))
                xp = inc(X, Counter())
                tape, run = {}, Run(tm)
                if ctx == "eater":                      # 1 D> 110 <X>, the 1 at cell -1
                    tape[-1] = 1; put(tape, 0, "110" + cells(X)); start, q = 0, "D"
                else:                                   # B> 1 110 <X>, B on cell -1
                    put(tape, -1, "1110" + cells(X)); start, q = -1, "B"
                try:
                    run.go(tape, start, q, 3, until_state="D")
                except RuntimeError:
                    continue
                got = "".join(str(tape.get(c, 0)) for c in range(-1, 3))
                steps_ok = run.steps == 6 * j + (ce if ctx == "eater" else cd)
                good += got == pre and read(tape, 3) == norm(cells(xp)) and steps_ok
            report(f"L entry ({'twin' if tm == TM else 'parent'}, {ctx}): -> {pre} D> <X+1> in 6j+"
                   f"{ce if ctx == 'eater' else cd} steps", good, len(xs))

    # 3. every clause on the real machine (section 3.1): pass list of R(tp,k,w) turns the input into <w'>
    cases = []
    for n in range(4 if FULL else 3):
        for p in itertools.product("L01S", repeat=n):
            w = "".join(p)
            for tp in ("x", "1"):
                for k in range(7):
                    cases.append((tp, k, w))
    for _ in range(3000 if FULL else 600):
        cases.append((rng.choice("x1"), rng.randint(0, 60), rword(rng)))
    good = tot = 0
    use = Counter()
    fread1 = 0
    for tp, k, w in cases:
        cnt = Counter()
        try:
            wp = R("x" if tp == "x" else "h1", k, w, cnt)
        except (Fail, RecursionError):
            continue
        tot += 1
        try:
            got, run = passes(TM, "x" if tp == "x" else "1", k, w)
        except RuntimeError:
            continue
        fread1 += run.fread[1]
        if got == norm(cells(wp)):
            good += 1
            use += cnt
    report("clauses on the real machine (all R derivations, exhaustive |w|<=%d k<=6 + random)" % (3 if FULL else 2),
           good, tot, "uses " + ", ".join(f"{a} {b}" for a, b in sorted(use.items())))
    report("F never reads 1 in those runs", int(fread1 == 0), 1, f"({fread1})")

    # 4. BigStep on the real machine for random words in the domain
    good = tot = 0
    zst = 0
    for _ in range(1500 if FULL else 300):
        w = rword(rng)
        try:
            cnt = Counter()
            wp = sweep(w, cnt)
        except (Fail, RecursionError):
            continue
        tot += 1
        got, run = sweep_tm(TM, w)
        good += got == norm(cells(wp))
        zst += cnt["Z"] > 0
    report("BigStep: S(w) -+-> S(sweep w) on the real machine, one cell left (random w)", good, tot,
           f"({zst} leave a stray)")

    # 5. from blank: boundaries follow csweep
    N = 10 ** 9 if FULL else 10 ** 8
    m = sim.parse(TM)
    bounds = []
    for ev in sim.events(m, max_steps=N, on_record=True, snapshot=True):
        if ev.side != "L":
            continue
        c = ev.config
        s = norm("".join(str(x) for x in c.tape[c.head:]))
        bounds.append((ev.step + 2, s[2:]))
    good = 0
    strays = Counter()
    use = Counter()
    par_ok = 0
    for (st, x), (_, y) in zip(bounds, bounds[1:]):
        w = parse(x)
        strays[w.count("S")] += 1
        cnt = Counter()
        try:
            good += norm(cells(sweep(w, cnt))) == y
            use += cnt
        except Fail:
            pass
        try:
            par_ok += norm(cells(sweep(w, parent=True))) == y
        except Fail:
            pass
    report(f"from blank to {N:.0e}: consecutive sweep boundaries follow sweep", good, len(bounds) - 1,
           f"(strays per boundary {dict(strays)}; last {bounds[-1][0]:,} = S({re.sub(r'L{5,}', lambda q: 'L^%d ' % len(q.group(0)), parse(bounds[-1][1]))}))")
    print("   clause uses:", ", ".join(f"{a} {b}" for a, b in sorted(use.items())), flush=True)
    print(f"   the parent's rules reproduce {par_ok}/{len(bounds) - 1} of these sweeps", flush=True)
    report("initial segment c0 -(5)-> S(1)", int(bounds[0] == (5, "1")), 1)

    # 6. abstract orbit
    w, S = "1", 3000 if FULL else 500
    maxs = 0
    okc = 0
    withs = right_after = maxu = 0
    ratio = None
    for s in range(S):
        try:
            w = csweep(w)
            okc += 1
        except (Fail, RecursionError):
            break
        maxs = max(maxs, w.count("S"))
        if "S" in w:
            nn = len(re.match("L*", w).group(0))
            u = w[nn:w.index("S")]
            withs += 1
            right_after += u == ""
            maxu = max(maxu, len(u))
            q = nn / 2 ** len(u)
            ratio = q if ratio is None else min(ratio, q)
    n = len(re.match("L*", w).group(0))
    report(f"abstract orbit from S(1): {S} sweeps stay in the domain", okc, S,
           f"(max strays {maxs}; ends at S(L^{n} {w[n:]}))")
    print(f"   boundaries L^n u S v: {withs} with a stray; u empty in {right_after} ({100 * right_after / max(withs, 1):.0f}%); "
          f"max |u| {maxu}; min n / 2^|u| = {ratio}", flush=True)

    # 7. the naive domain is NOT closed (the open part, proof.md section 3.3)
    L = 7 if FULL else 6
    dom = fail1 = 0
    surv = Counter()
    for n in range(1, L + 1):
        for p in itertools.product("L01S", repeat=n):
            w = "".join(p)
            if w.count("S") > 1 or parse(cells(w)) != w:
                continue
            chain = 0
            cur = w
            for _ in range(9):
                try:
                    cur = csweep(cur)
                    chain += 1
                except (Fail, RecursionError):
                    break
            if chain >= 1:
                dom += 1
                fail1 += chain == 1
            surv[chain] += 1
    print(f"   naive domain (|w|<={L}, <=1 stray): {dom} words have a sweep; {fail1} of them sweep to a word that has "
          f"none. Sweeps survived (cap 9): {dict(sorted(surv.items()))}", flush=True)
    report("evidence for a 3-sweep invariant: no word that survives 3 sweeps fails before 9",
           int(all(v == 0 for c, v in surv.items() if 3 <= c < 9)), 1)

    reduction_checks(rng, bounds)
    review3_checks(rng)

    if FULL:
        rocq_check()
    print("ALL PASS" if not FAILS else "FAILURES: " + "; ".join(FAILS))


# ---------------------------------------------------------------- 8. the reduction of proof.md section 3.3
def canon(w):
    return parse(cells(w))


def free(w):
    return "S" not in w


def lead(w):
    return len(w) - len(w.lstrip("L"))


def E(M, x):
    """E_M(x) := R(hx,M,x) as a word (raises Fail)."""
    return R("x", M, x, Counter())


def theta(x):
    """theta(x): least M with E_M(x) S-free after the canonical re-parse (x S-free and not in 0*).
    Binary search; it relies on monotonicity (Lemma D, checked exhaustively on small words in check 8d)."""
    hi = 1
    while not free(canon(E(hi, x))):
        hi *= 2
    lo = 0
    while lo < hi:
        m = (lo + hi) // 2
        if free(canon(E(m, x))):
            hi = m
        else:
            lo = m + 1
    return lo


def xof(w):
    """Boundary word w = L^n t  ->  (n, x) with x = the canonical form of Inc^n(t)."""
    n = lead(w)
    x = w[n:]
    for _ in range(n):
        x = inc(x, Counter())
    return n, canon(x)


def good(w):
    """G(w) (proof.md 3.3.3): x = Inc^n(t) is S-free and contains a 1 or an L."""
    n, x = xof(w)
    return free(x) and any(c in x for c in "1L")


def words(m, alpha, max_s=1):
    for ln in range(m + 1):
        for p in itertools.product(alpha, repeat=ln):
            w = "".join(p)
            if w.count("S") <= max_s:
                yield w


def reduction_checks(rng, bounds):
    # 8a. Lemma A: R(tp,k,r,r') and Inc(r',r'') give R(tp,k+1,r,r''), same tape
    tot = ok = 0
    for w in list(words(5 if FULL else 4, "L01S")) + [rword(rng) for _ in range(300 if FULL else 60)]:
        for tp in ("x", "h1"):
            for k in range(9):
                try:
                    r2 = inc(R(tp, k, w, Counter()), Counter())
                except (Fail, RecursionError):
                    continue
                tot += 1
                try:
                    ok += norm(cells(R(tp, k + 1, w, Counter()))) == norm(cells(r2))
                except (Fail, RecursionError):
                    pass
    report("8a Lemma A: R(tp,k,r,r') + Inc(r',r'') -> R(tp,k+1,r,r''), same tape (all r <= %d tokens, <= 1 S, k <= 8, "
           "+ random)" % (5 if FULL else 4), ok, tot)

    # 8b. Lemma B: S-free words are safe
    tot = ok = 0
    for w in words(8 if FULL else 6, "L01", 0):
        live = any(c in w for c in "1L")
        for k in range(11):
            tot += 1
            try:
                d = R("h1", k, w, Counter())
                good_d = d.count("S") <= 1
            except (Fail, RecursionError):
                good_d = False
            try:
                e = R("x", k, w, Counter())
                good_e = live and e.count("S") <= 1
            except (Fail, RecursionError):
                good_e = not live
            ok += good_d and good_e
    for w in words(6 if FULL else 5, "L01S"):
        tot += 1
        try:
            ok += inc(w, Counter()).count("S") <= 1
        except (Fail, RecursionError):
            pass
    report("8b Lemma B: on S-free r, R(h1,k,r) always exists, R(hx,k,r) exists iff r has a 1 or an L, each output has "
           "<= 1 S; Inc of a word with <= 1 S exists with <= 1 S (exhaustive)", ok, tot)

    # 8c/8d. Lemma C (G(w) -> sweep defined, = real machine) and Lemma D (G(sweep w) <-> E_{n+a}(x) S-free)
    ts = [t for t in words(5 if FULL else 4, "L01S") if not t.startswith("L") and canon(t) == t and t]
    ns = [0, 1, 2, 3, 5, 8, 13] if FULL else [0, 1, 3, 8]
    totc = okc = totd = okd = gtot = gclosed = 0
    example = None
    for t in ts:
        for n in ns:
            w = "L" * n + t
            if not good(w):
                continue
            gtot += 1
            totc += 1
            _, x = xof(w)
            try:
                out = "L" * n + E(0, x)
            except (Fail, RecursionError):
                continue
            got, _ = sweep_tm(TM, w)
            okc += out.count("S") <= 1 and got == norm(cells(out))
            w2 = canon(out)
            a = lead(E(0, x))
            totd += 1
            okd += good(w2) == free(canon(E(n + a, x)))
            gclosed += good(w2)
            if example is None and not good(w2):
                example = (w, w2)
    report("8c Lemma C: G(w) -> L^n E_0(x) exists, has <= 1 S, and equals the real machine's sweep (t <= %d tokens, "
           "n in %s)" % (5 if FULL else 4, ns), okc, totc)
    report("8d Lemma D: G(sweep w) <-> E_{n+a}(x) is S-free (same words)", okd, totd)
    # monotonicity of the threshold: E_M(x) S-free -> E_{M+1}(x) S-free
    tot = ok = 0
    for x in words(7 if FULL else 6, "L01", 0):
        if not any(c in x for c in "1L"):
            continue
        prev = False
        for M in range(48):
            cur = free(canon(E(M, x)))
            tot += 1
            ok += not (prev and not cur)
            prev = cur
    report("8d monotone threshold: E_M(x) S-free -> E_{M+1}(x) S-free (S-free x <= %d tokens, M < 48)"
           % (7 if FULL else 6), ok, tot)
    print(f"   G is not closed: of {gtot} words with G, {gtot - gclosed} sweep to a word without G; first: "
          f"{example[0]} -> {example[1]}", flush=True)
    report("8f G not closed: 1001 has G, its sweep LS001 does not (Inc(S001) = 0S01)",
           int(good("1001") and canon(sweep("1001")) == "LS001" and not good("LS001")
               and xof("LS001")[1] == "0S01"), 1)

    # 8e. G along the abstract orbit and on the from-blank boundaries; the quantitative picture
    S = 3000 if FULL else 500
    w = "1"
    okg = end1 = cons = 0
    worst = (0, None)
    vmax = 0
    prev = None
    for s in range(S):
        n, x = xof(w)
        g = free(x) and any(c in x for c in "1L")     # = good(w), without recomputing x
        okg += g
        if not g:
            break
        end1 += x.endswith("1")
        e0 = E(0, x)
        a = lead(e0)
        th = theta(x)
        assert free(canon(E(th, x))) and (th == 0 or not free(canon(E(th - 1, x))))
        if 2 * th / (n + a) > worst[0]:
            worst = (2 * th / (n + a), (s, n + a, th))
        v = sum(2 ** i for i, c in enumerate(y for y in x if y in "01") if c == "1")
        vmax = max(vmax, (v + lead(x)) / (n + a))
        if prev is not None:
            cons += v == prev[0] - prev[1]
        prev = (n + a, th)
        w = canon("L" * n + e0)
    report(f"8e G holds at every boundary of the abstract orbit ({S} sweeps)", okg, S)
    report(f"8e along the orbit, x_s ends with the token 1", end1, okg)
    report(f"8e along the orbit, 2 theta(x_s) <= n_(s+1) (max ratio {worst[0]:.3f} at sweep, n_(s+1), theta = "
           f"{worst[1]})", int(worst[0] <= 1), 1)
    report(f"8e along the orbit, V(x_s) + leading L's of x_s <= 1.25 n_(s+1) (max {vmax:.3f})", int(vmax <= 1.25), 1)
    print(f"   V(x_(s+1)) = N_s - theta(x_s) exactly on {cons} of {okg - 1} sweeps", flush=True)
    okb = sum(good(parse(x)) for _, x in bounds)
    report("8e G holds at every from-blank boundary", okb, len(bounds))

    # 8f. the bare candidate J0 := theta(x) <= N is not inductive
    x, N = "00010001", 5
    a = lead(E(0, x))
    xp = canon(E(N, x)[a:])
    a2 = lead(E(0, xp))
    report("8f J0 not inductive: x = 00010001, N = 5 has theta(x) = 4 <= N, x' = 100L, N' = 6, theta(x') = 7",
           int(theta(x) == 4 and a == 4 and xp == "100L" and a2 == 1 and theta(xp) == 7), 1)
    # 8f. one step of J0 with margin m = N - theta(x) >= 4 keeps theta(x') <= N'; but the margin can drop below m/2
    tot = ok = 0
    for x in words(7 if FULL else 6, "L01", 0):
        if not any(c in x for c in "1L"):
            continue
        th = theta(x)
        a = lead(E(0, x))
        for m in range(4, 40):
            xp = canon(E(th + m, x)[a:])
            tot += 1
            ok += theta(xp) <= th + m + lead(E(0, xp))
    report("8f J0 one step with margin: N - theta(x) >= 4 gives theta(x') <= N' (S-free x <= %d tokens, margin < 40)"
           % (7 if FULL else 6), ok, tot)
    x, N = "L11111L", 26
    a = lead(E(0, x))
    xp = canon(E(N, x)[a:])
    mp = N + lead(E(0, xp)) - theta(xp)
    report(f"8f margin can drop below half: x = L11111L, N = 26 (theta 0, margin 26) gives x' = {xp}, margin {mp}",
           int(theta(x) == 0 and 2 * mp < N), 1)
    # 8f. the strengthened candidate J2 := (2 theta(x) <= N and x ends with 1) is not inductive either
    x, N = "L1111", 60
    a = lead(E(0, x))
    xp = canon("".join(E(N, x))[a:])
    a2 = lead(E(0, xp))
    th, th2 = theta(x), theta(xp)
    report("8f J2 not inductive: x = L1111, N = 60 has 2 theta(x) = 0 <= N, x' = 1000001, N' = 61, theta(x') = 31",
           int(th == 0 and xp == "1000001" and a2 == 1 and th2 == 31 and 2 * th2 > N + a2), 1)

    # 8f. J3 := J2 and V(x) + leading L's <= 5N/4 is not inductive either (at small N)
    def val(y):
        return sum(2 ** i for i, c in enumerate(z for z in y if z in "01") if c == "1")

    def j3(n_, y):
        return y.endswith("1") and 2 * theta(y) <= n_ and 4 * (val(y) + lead(y)) <= 5 * n_

    x, N = "L11", 4
    a = lead(E(0, x))
    xp = canon(E(N, x)[a:])
    N2 = N + lead(E(0, xp))
    report("8f J3 not inductive: x = L11, N = 4 satisfies J3; x' = 111, N' = 5, V(x') = 7 > 25/4",
           int(j3(N, x) and xp == "111" and N2 == 5 and val(xp) == 7 and not j3(N2, xp)), 1)


# ---------------------------------------------------------------- 9. pass 3 (review of 2026-10-10): gaps and the V* reformulation
def vstar(y):
    """V*(y) (proof.md 3.3.5): little-endian binary value where an L standing after p digit tokens counts 2^p."""
    p = v = 0
    for ch in y:
        if ch == "L":
            v += 2 ** p
        else:
            v += 2 ** p if ch == "1" else 0
            p += 1
    return v


def val(y):
    return sum(2 ** i for i, c in enumerate(z for z in y if z in "01") if c == "1")


def theta_brute(x, cap):
    """theta(x) by linear search, without using monotonicity; None if above cap."""
    for M in range(cap + 1):
        if free(canon(E(M, x))):
            return M
    return None


def system_step(N, x):
    """(N, x) -> (N', x') of proof.md 3.3.4 (requires theta(x) <= N)."""
    a = lead(E(0, x))
    xp = canon(E(N, x))[a:]
    return N + lead(E(0, xp)), xp


def cfun(x):
    """c(x) := V*(y(x)) - theta(x), y(x) the threshold output without its first a(x) L's."""
    th = theta(x)
    return vstar(canon(E(th, x))[lead(E(0, x)):]) - th


def review3_checks(rng):
    # 9a. Lemma C needs the re-parse INSIDE the derivation: the 16 clauses alone fail on some G words
    T, NM = (6, 15) if FULL else (5, 8)
    ts = [t for t in words(T, "L01S") if t and not t.startswith("L") and canon(t) == t]
    gw = strict_fail = eq_ok = 0
    fails = []
    for t in ts:
        for n in range(NM + 1):
            w = "L" * n + t
            if not good(w):
                continue
            gw += 1
            try:
                R("x", 0, w, Counter())            # TwinSpec's 16 clauses (R0 chosen over Z at []; Z there only loses)
            except (Fail, RecursionError):
                strict_fail += 1
                fails.append(w)
            _, x = xof(w)
            out = "L" * n + E(0, x)                # RIncs_eq once, after the eater's L-run
            eq_ok += out.count("S") <= 1
    report(f"9a Lemma C with RIncs_eq: every G word (tails <= {T} tokens, n <= {NM}) sweeps, <= 1 S", eq_ok, gw,
           f"(without RIncs_eq the clauses fail on {strict_fail} of {gw}, first {fails[0] if fails else None})")
    real = sum(sweep_tm(TM, w)[0] == norm(cells("L" * lead(w) + E(0, xof(w)[1])))
               for w in fails[:80])
    report("9a those strict failures: the RIncs_eq sweep equals the real machine's", real, min(len(fails), 80))
    report("9a LS01: strict derivation fails, RIncs_eq gives LLLL01 (= scripts/TwinSpec.v Example sweep_LS01)",
           int("LS01" in fails and "L" + E(0, xof("LS01")[1]) == "LLLL01"
               and sweep_tm(TM, "LS01")[0] == norm(cells("LLLL01"))), 1)

    # 9b. Lemma D, B2 precondition: an S-free x(w') always contains a 1 or an L
    tot = ok = 0
    for t in ts:
        for n in range(NM + 1):
            w = "L" * n + t
            if not good(w):
                continue
            w2 = canon("L" * n + E(0, xof(w)[1]))
            x2 = xof(w2)[1]
            if free(x2):
                tot += 1
                ok += any(c in x2 for c in "1L") and lead(w2) >= 1
    report("9b Lemma D: when x(w') is S-free it contains a 1 or an L, and n' >= 1", ok, tot)

    # 9c. J3 with N >= N0 is refuted for every N0: (N, x) = (2^(m+1) - 4, 00001)
    ms = range(4, 21) if FULL else range(4, 13)
    ok = 0
    for m in ms:
        N, x = 2 ** (m + 1) - 4, "00001"
        th = theta_brute(x, 0)
        j3 = th == 0 and x.endswith("1") and 4 * (val(x) + lead(x)) <= 5 * N and 2 * th <= N
        N2, x2 = system_step(N, x)
        th2 = theta_brute(x2, 2 ** m) if m <= 10 else theta(x2)
        if m > 10:   # monotone (Lemma A + I1-I4): check both sides of the threshold explicitly
            th2 = th2 if free(canon(E(th2, x2))) and not free(canon(E(th2 - 1, x2))) else None
        ok += j3 and x2 == "1" + "0" * m + "1" and N2 == 2 ** (m + 1) - 3 and th2 == 2 ** m - 1 and 2 * th2 == N2 + 1
    report(f"9c J3 + (N >= N0) refuted: (2^(m+1)-4, 00001) is in J3, x' = 1 0^m 1, N' = 2^(m+1)-3, 2 theta' = N'+1 "
           f"(m = {ms.start}..{ms.stop - 1}; brute-force theta for m <= 10)", ok, len(ms))

    # 9d. the orbit's value bound V + lead <= 5N/4 is attained exactly at sweep 2 (N = 4, x = 101)
    w = "1"
    best = (0, None)
    for s in range(3000 if FULL else 500):
        n, x = xof(w)
        e0 = E(0, x)
        N = n + lead(e0)
        r = (val(x) + lead(x)) / N
        if r > best[0]:
            best = (r, (s, N, x))
        w = canon("L" * n + e0)
    report(f"9d orbit: max (V + lead)/N = {best[0]} at (sweep, N, x) = {best[1]}: the 5/4 constant has no slack",
           int(best[0] == 1.25 and best[1] == (2, 4, "101")), 1)

    # 9e. V*: Inc is +1 on V* for S-free words
    tot = ok = 0
    for x in words(8 if FULL else 6, "L01", 0):
        tot += 1
        ok += vstar(inc(x, Counter())) == vstar(x) + 1
    report("9e Inc adds exactly 1 to V* on S-free words (exhaustive)", ok, tot)

    # 9f. affine law: V*(drop_a(E_M(x))) = M + c(x) for every M >= theta(x)
    tot = ok = 0
    for x in words(7 if FULL else 6, "L01", 0):
        if not x.endswith("1"):
            continue
        th, a, c = theta(x), lead(E(0, x)), cfun(x)
        for M in range(th, th + 12):
            tot += 1
            ok += vstar(canon(E(M, x))[a:]) == M + c
    report("9f affine law: V*(x') = N + c(x) whenever theta(x) <= N (S-free x ending with 1, exhaustive, 12 values of N)",
           ok, tot)

    # 9g. Lemma T (CONJECTURE): 2 theta(x) <= V*(x) - 1 for S-free x ending with the token 1
    L_ = 11 if FULL else 8
    tot = ok = 0
    tight = 0
    for x in words(L_, "L01", 0):
        if not x.endswith("1"):
            continue
        tot += 1
        d = vstar(x) - 1 - 2 * theta(x)
        ok += d >= 0
        tight += d == 2
    xs_ = ["1" + "0" * m + "1" for m in range(1, 9)]
    report(f"9g Lemma T (conjecture): 2 theta(x) <= V*(x) - 1, S-free x ending with 1 (exhaustive <= {L_} tokens)", ok, tot,
           f"({tight} words have 2 theta = V* - 3)")
    report("9g T is tight on 1 0^m 1 (2 theta = V* - 3) and false for words ending with L (100000000L: theta 511, V* 513)",
           int(all(2 * theta(x) == vstar(x) - 3 for x in xs_) and theta("100000000L") == 511
               and vstar("100000000L") == 513), 1)
    tot = ok = 0
    for x in words(L_, "01", 0):
        if x.endswith("1"):
            tot += 1
            ok += cfun(x) <= len(x)
    report(f"9g pure binary x (no L) ending with 1: c(x) <= |x| (exhaustive <= {L_} tokens)", ok, tot)
    report("9g L-driven c: c(10000000L1) = 385 (V* 513), c(L0000LLL1) = 57 (V* 65), c(1000L11111) = 89 (Lambda 16)",
           int(cfun("10000000L1") == 385 and vstar("10000000L1") == 513 and cfun("L0000LLL1") == 57
               and vstar("L0000LLL1") == 65 and cfun("1000L11111") == 89), 1)

    # 9h. along the orbit: V*(x_(s+1)) = N_s + c(x_s) exactly; V*/N and c/N ranges
    w = "1"
    S = 3000 if FULL else 500
    aff = 0
    lo, hi, cmax, t_ok = 9.0, 0.0, 0.0, 0
    prev = None
    for s in range(S):
        n, x = xof(w)
        e0 = E(0, x)
        N = n + lead(e0)
        v = vstar(x)
        if prev is not None:
            aff += v == prev
        c = cfun(x)
        t_ok += 2 * theta(x) <= v - 1 or x == "1"
        if s >= 27:
            lo, hi, cmax = min(lo, v / N), max(hi, v / N), max(cmax, c / N)
        prev = N + c
        w = canon("L" * n + e0)
    report(f"9h orbit ({S} sweeps): V*(x_(s+1)) = N_s + c(x_s)", aff, S - 1)
    report(f"9h orbit: Lemma T's inequality at every x_s", t_ok, S)
    print(f"   orbit, s >= 27: {lo:.3f} <= V*(x_s)/N_s <= {hi:.3f}; c(x_s) <= {cmax:.3f} N_s", flush=True)

    # 9i. J4 := x ends with 1 and V*(x) <= 2N + 1 (which gives theta <= N under T) is NOT inductive ...
    N, x = 257, "10000000L1"
    N2, x2 = system_step(N, x)
    report("9i J4 not inductive: (257, 10000000L1) has V* = 513 <= 515, theta = 255; x' = 0100000LLL01, N' = 259, "
           "V*(x') = 642 > 519, but theta(x') = 94 <= N'",
           int(vstar(x) == 513 and theta(x) == 255 and (N2, x2) == (259, "0100000LLL01") and vstar(x2) == 642
               and theta(x2) == 94), 1)

    # 9j. ... but iterating the system from random J4 states never leaves the domain, and V*/N contracts
    starts = 400 if FULL else 40
    K = 25 if FULL else 12
    fails_ = tot = 0
    after2 = 0.0

    def rx():
        k = rng.random()
        if k < 0.3:
            y = "".join(rng.choice("L01") for _ in range(rng.randint(1, 18)))
        elif k < 0.6:
            y = "".join(rng.choice("0001LL") for _ in range(rng.randint(1, 18)))
        elif k < 0.8:
            y = "1" + "0" * rng.randint(3, 12) + "L" * rng.randint(1, 4) + rng.choice(["1", "11", "01"])
        else:
            y = "".join(rng.choice(["L", "LLLL", "0", "000", "1", "11"]) for _ in range(rng.randint(1, 10)))
        return y.rstrip("0").rstrip("L") + "1"

    while tot < starts:
        x = rx()
        v = vstar(x)
        lo_ = max(131, (v + 1) // 2)
        N = rng.randint(lo_, lo_ + max(3, v // 20)) if rng.random() < 0.6 else rng.randint(lo_, max(lo_, 4 * v))
        tot += 1
        for i in range(K):
            if theta(x) > N:
                fails_ += 1
                break
            if i >= 2:
                after2 = max(after2, vstar(x) / N)
            N, x = system_step(N, x)
    report(f"9j from {starts} random J4 states (N >= 131, biased to V* ~ 2N), {K} sweeps each: theta <= N throughout",
           starts - fails_, starts, f"(max V*/N from the third state on: {after2:.3f})")


ROCQ_SPEC = open(os.path.join(os.path.dirname(os.path.abspath(__file__)), "scripts", "TwinSpec.v")).read() \
    if os.path.exists(os.path.join(os.path.dirname(os.path.abspath(__file__)), "scripts", "TwinSpec.v")) else ""


def rocq_check():
    """Compile scripts/TwinSpec.v in a fresh temp dir; pass iff coqc exits 0, no escape hatch, 5 x closed."""
    banned = re.findall(r"\b(Admitted|admit|Axiom|Parameter|Hypothesis|Conjecture|Abort|bypass_check)\b|"
                        r"Unset Guard Checking", ROCQ_SPEC)
    with tempfile.TemporaryDirectory(prefix="twin2_rspec_") as d:
        with open(os.path.join(d, "TwinSpec.v"), "w") as f:
            f.write(ROCQ_SPEC)
        v = f"{ROOT}/vendor/busycoq/verify"
        cmd = (f"eval $(opam env --switch=bb6) && ulimit -s unlimited && (ulimit -v 8000000 2>/dev/null || true) && cd {d} && "
               f"coqc -Q {v} BusyCoq -Q {v}/BigInt BigInt -native-compiler yes TwinSpec.v")
        p = subprocess.run(["bash", "-c", cmd], capture_output=True, text=True, timeout=1200)
    out = p.stdout + p.stderr
    closed = out.count("Closed under the global context")
    ok = bool(ROCQ_SPEC) and p.returncode == 0 and not banned and closed == 5
    report("Rocq: Inc_spec + RIncs_spec (mutual, 16 clauses + RIncs_eq), end identities, BigStep, init, Example "
           "sweep_LS01 compile; 5 x Print Assumptions closed",
           int(ok), 1, f"(coqc exit {p.returncode}, closed {closed}/5, banned words {len(banned)})")
    if not ok:
        print(out[-3000:])


if __name__ == "__main__":
    if ROCQ_ONLY and not FULL:
        rocq_check()
        print("ALL PASS" if not FAILS else "FAILURES: " + "; ".join(FAILS))
    else:
        main()