#!/usr/bin/env python3
"""check_proof.py -- simulation evidence for proof.md of
1RB0LE_1RC0RC_1LD0RB_1LA0LE_1LD1LF_1LD---   (machine #11).

Self-contained: its own direct simulator, its own abstract (macro) model, and the invariant
INV = Common and (A or B or C) of proof.md section 3. Check 11 proves the 20 window lemmas of
proof.md 2.1 by finite runs; everything else samples.

  python3 check_proof.py [--seed S] [--orbit N] [--induct N] [--tm-steps T] [--quick]

Checks (each prints PASS/FAIL with counts):
  1 sweeps        B> (01)^n -> (10)^n B>  and  (01)^n <D -> <D (10)^n, 2n steps each
  2 rule M_k      random full K-configs (any bits, any u>=0/v>=1, drops included), exact tape,
                  exact step count, and the head never goes left of the 0^2 of the digit it flips
  3 rule O        random all-ones K-configs, exact tape and step count
  4 init          blank tape -> K([(3,1,1)];2,1) at step 142 -> K([(1,1,2)];8,5) at step 483
  5 closure       the real machine from the blank tape agrees with the abstract model at every
                  macro boundary up to --tm-steps (tape, state, head, step number)
  6 orbit         abstract orbit from step 483: INV holds at every macro state (--orbit steps)
  7 induction     random INV states of each class A/B/C (not just reachable ones): one macro step
                  is defined (no carry hits a v=0 digit, no LSB drop) and lands in INV again
  8 induction-TM  a sample of those random INV states, run on the real machine for the predicted
                  number of steps, equals the abstract successor (so merges are right too)
  9 drop algebra  the closed forms (F2/F3 of proof.md) for R, R1, S, T after a drop
 10 twins         the same table started in B / C reaches a class-A state (other two holdouts)
 11 lemmas        the 20 window lemmas S1-S6, T1-T14 of proof.md 2.1, by exact finite runs that stay
                  inside the window (a proof for each lemma; S1-S6 then follow by induction on n)
 12 M chain       the derivation of M in proof.md 2.2, link by link on the real machine: every
                  intermediate configuration, the step formula, the leftmost cell (exhaustive box + random)
 13 O chain       the same for O (proof.md 2.3)
"""
import sys, random, argparse

TM = "1RB0LE_1RC0RC_1LD0RB_1LA0LE_1LD1LF_1LD---"

# ---------------------------------------------------------------- simulator
def parse_tm(s):
    tab = []
    for grp in s.split('_'):
        row = []
        for i in (0, 3):
            t = grp[i:i + 3]
            row.append(None if t == '---' else (int(t[0]), 1 if t[1] == 'R' else -1, ord(t[2]) - 65))
        tab.append(row)
    return tab
TT = parse_tm(TM)

class Halt(Exception):
    pass

def run(cells, head, state, n, pad=64):
    """Run n steps. cells: str of 0/1; head: index into cells (may be len(cells) = first blank);
    state: 'A'..'F'. Returns (tape_str, head_index_in_tape_str, state, min_head_index_relative_to
    cells[0]). Raises Halt on '---'."""
    t = bytearray(pad) + bytearray(int(c) for c in cells) + bytearray(pad + n + 8)
    h = head + pad; st = ord(state) - 65; mn = h; tab = TT
    for _ in range(n):
        tr = tab[st][t[h]]
        if tr is None:
            raise Halt()
        t[h] = tr[0]; h += tr[1]; st = tr[2]
        if h < mn:
            mn = h
            if h < 1:
                return run(cells, head, state, n, pad * 8)
    # tape string and head index are in the same (padded) coordinates; min head is relative to cells[0]
    return ''.join('01'[x] for x in t), h, chr(st + 65), mn - pad

def same_cfg(res, cells, head, state):
    """res from run(); compare with configuration (cells, head, state) up to translation."""
    rs, rh, rst, _ = res
    if rst != state:
        return False
    f = rs.find('1')
    g = cells.find('1')
    if f < 0 or g < 0:
        return False
    if rh - f != head - g:
        return False
    return rs.strip('0') == cells.strip('0')

# ---------------------------------------------------------------- abstract model
# A K-state is (D, a, b); D = list of digits (u, beta, v) from the LEFT (j = 1 .. m), the leading
# "10 01" excluded. Tape: 0^inf 10 01 [(10)^u (00 if beta=0) (01)^v]_j (10)^a 0^3 (10)^b B> 0^inf
def render(D, a, b):
    s = '1001' + ''.join('10' * u + ('' if beta else '00') + '01' * v for (u, beta, v) in D)
    s += '10' * a + '000' + '10' * b
    return s, len(s), 'B'          # head on the first blank cell right of the text, state B

def mstep(D, a, b):
    """One macro step. Returns (name, D2, a2, b2, steps, info) or raises ValueError (no rule)."""
    m = len(D)
    if m >= 1 and all(d[1] == 1 for d in D):                          # rule O
        if any(d[0] < 1 or d[2] < 1 for d in D):
            raise ValueError('O: non-canonical digit')
        u = [d[0] for d in D]; v = [d[2] for d in D]
        D2 = [(0, 0, u[0])] + [(v[i - 1] - 1, 0, u[i]) for i in range(1, m)] + [(v[m - 1], 1, a)]
        steps = 89 - 4 * m + 8 * a + 20 * b + 4 * sum(u) + 8 * sum(v)
        return 'O', D2, b + 3, 1, steps, {}
    if m == 0:
        raise ValueError('empty digit list')
    k = 0
    while D[m - 1 - k][1] == 1:
        k += 1
    js = m - 1 - k                                                     # 0-based index of j*
    steps = 32 + 8 * a + 8 * b
    D2 = list(D)
    for i in range(js + 1, m):                                         # 1 -> 0
        u, _, v = D[i]
        if v < 1 or u < 1:
            raise ValueError('M: bit-1 digit with u or v = 0')
        D2[i] = (u, 0, v - 1); steps += 8 * u + 4 * v
    u, _, v = D[js]
    if v < 1:
        raise ValueError('M%d: carry reaches a bit-0 digit with v = 0' % k)   # the failure mode
    steps += 4 * v
    info = {'k': k, 'jstar': js + 1, 'drop': None}
    if v - 1 >= 1:
        D2[js] = (u + 1, 1, v - 1)
        a2 = a + 1
    elif js == m - 1:                                                  # LSB drop: (10)^(u+1) joins a
        del D2[js]; a2 = a + 1 + u + 1; info['drop'] = 'LSB'
    else:                                                              # drop: joins the next digit
        uu, bb, vv = D2[js + 1]
        D2[js + 1] = (uu + u + 1, bb, vv); del D2[js]; a2 = a + 1
        info['drop'] = 'top' if js == 0 else 'mid'
    return 'M%d' % k, D2, a2, b + 2, steps, info

# ---------------------------------------------------------------- invariant
def stats(D):
    """per digit j=1..m: (u, beta, v, R, R1, ubar, vbar); plus V."""
    out = []; U = 0
    for j, (u, beta, v) in enumerate(D, start=1):
        U = 2 * U + beta
        R = 2 ** j - 1 - U
        R1 = (R + 1 - beta) // 2
        out.append((u, beta, v, R, R1, u + R1, v - R))
    return out, U

def inv(D, a, b):
    """Returns 'A', 'B' or 'C' if INV holds, else None. proof.md section 3."""
    m = len(D)
    if m < 1 or a < 1 or b % 2 == 0:
        return None
    st, V = stats(D)
    Rtot = 2 ** m - 1 - V
    S = a + Rtot - 2 ** (m + 1) + 2
    T = b - 2 * V + 1
    if (a + Rtot) % 2 or S < 0 or T < 0:
        return None
    for (u, beta, v, R, R1, ub, vb) in st:
        if v < 1 or (beta == 1 and u < 1) or u < 0:
            return None
    if st[-1][6] < 1:                                   # LSB keeps spare fuel
        return None
    top = D[0]
    par = all((v + beta) % 2 == 1 for (u, beta, v, R, R1, ub, vb) in st)
    if (top[1] == 0 and top[0] == 0 and top[2] >= 2) or (top[1] == 1 and top[0] == 1):
        if S >= 1 and all(vb >= 1 and ub % 2 == 1 for (u, beta, v, R, R1, ub, vb) in st):
            return 'A'
        return None
    if top == (0, 0, 1):
        if par and all(st[j - 1][5] >= 2 ** (j - 1) for j in range(2, m + 1)):
            return 'B'
        return None
    if par and S >= 1 and T >= 2 and all(st[j - 1][5] >= 2 ** j for j in range(1, m + 1)):
        return 'C'
    return None

# ---------------------------------------------------------------- random INV states
def rnd(r, lo, big):
    x = r.random()
    if x < 0.45: return lo + r.randint(0, 2)
    if x < 0.85: return lo + r.randint(0, 40)
    return lo + r.randint(0, big)

def rand_inv_state(r, cls, m, big=300):
    """A random state satisfying INV in class cls (rejection on the few global constraints)."""
    for _ in range(200):
        bits = [r.randint(0, 1) for _ in range(m)]
        if cls == 'B':
            bits[0] = 0
        if m >= 1 and all(bits) and r.random() < 0.5:
            pass                                          # keep some all-ones (pre-O) states
        st0 = []; U = 0
        for j in range(1, m + 1):
            U = 2 * U + bits[j - 1]; R = 2 ** j - 1 - U; R1 = (R + 1 - bits[j - 1]) // 2
            st0.append((R, R1))
        D = []
        ok = True
        for j in range(1, m + 1):
            beta = bits[j - 1]; R, R1 = st0[j - 1]
            umin = 1 if beta else 0
            if cls == 'A':
                if j == 1:
                    u = 1 if beta else 0
                    v = rnd(r, 1 if beta else 2, big)
                else:
                    ub = R1 + umin + rnd(r, 0, big); ub += (ub % 2 == 0)       # odd ubar
                    u = ub - R1
                    v = R + 1 + rnd(r, 0, big)
            else:
                if cls == 'B' and j == 1:
                    u, v = 0, 1
                else:
                    need = 2 ** (j - 1) if cls == 'B' else 2 ** j
                    ub = max(need, R1 + umin) + rnd(r, 0, big)
                    u = ub - R1
                    v = rnd(r, 1, big)
                    if (v + beta) % 2 == 0: v += 1                              # safe parity
                    if j == m and v - R < 1:
                        v = R + 1 + (R + 1 + beta + 1) % 2                      # LSB fuel, parity
            if u < umin: ok = False
            D.append((u, beta, v))
        if not ok:
            continue
        st, V = stats(D)
        Rtot = 2 ** m - 1 - V
        Tmin = 2 if cls == 'C' else 0
        b = 2 * V - 1 + Tmin + 2 * rnd(r, 0, big // 4)
        if b < 1: b = 1
        if b % 2 == 0: b += 1
        Smin = 0 if cls == 'B' else 1
        a = max(1, 2 ** (m + 1) - 2 + Smin - Rtot) + rnd(r, 0, big)
        if (a + Rtot) % 2: a += 1
        c = inv(D, a, b)
        if c == cls:
            return D, a, b
    return None

# ---------------------------------------------------------------- checks
def report(name, ok, tot, extra=''):
    print('%-14s %s  %d/%d %s' % (name, 'PASS' if ok == tot else 'FAIL', ok, tot, extra))
    return ok == tot

def check_sweeps():
    ok = tot = 0
    for n in list(range(0, 12)) + [50, 333, 1000]:
        for (f, fh, fs, t, th) in [
            ('1' + '01' * n + '1', 1, 'B', '1' + '10' * n + '1', 1 + 2 * n),          # B> (01)^n
            ('11' + '01' * n + '1', 1 + 2 * n, 'D', '11' + '10' * n + '1', 1)]:       # (01)^n <D
            res = run(f, fh, fs, 2 * n)
            tot += 1; ok += same_cfg(res, t, th, fs)
    return report('1 sweeps', ok, tot)

def rand_digit(r, beta, big):
    return (rnd(r, 1 if beta else 0, big), beta, rnd(r, 1, big))

def check_M(r, N, big):
    ok = tot = 0; maxk = 0; drops = 0; leftok = 0
    for it in range(N):
        k = it % 7
        pre = [rand_digit(r, r.randint(0, 1), big) for _ in range(r.randint(0, 3))]
        flip = rand_digit(r, 0, big)
        if r.random() < 0.25: flip = (0, 0, flip[2])        # u = 0 (fresh digit after O)
        if r.random() < 0.25: flip = (flip[0], 0, 1)        # v = 1: the flipped digit drops
        ones = [rand_digit(r, 1, big) for _ in range(k)]
        D = pre + [flip] + ones
        a = rnd(r, 1, big); b = rnd(r, 0, big)
        name, D2, a2, b2, steps, info = mstep(D, a, b)
        s, h, q = render(D, a, b)
        res = run(s, h, q, steps)
        t, th, tq = render(D2, a2, b2)
        good = same_cfg(res, t, th, tq)
        # leftmost cell visited must be inside the flipped digit's 0^2 (or to its right)
        start00 = 4 + sum(2 * u + (0 if be else 2) + 2 * v for (u, be, v) in pre) + 2 * flip[0]
        lgood = res[3] >= start00
        tot += 1; ok += good and lgood; leftok += lgood
        drops += info['drop'] is not None; maxk = max(maxk, k)
        if not (good and lgood):
            print('  M FAIL', D, a, b, good, lgood)
    return report('2 rule M_k', ok, tot, '(k=0..%d, %d with a drop/merge, head never left of the flipped 0^2: %d/%d)'
                  % (maxk, drops, leftok, tot))

def check_O(r, N, big):
    ok = tot = 0
    for it in range(N):
        n = 1 + it % 7
        D = [rand_digit(r, 1, big) for _ in range(n)]
        a = rnd(r, 1, big); b = rnd(r, 0, big)
        name, D2, a2, b2, steps, _ = mstep(D, a, b)
        s, h, q = render(D, a, b)
        res = run(s, h, q, steps)
        t, th, tq = render(D2, a2, b2)
        good = same_cfg(res, t, th, tq)
        tot += 1; ok += good
        if not good: print('  O FAIL', D, a, b)
    return report('3 rule O', ok, tot, '(n=1..7)')

def check_init():
    ok = 0
    res = run('', 0, 'A', 142)
    t, th, tq = render([(3, 1, 1)], 2, 1); ok += same_cfg(res, t, th, tq)
    res = run('', 0, 'A', 483)
    t, th, tq = render([(1, 1, 2)], 8, 5); ok += same_cfg(res, t, th, tq)
    D, a, b, tt = [(3, 1, 1)], 2, 1, 142
    for _ in range(3):
        _, D, a, b, st, _ = mstep(D, a, b); tt += st
    ok += (tt == 483 and (D, a, b) == ([(1, 1, 2)], 8, 5))
    ok += inv(D, a, b) == 'A'
    return report('4 init', ok, 4, '(blank -142-> K0, -483-> K*=K([(1,1,2)];8,5) = 3 macro steps from K0, K* in class A)')

def check_closure(Tmax):
    """real TM vs abstract model from the blank tape."""
    s, h, q = '', 0, 'A'
    res = run(s, h, q, 142); cur = res
    D, a, b, t = [(3, 1, 1)], 2, 1, 142
    n = ok = 0
    while True:
        _, D2, a2, b2, st, _ = mstep(D, a, b)
        if t + st > Tmax:
            break
        cs, ch, cq, _ = cur
        f = cs.find('1'); l = cs.rfind('1')
        cur = run(cs[f:l + 1], ch - f, cq, st)
        D, a, b, t = D2, a2, b2, t + st
        tt, th, tq = render(D, a, b)
        n += 1; ok += same_cfg(cur, tt, th, tq)
    return report('5 closure', ok, n, '(every macro step of the real run to step %d; last macro state m=%d)' % (t, len(D)))

def check_orbit(N):
    D, a, b, t = [(1, 1, 2)], 8, 5, 483
    from collections import Counter
    cls = Counter(); drops = Counter(); ok = 0
    for i in range(N):
        c = inv(D, a, b)
        cls[c] += 1; ok += c is not None
        if c is None:
            print('  orbit FAIL at macro step', i, D, a, b); break
        name, D, a, b, st, info = mstep(D, a, b); t += st
        if info.get('drop'): drops[info['drop']] += 1
    return report('6 orbit', ok, N, '(classes %s; drops %s; t=%.3g, %d digits)' % (dict(cls), dict(drops), t, len(D)))

def check_induct(r, N):
    from collections import Counter
    ok = tot = 0; trans = Counter(); drops = Counter(); fails = []
    for it in range(N):
        cls = 'ABC'[it % 3]
        m = 1 + r.randint(0, 9) if r.random() < 0.8 else 1 + r.randint(0, 2)
        if cls == 'B' and m < 2: m = 2
        big = r.choice([5, 40, 300, 5000])
        st = rand_inv_state(r, cls, m, big)
        if st is None:
            continue
        D, a, b = st
        tot += 1
        try:
            name, D2, a2, b2, steps, info = mstep(D, a, b)
        except ValueError as e:
            fails.append((cls, D, a, b, str(e))); continue
        c2 = inv(D2, a2, b2)
        good = c2 is not None and info.get('drop') != 'LSB' and steps >= 1
        ok += good
        trans[(cls, c2)] += 1
        if info.get('drop'): drops[info['drop']] += 1
        if not good: fails.append((cls, D, a, b, name, D2, a2, b2))
    for f in fails[:5]: print('  induct FAIL', f)
    return report('7 induction', ok, tot, '(transitions %s; drops %s)' % (dict(sorted(trans.items(), key=str)), dict(drops)))

def check_induct_tm(r, N):
    ok = tot = 0
    for it in range(N):
        cls = 'ABC'[it % 3]
        m = 1 + r.randint(0, 5)
        if cls == 'B' and m < 2: m = 2
        st = rand_inv_state(r, cls, m, 12)
        if st is None: continue
        D, a, b = st
        if sum(u + v for (u, _, v) in D) + a + b > 400: continue
        name, D2, a2, b2, steps, info = mstep(D, a, b)
        s, h, q = render(D, a, b)
        res = run(s, h, q, steps)
        t, th, tq = render(D2, a2, b2)
        tot += 1; ok += same_cfg(res, t, th, tq)
    return report('8 induction-TM', ok, tot, '(random INV states, one macro step on the real machine)')

def check_drop_algebra(r, N):
    ok = tot = 0
    for _ in range(N):
        m = r.randint(2, 12); js = r.randint(1, m - 1)                  # j* in 1..m-1 (1-based)
        bits = [r.randint(0, 1) for _ in range(js - 1)] + [0] + [1] * (m - js)
        H = int(''.join(map(str, bits[:js - 1])) or '0', 2)
        V = int(''.join(map(str, bits)), 2)
        # after the carry, without the drop: bits[:js-1], 1, 0...0 ; with the drop: bits[:js-1], 0..0 (m-1)
        nd = bits[:js - 1] + [1] + [0] * (m - js)
        dr = bits[:js - 1] + [0] * (m - js)
        def RR1(bb, j):
            U = int(''.join(map(str, bb[:j])), 2); R = 2 ** j - 1 - U
            return R, (R + 1 - bb[j - 1]) // 2
        good = True
        for j in range(js + 1, m + 1):
            Rn, R1n = RR1(nd, j); Rp, R1p = RR1(dr, j - 1)
            good &= (Rn - Rp == 2 ** (j - js) * (2 ** (js - 1) - 1 - H))
            if j >= js + 2: good &= (2 * R1p == R1n + 2 ** (j - js - 1))
            else:           good &= (R1p == RR1(bits, js)[1])            # absorber: R1' = R1 of dropper
        Vp = int(''.join(map(str, dr)) or '0', 2)
        good &= (V + 1 - Vp == (H + 1) * 2 ** (m - js))                    # T grows by 2(H+1)2^(m-j*)
        Rtot0 = 2 ** m - 1 - V; Rtotp = 2 ** (m - 1) - 1 - Vp
        dS = (1 + Rtotp - 2 ** m + 2) - (Rtot0 - 2 ** (m + 1) + 2)
        good &= (dS == 2 ** (m - 1) + (H + 1) * 2 ** (m - js))
        good &= ((Rtotp - (Rtot0 - 1)) % 2 == 0)                           # a + Rtot parity kept
        tot += 1; ok += good
    return report('9 drop algebra', ok, tot, '(m=2..12, every j*)')

def check_twins():
    """Same table started in B / C on the blank tape (the two other holdouts with this table)."""
    ok = 0
    for st, t, D, a, b in [('B', 732, [(1, 1, 10)], 8, 1), ('C', 687, [(1, 1, 9)], 8, 1)]:
        res = run('', 0, st, t)
        ok += same_cfg(res, *render(D, a, b)) and inv(D, a, b) == 'A'
    return report('10 twins', ok, 2, '(start B: step 732 -> K([(1,1,10)];8,1); start C: step 687 -> K([(1,1,9)];8,1); both class A)')

# ---------------------------------------------------------------- M and O derived from lemmas (proof.md 2.1-2.3)
# Window lemmas: (lhs, rhs, steps). 'X>' = head in state X on the next cell, '<X' = head in state X on the
# previous cell. check_lemmas runs each lhs window exactly and requires that the head reads no cell outside the
# window before the last step. A lemma checked this way holds in every tape context, so for T1..T14 the check is
# a proof. S1..S6 are checked on one block; induction on n gives every n (proof.md 2.1).
LEM = {
    'S1': ('B>01', '10B>', 2), 'S2': ('C>10', '01C>', 2), 'S3': ('01<D', '<D10', 2), 'S4': ('10<E', '<E01', 2),
    'S5': ('C>1101', '10C>11', 6), 'S6': ('0B>10', '100B>', 6),
    'T1': ('B>00', '<E01', 3), 'T2': ('000<E', '100B>', 5), 'T3': ('100<E', '<E011', 3),
    'T4': ('01<E', '<D11', 2), 'T5': ('10<D', '<E01', 2), 'T6': ('00<D', '10C>', 3), 'T7': ('C>11', '00C>', 2),
    'T8': ('C>11011', '10000B>', 9), 'T9': ('0001001<E', '1001000B>', 13), 'T10': ('B>11', '00B>', 2),
    'T11': ('B>0', '1C>', 1), 'T12': ('C>11000', '<E01101', 13), 'T13': ('101<E', '<E011', 3),
    'T14': ('B>1101', '0010B>', 4),
}

def parse_window(w):
    i = w.find('>')
    if i > 0 and w[i - 1] in 'ABCDEF':
        return w[:i - 1] + w[i + 1:], i - 1, w[i - 1]
    i = w.find('<')
    return w[:i] + w[i + 2:], i - 1, w[i + 1]

def check_lemmas():
    ok = 0; fsteps = []
    for name, (lhs, rhs, n) in LEM.items():
        c, h, st = parse_window(lhs); c2, h2, st2 = parse_window(rhs)
        t = [int(x) for x in c]; s = ord(st) - 65; good = len(c) == len(c2)
        for _ in range(n):
            if not (0 <= h < len(t)) or TT[s][t[h]] is None:
                good = False; break
            if s == 5: fsteps.append((name, t[h]))
            w, d, s = TT[s][t[h]]; t[h] = w; h += d
        good = good and ''.join(map(str, t)) == c2 and h == h2 and s == ord(st2) - 65
        ok += good
        if not good: print('  lemma FAIL', name, lhs, rhs)
    fok = all(sym == 0 for (_, sym) in fsteps)
    return report('11 lemmas', ok + fok, len(LEM) + 1,
                  '(S1-S6 one block, T1-T14 exact; F is entered only in %s and reads 0 each time)'
                  % ','.join(sorted({nm for nm, _ in fsteps}, key=lambda x: int(x[1:]))))

def _dig(D):
    return ''.join('10' * u + ('' if be else '00') + '01' * v for (u, be, v) in D)

def M_chain(pre, us, vs, ones, a, b):
    """M on K(pre, (us,0,vs), (u,1,v) for (u,v) in ones; a, b) as the chain of proof.md 2.2.
    Returns cfg0, links, index of the first cell of d_{j*}'s 00. cfg = (left, mode, state, right) as strings,
    mode '>' = head on right[0], '<' = head on left[-1]; link = (lemma, repetitions, cfg after)."""
    L = '1001' + _dig(pre); P = '10' * us; i00 = len(L) + len(P)
    X = L + P + '00' + '01' * vs + ''.join('10' * u + '01' * v for (u, v) in ones)
    ln = []
    def A(name, n, l, mo, st, r): ln.append((name, n, (l, mo, st, r)))
    c0 = (X + '10' * a + '000' + '10' * b, '>', 'B', '')
    # phase 1: right-edge preamble with a dip into the 0^3 gap
    A('T1', 1, X + '10' * a + '000' + '10' * b, '<', 'E', '01')
    A('S4', b, X + '10' * a + '000', '<', 'E', '01' * (b + 1))
    A('T2', 1, X + '10' * a + '100', '>', 'B', '01' * (b + 1))
    A('S1', b + 1, X + '10' * (a + 1) + '0' + '10' * (b + 1), '>', 'B', '')
    # phase 2: one long left sweep to the 00 of d_{j*}
    A('T1', 1, X + '10' * (a + 1) + '0' + '10' * (b + 1), '<', 'E', '01')
    A('S4', b + 1, X + '10' * (a + 1) + '0', '<', 'E', '01' * (b + 2))
    A('T3', 1, X + '10' * a, '<', 'E', '011' + '01' * (b + 2))
    Rt = '01' * a + '011' + '01' * (b + 2)
    A('S4', a, X, '<', 'E', Rt)
    Y = L + P + '00' + '01' * vs
    for i in range(len(ones) - 1, -1, -1):                  # Lb on each bit-1 digit, right to left
        u, v = ones[i]
        Xi = Y + ''.join('10' * uu + '01' * vv for (uu, vv) in ones[:i])
        A('T4', 1, Xi + '10' * u + '01' * (v - 1), '<', 'D', '11' + Rt)
        A('S3', v - 1, Xi + '10' * u, '<', 'D', '10' * (v - 1) + '11' + Rt)
        A('T5', 1, Xi + '10' * (u - 1), '<', 'E', '01' + '10' * (v - 1) + '11' + Rt)
        Rt = '01' * u + '10' * (v - 1) + '11' + Rt
        A('S4', u - 1, Xi, '<', 'E', Rt)
    # Lf: flip d_{j*}
    A('T4', 1, L + P + '00' + '01' * (vs - 1), '<', 'D', '11' + Rt)
    A('S3', vs - 1, L + P + '00', '<', 'D', '10' * (vs - 1) + '11' + Rt)
    A('T6', 1, L + P + '10', '>', 'C', '10' * (vs - 1) + '11' + Rt)
    Z = L + P + '10' + '01' * (vs - 1)
    A('S2', vs - 1, Z, '>', 'C', '11' + Rt)
    # phase 3: rightward pass with 2-cell back-bounces (Rb on each swept digit), then the a block and the edge
    for i, (u, v) in enumerate(ones):
        rest = ''.join('01' * uu + '10' * (vv - 1) + '11' for (uu, vv) in ones[i + 1:]) + '01' * a + '011' + '01' * (b + 2)
        A('S5', u, Z + '10' * u, '>', 'C', '11' + '10' * (v - 1) + '11' + rest)
        A('T7', 1, Z + '10' * u + '00', '>', 'C', '10' * (v - 1) + '11' + rest)
        Z = Z + '10' * u + '00' + '01' * (v - 1)
        A('S2', v - 1, Z, '>', 'C', '11' + rest)
    A('S5', a, Z + '10' * a, '>', 'C', '11' + '011' + '01' * (b + 2))
    A('T8', 1, Z + '10' * a + '10000', '>', 'B', '01' * (b + 2))
    A('S1', b + 2, Z + '10' * (a + 1) + '000' + '10' * (b + 2), '>', 'B', '')
    return c0, ln, i00

def O_chain(ones, a, b):
    """O on the all-ones K((u,1,v) for (u,v) in ones; a, b) as the chain of proof.md 2.3. Left strings carry a
    3-cell '000' margin before the leading 10 01, so that T9 keeps the coordinates."""
    m = len(ones)
    X = '000' + '1001' + ''.join('10' * u + '01' * v for (u, v) in ones)
    ln = []
    def A(name, n, l, mo, st, r): ln.append((name, n, (l, mo, st, r)))
    c0 = (X + '10' * a + '000' + '10' * b, '>', 'B', '')
    A('T1', 1, X + '10' * a + '000' + '10' * b, '<', 'E', '01')                        # phase 1
    A('S4', b, X + '10' * a + '000', '<', 'E', '01' * (b + 1))
    A('T2', 1, X + '10' * a + '100', '>', 'B', '01' * (b + 1))
    A('S1', b + 1, X + '10' * (a + 1) + '0' + '10' * (b + 1), '>', 'B', '')
    A('T1', 1, X + '10' * (a + 1) + '0' + '10' * (b + 1), '<', 'E', '01')              # phase 2
    A('S4', b + 1, X + '10' * (a + 1) + '0', '<', 'E', '01' * (b + 2))
    A('T3', 1, X + '10' * a, '<', 'E', '011' + '01' * (b + 2))
    Rt = '01' * a + '011' + '01' * (b + 2)
    A('S4', a, X, '<', 'E', Rt)
    for i in range(m - 1, -1, -1):
        u, v = ones[i]
        Xi = '000' + '1001' + ''.join('10' * uu + '01' * vv for (uu, vv) in ones[:i])
        A('T4', 1, Xi + '10' * u + '01' * (v - 1), '<', 'D', '11' + Rt)
        A('S3', v - 1, Xi + '10' * u, '<', 'D', '10' * (v - 1) + '11' + Rt)
        A('T5', 1, Xi + '10' * (u - 1), '<', 'E', '01' + '10' * (v - 1) + '11' + Rt)
        Rt = '01' * u + '10' * (v - 1) + '11' + Rt
        A('S4', u - 1, Xi, '<', 'E', Rt)
    A('T9', 1, '1001000', '>', 'B', Rt)                                                # phase 3: far-left turn
    Z = '1001000'
    for i, (u, v) in enumerate(ones):                                                  # phase 4: Ro on each digit
        rest = ''.join('01' * uu + '10' * (vv - 1) + '11' for (uu, vv) in ones[i + 1:]) + '01' * a + '011' + '01' * (b + 2)
        A('S1', u, Z + '10' * u, '>', 'B', '10' * (v - 1) + '11' + rest)
        A('S6', v - 1, Z + '10' * (u - 1) + '1' + '10' * (v - 1) + '0', '>', 'B', '11' + rest)
        Z = Z + '10' * (u - 1) + '1' + '10' * (v - 1) + '000'
        A('T10', 1, Z, '>', 'B', rest)
    A('S1', a, Z + '10' * a, '>', 'B', '011' + '01' * (b + 2))                         # phase 5: right end
    A('T11', 1, Z + '10' * a + '1', '>', 'C', '11' + '01' * (b + 2))
    A('S5', b + 2, Z + '10' * a + '1' + '10' * (b + 2), '>', 'C', '11')
    A('T12', 1, Z + '10' * a + '1' + '10' * (b + 2), '<', 'E', '01101')
    A('S4', b + 2, Z + '10' * a + '1', '<', 'E', '01' * (b + 2) + '01101')
    A('T13', 1, Z + '10' * (a - 1), '<', 'E', '011' + '01' * (b + 2) + '01101')        # phase 6: back over a
    A('S4', a - 1, Z, '<', 'E', '01' * (a - 1) + '011' + '01' * (b + 2) + '01101')
    Zp = Z[:-3] + '100'                                                                # phase 7: turn, final run
    A('T2', 1, Zp, '>', 'B', '01' * (a - 1) + '011' + '01' * (b + 2) + '01101')
    A('S1', a, Zp + '10' * a, '>', 'B', '1' + '01' * (b + 2) + '01101')
    A('S6', b + 3, Zp + '10' * (a - 1) + '1' + '10' * (b + 3) + '0', '>', 'B', '1101')
    A('T14', 1, Zp + '10' * (a - 1) + '1' + '10' * (b + 3) + '0' + '0010', '>', 'B', '')
    return c0, ln

def run_chain(c0, ln, pad=64):
    """Run the real machine link by link from c0; after each link the whole tape, head and state must equal
    the link's configuration. Returns (ok, steps, leftmost head index, final nonblank tape)."""
    l, mo, st, r = c0
    n = sum(LEM[nm][2] * k for (nm, k, _) in ln)
    t = bytearray(pad) + bytearray(int(x) for x in l + r) + bytearray(pad + n)
    h = pad + (len(l) if mo == '>' else len(l) - 1); s = ord(st) - 65; mn = h; tot = 0
    for (nm, k, (l2, mo2, st2, r2)) in ln:
        for _ in range(LEM[nm][2] * k):
            tr = TT[s][t[h]]
            if tr is None:
                return False, tot, None, None
            t[h] = tr[0]; h += tr[1]; s = tr[2]; mn = min(mn, h); tot += 1
        want = bytearray(pad) + bytearray(int(x) for x in l2 + r2)
        want += bytearray(len(t) - len(want))
        if t != want or h != pad + (len(l2) if mo2 == '>' else len(l2) - 1) or s != ord(st2) - 65:
            return False, tot, None, None
    return True, tot, mn - pad, ''.join(map(str, t)).strip('0')

def _M_case(pre, us, vs, ones, a, b):
    c0, ln, i00 = M_chain(pre, us, vs, ones, a, b)
    good, steps, mn, end = run_chain(c0, ln)
    D = pre + [(us, 0, vs)] + [(u, 1, v) for (u, v) in ones]
    _, D2, a2, b2, st, _ = mstep(D, a, b)
    return good and steps == st and mn == i00 and end == render(D2, a2, b2)[0].strip('0')

def _O_case(ones, a, b):
    c0, ln = O_chain(ones, a, b)
    good, steps, mn, end = run_chain(c0, ln)
    _, D2, a2, b2, st, _ = mstep([(u, 1, v) for (u, v) in ones], a, b)
    return good and steps == st and mn == 0 and end == render(D2, a2, b2)[0].strip('0')

def check_M_chain(r, quick):
    """box: every link of the 2.2 chain on the real machine; the end is mstep's K-config, the total is the
    step formula, and the leftmost cell read is the first 0 of d_{j*}'s 00. a = 0 is included."""
    import itertools
    ok = tot = 0
    pres = [[], [(1, 1, 1)], [(0, 0, 2)], [(2, 1, 3), (1, 0, 1)]]
    kmax = 2 if quick else 3
    for pre in pres:
        for us in range(3):
            for vs in range(1, 4):
                for k in range(kmax + 1):
                    for ones in itertools.product([(u, v) for u in (1, 2) for v in (1, 2)], repeat=k):
                        for a in range(4):
                            for b in range(3):
                                tot += 1; g = _M_case(pre, us, vs, list(ones), a, b); ok += g
                                if not g: print('  M chain FAIL', pre, us, vs, ones, a, b)
    box = tot
    for _ in range(100 if quick else 500):
        pre = [rand_digit(r, r.randint(0, 1), 30) for _ in range(r.randint(0, 3))]
        ones = [(rnd(r, 1, 60), rnd(r, 1, 60)) for _ in range(r.randint(0, 7))]
        args = (pre, rnd(r, 0, 60), rnd(r, 1, 60), ones, rnd(r, 0, 300), rnd(r, 0, 300))
        tot += 1; g = _M_case(*args); ok += g
        if not g: print('  M chain FAIL', args)
    return report('12 M chain', ok, tot, '(box: %d configs, u*<=2, v*<=3, k<=%d, u,v<=2, a<=3, b<=2, 4 left contexts; '
                  'plus %d random, k<=7, values up to 360)' % (box, kmax, tot - box))

def check_O_chain(r, quick):
    import itertools
    ok = tot = 0
    mmax = 3 if quick else 4
    for m in range(1, mmax + 1):
        rng = (1, 2, 3) if m <= 3 else (1, 2)
        for ones in itertools.product([(u, v) for u in rng for v in rng], repeat=m):
            for a in range(1, 4):
                for b in range(3):
                    tot += 1; g = _O_case(list(ones), a, b); ok += g
                    if not g: print('  O chain FAIL', ones, a, b)
    box = tot
    for _ in range(60 if quick else 300):
        ones = [(rnd(r, 1, 60), rnd(r, 1, 60)) for _ in range(r.randint(1, 8))]
        args = (ones, rnd(r, 1, 300), rnd(r, 0, 300))
        tot += 1; g = _O_case(*args); ok += g
        if not g: print('  O chain FAIL', args)
    return report('13 O chain', ok, tot, '(box: %d configs, m<=%d, u,v<=3 (<=2 at m=4), a<=3, b<=2; plus %d random, m<=8)'
                  % (box, mmax, tot - box))


if __name__ == '__main__':
    ap = argparse.ArgumentParser()
    ap.add_argument('--seed', type=int, default=11)
    ap.add_argument('--orbit', type=float, default=3e5)
    ap.add_argument('--induct', type=float, default=1.5e5)
    ap.add_argument('--tm-steps', type=float, default=3e7)
    ap.add_argument('--quick', action='store_true')
    A = ap.parse_args()
    r = random.Random(A.seed)
    if A.quick:
        A.orbit, A.induct, A.tm_steps = 2e4, 2e4, 2e6
    res = [check_sweeps(),
           check_M(r, 700 if not A.quick else 140, 300),
           check_O(r, 350 if not A.quick else 70, 300),
           check_init(),
           check_closure(int(A.tm_steps)),
           check_orbit(int(A.orbit)),
           check_induct(r, int(A.induct)),
           check_induct_tm(r, 3000 if not A.quick else 300),
           check_drop_algebra(r, 20000),
           check_twins(),
           check_lemmas(),
           check_M_chain(r, A.quick),
           check_O_chain(r, A.quick)]
    print('ALL PASS' if all(res) else 'SOME CHECKS FAILED')
    sys.exit(0 if all(res) else 1)