#!/usr/bin/env python3
"""check_proof.py for W7-5, table 1RB0RD_0LC0RA_1LA1LC_1RE0LF_0RA1RF_1LB--- (740 lines 274 and 514).

Validates every claim of proof.md.  Usage:
    python3 -B check_proof.py              # all sections (about a minute)
    python3 -B check_proof.py A C F        # selected sections

Sections (proof.md section in brackets)
  A  table facts: both strings byte-exact on the 740 list, F1 is the only undefined transition, line 514 is
     this table mirrored and renamed (12/12 transitions), started in F                               [1, 7.2]
  B  executor soundness: every pass lemma re-run by plain simulation on concrete instances with random
     context cells (never read), exact step counts; the D-record lemmas reach no earlier D-record; S0
     cell by cell (it visits only the four cells of 1010)                                                [3]
  C  the 14 pass lemmas proved for all values of the block variables by symbolic execution (symx.py)   [3]
  D  sweep theorem: the abstract sweep (sweep.py, built from the lemmas) = the machine, configuration and
     step count, on random and on all small in-domain configurations                                    [4]
  E  counter and consumption laws of the sweep (phases, wall, consumers, which field pays, r and the
     lemma that ends the pass), and the O2 gap telescoping against the machine                          [4.3, 4.4]
  F  the epoch theorem: arithmetic (closed forms, the inequality) and full abstract epochs from
     adversarial birth configurations, k = 2..7, with the number of block groups at the O2 sweep        [5.2]
  G  the two openings (274: blank -> B(3; 0,41) at step 7,910; 514: blank -> mirror of B(2;7) at step 792), and
     the abstract orbit = every real D-record to 1e9 on both starts, plus the first D-record after 1e11  [7]
  H  agreement with the observe-stage model (scripts/model.py sweep + step rule R3), informational       [8]
"""
import sys, os, random, itertools, time, collections
HERE = os.path.dirname(os.path.abspath(__file__))
REPO = os.path.abspath(os.path.join(HERE, '..', '..'))
sys.path.insert(0, os.path.join(HERE, 'proof_scripts'))
import symx
import sweep as sw
import tm as T

TM = '1RB0RD_0LC0RA_1LA1LC_1RE0LF_0RA1RF_1LB---'
TM514 = '1RB---_0RC0LD_1RD1RC_1LB0LE_1LF0RA_0LD1LA'
FAIL = []
X, P0 = '10101010', '11001100'


def check(cond, msg):
    print(('  ok   ' if cond else '  FAIL ') + msg, flush=True)
    if not cond:
        FAIL.append(msg)


# ------------------------------------------------------------------ the lemmas of proof.md section 3
LEMMAS = [
    ('S0', '100 D> 0^inf', '100 <B 1010 0^inf', '9'),
    ('L1', '100 <B', '<B 100', '5'),
    ('L2', '10000 <B (X)^j 100', '<B (X)^j+1', '19+20j'),
    ('L3', '11001100 <B', '110000 10 A>', '9'),
    ('L4', '1001100 <B', '10000 10 A>', '9'),
    ('L5', '110000 <B (X)^j 100 (X)^i 100', '1001100 (P0)^i+j 100 10 A>', '59+56j+36i'),
    ('L5e', '110000 <B (X)^j 100 (X)^i 1010 0^inf', '1001100 (P0)^i+j 100 100 D> 0^inf', '62+56j+36i'),
    ('L6', '1100 <B', '00 10 A>', '9'),
    ('L7', '0^inf <B (X)^j 100', '0^inf 1100 (P0)^j 100 1 E>', '38+56j'),
    ('RA1', '10 A> (X)^j 100', '(P0)^j 100 10 A>', '3+36j'),
    ('RA2', '10 A> 1010 0^inf', '100 100 D> 0^inf', '6'),
    ('RA3', '10 A> (X)^j 1010 0^inf', '(P0)^j 100 100 D> 0^inf', '6+36j'),
    ('RE1', '100 1 E> (X)^j 100', '(P0)^j 100 100 1 E>', '5+36j'),
    ('RE2', '100 100 1 E> 1010 0^inf', '10000 100 100 100 D> 0^inf', '49'),
]


def subst_words(s):
    return s.replace('(X)', f'({X})').replace('(P0)', f'({P0})')


def concrete(s, env):
    """notation -> (cells, head index, state, left end '0'/'?', right end '0'/'?').  Head index may be -1 or
    len(cells) (on a context cell)."""
    s = s.strip()
    lend = rend = '?'
    if s.startswith('0^inf'):
        lend, s = '0', s[5:].strip()
    if s.endswith('0^inf'):
        rend, s = '0', s[:-5].strip()
    cells, head, q = '', None, None
    for w in s.split():
        if len(w) == 2 and w[1] == '>' and w[0] in 'ABCDEF':
            head, q = len(cells), w[0]
        elif len(w) == 2 and w[0] == '<' and w[1] in 'ABCDEF':
            head, q = len(cells) - 1, w[1]
        elif w.startswith('('):
            word, _, e = w[1:].partition(')^')
            n = symx.parse_expr(e)
            val = n.get('', 0) + sum(c * env[v] for v, c in n.items() if v)
            cells += word * val
        else:
            cells += w
    return cells, head, q, lend, rend


def run_concrete(lhs, rhs, steps, env, rnd, drec_end):
    """Plain simulation of one lemma instance with random context; returns '' or an error string."""
    c1, h1, q1, le, re = concrete(lhs, env)
    c2, h2, q2, le2, re2 = concrete(rhs, env)
    PADL = PADR = 40
    ctxL = [0] * PADL if le == '0' else [rnd.randint(0, 1) for _ in range(PADL)]
    ctxR = [0] * PADR if re == '0' else [rnd.randint(0, 1) for _ in range(PADR)]
    t = ctxL + [int(c) for c in c1] + ctxR
    pos = PADL + h1
    q = q1
    lo_ok = 0 if le == '0' else PADL          # cells < lo_ok are context (must not be read)
    hi_ok = len(t) if re == '0' else PADL + len(c1)
    hi_vis = pos if drec_end or lhs.startswith('100 D>') else None
    if lhs.strip().endswith('0^inf'):
        hi_vis = PADL + len(c1) - 1 if not lhs.startswith('100 D>') else pos
    early = 0
    for n in range(steps):
        if not lo_ok <= pos < hi_ok:
            return f'reads context at step {n}'
        tr = T.T274[q, t[pos]]
        if tr is None:
            return 'halts'
        a, mv, q = tr
        t[pos] = a
        pos += mv
        if hi_vis is not None and pos > hi_vis:
            hi_vis = pos
            if q == 'D' and n < steps - 1:
                early += 1
    if early:
        return f'{early} earlier D-record(s)'
    # compare: windows are aligned at their left ends, except when the left end is 0^inf (L7), where the
    # window grows to the left and the right ends are aligned
    off = PADL + (len(c1) - len(c2) if le == '0' else 0)
    want = [int(c) for c in c2]
    got = t[off:off + len(c2)]
    if got != want:
        return 'tape differs'
    if le == '0':
        if any(t[:off]):
            return 'nonzero left of the window'
    elif t[:PADL] != ctxL:
        return 'left context changed'
    if re == '0':
        if any(t[off + len(c2):]):
            return 'nonzero beyond the window'
    elif le == '0':
        if t[off + len(c2):] != ctxR:
            return 'right context changed'
    elif t[off + len(c2):] != ctxR[len(c2) - len(c1):]:
        return 'right context changed'
    if pos != off + h2 or q != q2:
        return f'head {pos - off}/{q} vs {h2}/{q2}'
    if drec_end and not (q == 'D' and pos == hi_vis):
        return 'end is not a D-record'
    return ''


def sec_A():
    print('A  table facts')
    lines = open(os.path.join(REPO, 'data/holdouts/BB6_holdouts_740.txt')).read().split('\n')
    check(lines[273] == TM, f'line 274 of the 740 list is {TM} byte for byte')
    check(lines[513] == TM514, f'line 514 of the 740 list is {TM514} byte for byte')
    check(lines.count(TM) == 1 and lines.count(TM514) == 1, 'each string occurs exactly once on the list')
    t, t5 = T.table(TM), T.table(TM514)
    undef = [k for k, v in t.items() if v is None]
    check(undef == [('F', 1)], f'the only undefined transition of 274 is F1 ({undef})')
    entering_F = sorted(k for k, v in t.items() if v and v[2] == 'F')
    check(entering_F == [('D', 1), ('E', 1)], f'F is entered only by D1 = 0LF and E1 = 1RF ({entering_F})')
    phi = dict(A='F', B='B', C='C', D='A', E='D', F='E')     # 514 state -> 274 state
    inv = {v: k for k, v in phi.items()}
    ok = 0
    for (s, b), v in t5.items():
        w = t[phi[s], b]
        if v is None and w is None:
            ok += 1
        elif v and w and v[0] == w[0] and v[1] == -w[1] and v[2] == inv[w[2]]:
            ok += 1
    check(ok == 12, f'514 = 274 mirrored with states renamed by phi = A->F,B->B,C->C,D->A,E->D,F->E: {ok}/12 '
          'transitions (so 514 from blank = the mirror of 274\'s table from 0^inf F> 0^inf, phi on states)')
    others = [m for m in itertools.permutations('ABCDEF') for mir in (0, 1)
              if all((t5[s, b] is None and t[m['ABCDEF'.index(s)], b] is None) or
                     (t5[s, b] and t[m['ABCDEF'.index(s)], b] and t5[s, b][0] == t[m['ABCDEF'.index(s)], b][0]
                      and t5[s, b][1] == (-1 if mir else 1) * t[m['ABCDEF'.index(s)], b][1]
                      and t5[s, b][2] == 'ABCDEF'[m.index(t[m['ABCDEF'.index(s)], b][2])])
                     for s in 'ABCDEF' for b in (0, 1))]
    check(len(others) == 1, f'phi with mirror is the only map among 720 x 2 relating the two tables ({len(others)})')


def sec_B(n_inst=6):
    print('B  executor soundness: every lemma by plain simulation on concrete instances')
    rnd = random.Random(7)
    for nm, l, r, s in LEMMAS:
        lhs, rhs = subst_words(l), subst_words(r)
        vars_ = sorted({v for v in 'ij' if f')^{v}' in lhs or f')^{v}+' in lhs or f'+{v}' in lhs})
        errs, cnt = [], 0
        for vals in itertools.product(range(n_inst), repeat=len(vars_)):
            env = dict(zip(vars_, vals))
            e = symx.parse_expr(s)
            steps = e.get('', 0) + sum(c * env[v] for v, c in e.items() if v)
            for _ in range(3):
                msg = run_concrete(lhs, rhs, steps, env, rnd, drec_end='D> 0^inf' in r)
                cnt += 1
                if msg:
                    errs.append((env, msg))
        check(not errs, f'{nm}: {cnt} concrete instances ({"x".join(vars_) or "no variables"} in 0..{n_inst - 1}), '
              f'3 random contexts each: exact, context never read{", no earlier D-record" if "D> 0^inf" in r else ""}'
              + (f'  {errs[:2]}' if errs else ''))
    # proof.md 3 (Review log, gap 1): S0 visits exactly the four cells of 1010, so every cell to their right is
    # still unvisited and blank when the sweep's last lemma (RA2, RA3, L5e or RE2) reaches it
    t = [1, 0, 0] + [0] * 12
    pos, q, h0, read = 3, 'D', 3, []
    for _ in range(9):
        a, mv, q = T.T274[q, t[pos]]
        read.append(pos - h0)
        t[pos] = a
        pos += mv
    check(sorted(set(read)) == [0, 1, 2, 3] and (pos - h0, q) == (-1, 'B') and t[h0:h0 + 4] == [1, 0, 1, 0]
          and not any(t[h0 + 4:]),
          f'S0 cell by cell: it reads only the four cells of 1010 (offsets {read} from the D-record cell), writes '
          f'1010 there and ends in B on the last cell of the 100 before them; the cell after 1010 is never visited '
          f'by S0 and is blank')


SHIFTS = [   # (name, object, head offset in the object, state, direction, X -> word, steps per copy)
    ('Ch1', '010', 0, 'A', 'R', X, 8),
    ('Ch2', '101010', 5, 'B', 'L', X, 12),
    ('Ch3', '100101010', 1, 'F', 'R', '10011001', 36),
    ('Ch4', '1101010', 0, 'B', 'R', '10011001', 36),
    ('Ch5', '010101010', 1, 'A', 'R', '10011001', 36),
    ('Ch6', '00110101010', 3, 'E', 'R', '10011001', 36),
]


def sec_C():
    print('C  the shift rules and the pass lemmas, for all values of the variables')
    rnd = random.Random(13)
    for nm, o, d, q, dr, wp, per in SHIFTS:
        bad = []
        for n in range(0, 9):
            for _ in range(3):
                ctxL = [rnd.randint(0, 1) for _ in range(30)]
                ctxR = [rnd.randint(0, 1) for _ in range(30)]
                if dr == 'R':      # o X^n  -->  wp^n o, head at offset d of o
                    win = [int(c) for c in o + X * n]
                    want = [int(c) for c in wp * n + o]
                    pos0, pos1 = d, len(X) * n + d
                else:              # X^n o  -->  o wp^n
                    win = [int(c) for c in X * n + o]
                    want = [int(c) for c in o + wp * n]
                    pos0, pos1 = len(X) * n + d, d
                t = ctxL + win + ctxR
                pos, st = 30 + pos0, q
                for _k in range(per * n):
                    if not 30 <= pos < 30 + len(win):
                        bad.append((n, 'context read'))
                        break
                    a, mv, st = T.T274[st, t[pos]]
                    t[pos] = a
                    pos += mv
                if t[30:30 + len(win)] != want or pos != 30 + pos1 or st != q or t[:30] != ctxL or t[30 + len(win):] != ctxR:
                    bad.append((n, 'result'))
        lhs = (f'{o[:d]} {q}> {o[d:]} X^n' if dr == 'R' else f'X^n {o[:d]} {q}> {o[d:]}')
        rhs = (f'{"X" if wp == X else "(" + wp + ")"}^n {o[:d]} {q}> {o[d:]}' if dr == 'R'
               else f'{o[:d]} {q}> {o[d:]} {"X" if wp == X else "(" + wp + ")"}^n')
        check(not bad, f'{nm}: {lhs}  -->({per}n)  {rhs}: n = 0..8 by plain simulation, random context never read '
              f'(the executor proves the rule for all n from n = 1 by induction)' + (f' {bad[:2]}' if bad else ''))
    ok = all('1001100' + P0 * n + '1' == '10011001' * (n + 1) for n in range(40))
    check(ok, 'word identity P2 P0^n 1 = (10011001)^(n+1)  (n = 0..39; by induction it holds for all n)')
    print('   pass lemmas (symbolic execution):')
    for nm, l, r, s in LEMMAS:
        res = symx.lemma(subst_words(l), subst_words(r), s)
        ok = all(x[1] for x in res)
        chains = sorted({ln.split('  [')[0] for x in res for ln in x[3] if ln.startswith('chain')})
        check(ok, f'{nm}: {l}  -->({s})  {r}   [{len(res)} branch(es); chains: {len(chains)}]'
              + ('' if ok else f'  {[x[2] for x in res if not x[1]][:1]}'))


def tm_sweep(cfg, maxsteps=10 ** 6):
    s, h = sw.encode(cfg)
    res = T.to_drecord(s, h, maxsteps=maxsteps)
    if res[0] in ('HALT', 'MAX'):
        return None, res
    n, ts, hh, _, _ = res
    return T.decode(ts[:hh]), n


def rand_cfg(rnd, kmax=6, fmax=12, rmax=15):
    k = rnd.randint(1, kmax)
    m = rnd.random()
    if m < 0.3:
        p = [3] * k
    elif m < 0.6:
        j = rnd.randint(0, k - 1)
        p = [rnd.randint(0, 3) for _ in range(j)] + [3] * (k - j)
        if j:
            p[j - 1] = rnd.randint(0, 2)
    else:
        p = [rnd.randint(0, 3) for _ in range(k)]
    f = [rnd.choice([0, 0, 0, 1, 2, rnd.randint(0, fmax)]) for _ in range(k - 1)]
    return (rnd.randint(0, 1), p, f, rnd.randint(1, rmax))


def sec_D(n_rand=6000):
    print('D  sweep theorem: abstract sweep (pass lemmas composed) = the machine')
    rnd = random.Random(11)
    c = collections.Counter()
    bad = []
    outside = collections.Counter()
    for _ in range(n_rand):
        cfg = rand_cfg(rnd)
        try:
            nxt, steps, ev, _ = sw.sweep(cfg)
        except sw.Domain:
            outside['domain'] += 1
            continue
        got, n = tm_sweep(cfg)
        ok = got is not None and tuple(got) == nxt and n == steps
        c[ev] += 1
        if not ok:
            bad.append((cfg, ev))
    check(not bad, f'random in-domain configurations: {sum(c.values())} sweeps {dict(c)} exact (configuration and '
          f'steps); {outside["domain"]} random configurations were outside the domain' + (f' {bad[:2]}' if bad else ''))
    # exhaustive small
    c = collections.Counter()
    bad = []
    nout = 0
    for W in (0, 1):
        for k in range(1, 5):
            for p in itertools.product(range(4), repeat=k):
                for f in itertools.product(range(3), repeat=k - 1):
                    for r in range(1, 6):
                        cfg = (W, list(p), list(f), r)
                        try:
                            nxt, steps, ev, _ = sw.sweep(cfg)
                        except sw.Domain:
                            nout += 1
                            continue
                        got, n = tm_sweep(cfg)
                        c[ev] += 1
                        if not (got is not None and tuple(got) == nxt and n == steps):
                            bad.append(cfg)
    check(not bad, f'all configurations with k <= 4, fields 0..2, r 1..5, both walls: {sum(c.values())} in-domain '
          f'sweeps exact {dict(c)}; {nout} outside the domain (not claimed)' + (f' {bad[:3]}' if bad else ''))


def o2_by_gaps(cfg):
    """proof.md 4.3 (E-pass fields), independent of sweep.py.  cfg: no wall, every phase 3.  Wraps: each digit,
    right to left, takes one s from the first nonempty field at or right of it (h = fields after the wraps,
    h[k-1] = the rest of r).  T = X s^h0 X s^h1 ... X s^h(k-1) end, so the maximal block groups are separated by
    the nonzero h, and the gaps g_1..g_m are the nonzero h left of the last digit, then h[k-1].  L7 keeps the s
    after group 1, each RE1 on a group i >= 2 moves one s from gap i-1 to gap i, RE2 takes two s from gap m:
    g_1 - 1, g_2, ..., g_(m-1), g_m - 1  (m >= 2)   or   g_1 - 2  (m = 1).
    Domain: (D1) every wrap finds an s; (D2) g_m >= 1, and g_1 >= 2 if m = 1.
    Returns (next configuration, m, h) or None outside the domain."""
    W, p, f, r = cfg
    assert W == 0 and all(x == 3 for x in p)
    k = len(p)
    h = list(f) + [r]
    for d in range(k - 1, -1, -1):
        x = next((x for x in range(d, k) if h[x] > 0), None)
        if x is None:
            return None
        h[x] -= 1
    gi = [x for x in range(k - 1) if h[x] > 0] + [k - 1]
    g = [h[x] for x in gi]
    m = len(g)
    if g[-1] < 1 or (m == 1 and g[0] < 2):
        return None
    g2 = [g[0] - 2] if m == 1 else [g[0] - 1] + g[1:-1] + [g[-1] - 1]
    h2 = list(h)
    for x, v in zip(gi, g2):
        h2[x] = v
    return (1, [0] * k + [3], h2, 3), m, h


def sec_E(n=4000):
    print('E  counter and consumption laws of the sweep')
    rnd = random.Random(5)
    bad = collections.Counter()
    ends = collections.Counter()
    cnt = 0
    for _ in range(n):
        cfg = rand_cfg(rnd, kmax=8)
        W, p, f, r = cfg
        try:
            (W2, p2, f2, r2), steps, ev, info = sw.sweep(cfg)
        except sw.Domain:
            continue
        cnt += 1
        k = len(p)
        V, V2 = sw.value(p), sw.value(p2)
        stop = info['stop']
        if ev in ('R1', 'R2'):
            if not (len(p2) == k and V2 == V + 1 and W2 == W):
                bad['increment'] += 1
            want = {stop} if p[stop] == 1 else set()
            want |= set(range(stop + 1, k))
        elif ev == 'O1':
            if not (W == 1 and W2 == 0 and p2 == [0] * k and all(x == 3 for x in p)):
                bad['O1'] += 1
            want = set(range(k))
        else:
            if not (W == 0 and W2 == 1 and p2 == [0] * k + [3] and all(x == 3 for x in p)):
                bad['O2'] += 1
            want = set(range(k))
        if {d for d, _ in info['cons']} != want or len(info['cons']) != len(want):
            bad['consumers'] += 1
        # which field pays: the first field at or right of the digit that still has an s at that moment
        fields = list(f) + [r]
        for d, fld in info['cons']:          # in the order the left pass makes them (right to left)
            first = next((x for x in range(d, k) if fields[x] > 0), None)
            if first != fld:
                bad['payer'] += 1
            fields[fld] -= 1
        # fields after a non-O2 sweep: the stop digit and the digits left of it keep their fields minus what
        # was taken; the order of digits is kept; r gains 2
        if ev != 'O2':
            want_f = fields[:-1]
            if f2 != want_f or r2 != fields[-1] + 2:
                bad['fields'] += 1
            ends[info['end']] += 1
            if info['end'] not in ('RA2', 'RA3', 'L5e') or (fields[-1] == 0) != (info['end'] != 'RA2'):
                bad['ending'] += 1
        else:
            # O2 (proof.md 4.4): the first nonzero field and the rest of r each lose one s (two if they are the
            # same field); the rest of r becomes the newborn's left field; r = 3; sea count + 1
            nz = [x for x in range(k) if fields[x] > 0]
            want_f = list(fields)
            want_f[nz[0]] -= 1
            want_f[-1] -= 1
            if r2 != 3 or f2 != want_f:
                bad['O2 fields'] += 1
            if sum(f2) + r2 != sum(fields) + 1:
                bad['O2 sea count'] += 1
        # ids: the digits keep their order
        ids = [x for x in info['ids'] if x != 'new']
        if ids != list(range(k)):
            bad['order'] += 1
    check(not bad, f'{cnt} in-domain random sweeps: V -> V+1 (R1/R2), O1: wall off and all phases 0, O2: wall on, '
          'all 0 and a newborn 3; consumers = {stop digit if its phase is 1} + wrapped digits; each takes its s '
          'from the first nonempty field at or right of it (right to left); other fields unchanged; r += 2 '
          '(non-O2); O2: first nonzero field and rest of r lose one s each, rest of r becomes the newborn\'s '
          'left field, r = 3, sea count + 1; digit order kept'
          + (f' {dict(bad)}' if bad else ''))
    # proof.md 4.4 (Review log, gap 3): r' = (rest of r after consumption) + 2 for every ending of an A-pass,
    # on the random sweeps above and on every small configuration (no machine needed: section D ties sweep.py
    # to the machine)
    for W in (0, 1):
        for k in range(1, 5):
            for p in itertools.product(range(4), repeat=k):
                for f in itertools.product(range(3), repeat=k - 1):
                    for r in range(1, 6):
                        try:
                            (W2, p2, f2, r2), steps, ev, info = sw.sweep((W, list(p), list(f), r))
                        except sw.Domain:
                            continue
                        if ev == 'O2':
                            continue
                        rest = r - sum(1 for _, fld in info['cons'] if fld == k - 1)
                        ends[info['end']] += 1
                        if r2 != rest + 2 or (rest == 0) != (info['end'] != 'RA2'):
                            bad['ending small'] += 1
    check(not bad and set(ends) == {'RA2', 'RA3', 'L5e'},
          f'non-O2 sweeps by their last lemma {dict(ends)}: r\' = rest of r + 2 in every one; RA2 ends exactly the '
          f'sweeps whose rest of r is >= 1, RA3 and L5e exactly those whose rest is 0 (r\' = 2)'
          + (f' {dict(bad)}' if bad else ''))
    # proof.md 4.3 (Review log, gap 2): the E-pass gap telescoping, checked against the machine itself
    rnd = random.Random(17)
    cfgs = [(0, [3] * k, list(f), r) for k in range(1, 6) for f in itertools.product(range(4), repeat=k - 1)
            for r in range(1, 9)]
    for _ in range(1500):
        k = rnd.randint(1, 8)
        cfgs.append((0, [3] * k, [rnd.choice([0, 0, 0, 1, 2, rnd.randint(0, 12)]) for _ in range(k - 1)],
                     rnd.randint(1, 15)))
    cm = collections.Counter()
    bad = collections.Counter()
    for cfg in cfgs:
        g = o2_by_gaps(cfg)
        try:
            nxt, steps, ev, info = sw.sweep(cfg)
        except sw.Domain:
            nxt = None
        if (g is None) != (nxt is None):
            bad['domain'] += 1
            continue
        if g is None:
            cm['outside'] += 1
            continue
        want, m, h = g
        cm['m = 1' if m == 1 else 'm >= 2'] += 1
        got, n_ = tm_sweep(cfg)
        if got is None or tuple(got) != want:
            bad['machine'] += 1
        if nxt != want or ev != 'O2':
            bad['sweep.py'] += 1
        law = list(h)                      # proof.md 4.4 wording: first nonzero field after the wraps, rest of r
        law[next(x for x in range(len(h)) if h[x] > 0)] -= 1
        law[-1] -= 1
        if law != want[2]:
            bad['law'] += 1
    check(not bad and cm['m = 1'] and cm['m >= 2'],
          f'O2 fields by gaps (4.3: L7 keeps the s after group 1, each RE1 on group i >= 2 moves one s from gap i-1 '
          f'to gap i, RE2 takes two from gap m -> g_1 - 1, g_2..g_(m-1), g_m - 1, or g_1 - 2 if m = 1; domain g_m >= 1 '
          f'and g_1 >= 2 if m = 1) = the machine\'s next D-record = sweep.py = the 4.4 law, on all O2 '
          f'configurations with k <= 5, fields 0..3, r 1..8 and 1,500 random ones up to k = 8: {dict(cm)}; '
          f'the gap domain equals sweep.py\'s domain' + (f' {dict(bad)}' if bad else ''))


def epoch_consumption(k):
    """consumptions of the digit at position e (e = 0 the last digit) over one epoch: counter values
    v in (3, 2*4^k] reached by an increment, digit e goes 1 -> 2 or 3 -> 0."""
    c = [0] * k
    for v in range(4, 2 * 4 ** k + 1):
        for e in range(k):
            d_new = (v // 4 ** e) % 4
            d_old = ((v - 1) // 4 ** e) % 4
            if d_new != d_old and d_old in (1, 3):
                c[e] += 1
    return c


def r_closed(t):
    return 3 + 2 * t - (t + 1) // 2


def sec_F(kmax=7):
    print('F  the epoch theorem')
    ok = all(epoch_consumption(k) == [4 ** k - 1] + [4 ** (k - e) for e in range(1, k)] for k in range(1, 7))
    check(ok, 'one epoch (counter 3 -> 2*4^k, 2*4^k - 3 sweeps): the last digit consumes 4^k - 1 times, the digit '
          'at position e >= 1 consumes 4^(k-e) times (k = 1..6, by counting); digits 1..k-1 together (4^k - 4)/3')
    ok = all(sum(4 ** (k - e) for e in range(1, k)) * 3 == 4 ** k - 4 for k in range(1, 60))
    check(ok, 'sum_{e=1}^{k-1} 4^(k-e) = (4^k - 4)/3 (k = 1..59)')
    ok = all(r_closed(2 * 4 ** k - 4) == 3 * 4 ** k - 3 for k in range(1, 60))
    check(ok, 'r(t) = 3 + 2t - ceil(t/2) gives r(2*4^k - 4) = 3*4^k - 3 at the O2 sweep (k = 1..59)')
    ok = all(3 * (3 * 4 ** k - 6) - (4 ** (k + 1) - 4) == 5 * 4 ** k - 14 > 0 for k in range(1, 200))
    check(ok, '3*(3*4^k - 6) - (4^(k+1) - 4) = 5*4^k - 14 > 0 for every k >= 1 (identity and sign, k = 1..199)')
    ok = all(3 * (3 * 4 ** (k - 1) - 5) - (4 ** k - 4) == 5 * 4 ** (k - 1) - 11 > 0 for k in range(2, 200))
    check(ok, '3*(3*4^(k-1) - 5) - (4^k - 4) = 5*4^(k-1) - 11 > 0 for every k >= 2 (so from the second epoch on '
          'the newborn field is strictly above the minimum)')
    rnd = random.Random(3)
    for k in range(2, kmax + 1):
        for mode in ('zeros', 'random', 'big'):
            fmin = (4 ** k - 4) // 3
            if mode == 'zeros':
                f = [0] * (k - 2) + [fmin]
            elif mode == 'random':
                f = [rnd.randint(0, 3) for _ in range(k - 2)] + [fmin]
            else:
                f = [rnd.randint(0, 50) for _ in range(k - 2)] + [fmin + rnd.randint(0, 20)]
            cfg = (1, [0] * (k - 1) + [3], f, 3)
            N = 2 * 4 ** k - 3
            bad = []
            rem = sum(4 ** (k - e) for e in range(1, k))      # remaining consumption of digits 1..k-1
            m_o2 = None
            for t in range(N):
                W, p, ff, r = cfg
                if t == N - 1:
                    g = o2_by_gaps(cfg)
                    if g is None:
                        bad.append(('O2 outside the gap domain', t))
                        break
                    gap_next, m_o2, _ = g
                V = (3 + t) % 4 ** k
                if p != [(V // 4 ** (k - 1 - i)) % 4 for i in range(k)] or W != (1 if t <= 4 ** k - 4 else 0) \
                        or r != r_closed(t) or ff[-1] < rem:
                    bad.append(('I', t))
                    break
                try:
                    cfg, steps, ev, info = sw.sweep(cfg)
                except sw.Domain as ex:
                    bad.append(('domain', t, str(ex)))
                    break
                rem -= sum(1 for d, _ in info['cons'] if d < k - 1)
                if any(fld == k - 1 for d, fld in info['cons'] if d < k - 1):
                    bad.append(('reached r', t))
                    break
            W, p, ff, r = cfg
            want = (3 * 4 ** k - 5,) if f[-1] > fmin else (3 * 4 ** k - 5, 3 * 4 ** k - 6)
            end_ok = not bad and W == 1 and p == [0] * k + [3] and r == 3 and ff[-1] in want and rem == 0 \
                and cfg == gap_next and (m_o2 == 1) == (ff[-1] == 3 * 4 ** k - 6)
            check(end_ok, f'k = {k}, other fields {mode}, f_(k-1) = {f[-1]} (minimum {fmin}): {N} sweeps, each in the '
                  f'domain with I(t); the O2 sweep has m = {m_o2} block group(s) (m = 1 exactly when -6) and agrees '
                  f'with the gap telescoping; ends at B({k + 1}) with f_k = {ff[-1]} (= 3*4^k - {3 * 4 ** k - ff[-1]}; '
                  f'{"must be -5: start above the minimum" if f[-1] > fmin else "-5 or -6 allowed: start at the minimum"})' +
                  (f'  {bad[:1]}' if bad else ''))


def drec_lines(path):
    """(step, tokens of the tape from its first 1 up to the head cell, head cell excluded) per D-record line.
    The C core writes the head cell into the last token (1.3 = 100 + head cell)."""
    for line in open(path):
        f = line.rstrip('\n').split('\t')
        if len(f) < 4 or not f[0].isdigit():
            continue
        toks = []
        for w in f[3].split():
            ab, _, kk = w.partition('^')
            a, b = map(int, ab.split('.'))
            toks.append([a, b, int(kk) if kk else 1])
        if toks[-1][:2] != [1, 3]:
            yield int(f[0]), None
            continue
        if toks[-1][2] > 1:
            toks[-1][2] -= 1
            toks.append([1, 2, 1])
        else:
            toks[-1] = [1, 2, 1]
        yield int(f[0]), toks


def sec_G(far=True):
    print('G  openings and the real orbit')
    # 274: blank -> D-records up to the second birth (direct simulation)
    t = T.T274
    recs = []
    tape = bytearray(8000)
    pos = 4000
    hi = pos
    q = 'A'
    for n in range(1, 7911):
        a, mv, q = t[q, tape[pos]]
        tape[pos] = a
        pos += mv
        if pos > hi:
            hi = pos
            if q == 'D':
                s = ''.join(map(str, tape[:pos]))
                recs.append((n, T.decode(s[s.find('1'):])))
    first_in = next((r for r in recs if r[1] is not None), None)
    check(first_in == (108, (0, [1], [], 3)), f'274 from blank: the first D-record in the language is at step 108: '
          f'C(no wall; P1; r = 3) ({first_in})')
    b2 = next((r for r in recs if r[0] == 477), None)
    check(b2 == (477, (1, [0, 3], [3], 3)), f'274: the first birth, step 477, is B(2; f1 = 3), one below the minimum '
          f'(4^2 - 4)/3 = 4, so the opening runs on to the next birth ({b2})')
    b3 = recs[-1]
    check(b3 == (7910, (1, [0, 0, 3], [0, 41], 3)), f'274: direct simulation, step 7,910 is a D-record with '
          f'B(3; f = 0, 41) = 0^inf 1100 P0 P0 s^41 P3 s^3 D> ({b3}); 41 >= (4^3 - 4)/3 = 20')
    seq = [r for r in recs if r[0] >= 108]
    cfg, st, path = seq[0][1], 108, []
    for _ in range(len(seq) - 1):
        cfg, n, ev, _ = sw.sweep(cfg)
        st += n
        path.append(st)
    check(path == [r[0] for r in seq[1:]] and [r[1] for r in seq[1:]][-1] == cfg and len(path) == 32,
          f'274: steps 108 -> 7,910 are 32 sweeps of section 4 (every D-record in between matches)')
    # 514: mirrored run of TM514 from blank; records = new leftmost cells in state E (= 274's D under phi)
    t5 = T.table(TM514)
    tape = bytearray(6000)
    pos = 3000
    lo = pos
    q = 'A'
    rec514 = None
    for n in range(1, 2000):
        a, mv, q = t5[q, tape[pos]]
        tape[pos] = a
        pos += mv
        if pos < lo:
            lo = pos
            if q == 'E':
                s = ''.join(map(str, tape[pos + 1:]))[::-1]
                s = s[s.find('1'):]
                d = T.decode(s)
                if d is not None and n >= 792:
                    rec514 = (n, d)
                    break
    check(rec514 == (792, (1, [0, 3], [7], 3)), f'514 from blank (TM514 itself, read mirrored): at step 792 the head '
          f'is in state E (= D of the table) on a new leftmost cell and the tape is B(2; f1 = 7) = '
          f'0^inf 1100 P0 s^7 P3 s^3 D> mirrored ({rec514}); 7 >= 4')
    # same via 274's table from F
    tape = bytearray(6000)
    pos = 3000
    hi = pos
    q = 'F'
    recF = None
    for n in range(1, 2000):
        a, mv, q = t[q, tape[pos]]
        tape[pos] = a
        pos += mv
        if pos > hi:
            hi = pos
            if q == 'D' and n == 792:
                s = ''.join(map(str, tape[:pos]))
                recF = T.decode(s[s.find('1'):])
    check(recF == (1, [0, 3], [7], 3), f'274\'s table from 0^inf F> 0^inf: D-record at step 792 = B(2; 7) ({recF})')
    # abstract orbit vs every real D-record to 1e9
    for which, start, path in (('274', (108, (0, [1], [], 3)), 'work/drec_274_1e9.txt'),
                               ('514', (792, (1, [0, 3], [7], 3)), 'work/drec_514_1e9.txt')):
        st, cfg = start
        n_ok = n_all = 0
        evs = collections.Counter()
        first_bad = None
        it = drec_lines(os.path.join(HERE, path))
        for s_real, tape_s in it:
            if s_real < st:
                continue
            if s_real == st:
                d = T.decode_tokens(tape_s) if tape_s else None
                if d != cfg:
                    first_bad = ('start', s_real)
                    break
                continue
            cfg, n, ev, _ = sw.sweep(cfg)
            st += n
            evs[ev] += 1
            n_all += 1
            d = T.decode_tokens(tape_s) if tape_s else None
            if st == s_real and d == cfg:
                n_ok += 1
            elif first_bad is None:
                first_bad = (s_real, st)
        check(n_ok == n_all and n_all > 18000 and first_bad is None,
              f'{which}: abstract orbit = every real D-record from step {start[0]} to 1e9: {n_ok}/{n_all} sweeps '
              f'{dict(evs)}, configuration and step' + (f' {first_bad}' if first_bad else ''))
    if far:
        for which, start, path in (('274', (108, (0, [1], [], 3)), 'work/drec_274_first_after_1e11.txt'),
                                   ('514', (792, (1, [0, 3], [7], 3)), 'work/drec_514_first_after_1e11.txt')):
            target = next(drec_lines(os.path.join(HERE, path)))
            st, cfg = start
            while st < target[0]:
                cfg, n, ev, _ = sw.sweep(cfg)
                st += n
            check(st == target[0] and T.decode_tokens(target[1]) == cfg,
                  f'{which}: abstract orbit hits the C core\'s first D-record after 1e11 (step {target[0]:,}), '
                  f'configuration and step')


def sec_H(n=3000):
    print('H  agreement with the observe-stage model (informational)')
    sys.path.insert(0, os.path.join(HERE, 'scripts'))
    try:
        import model
        from rules_random import build
    except Exception as ex:
        print('  skip ', ex)
        return
    rnd = random.Random(9)
    c = collections.Counter()
    for _ in range(n):
        cfg = rand_cfg(rnd)
        W, p, f, r = cfg
        try:
            nxt, steps, ev, _ = sw.sweep(cfg)
        except sw.Domain:
            continue
        st = build(W, list(p), list(f), r - 1)
        pred, ev2 = model.sweep(st)
        ff, R = model.fields(pred)
        mp = [d[1] for d in pred['D']]
        mf = ff[1:] if pred['w'] else ff
        same = (pred['w'], mp, mf, R + 1) == nxt and ev2 == ev and model.sweep_steps(st, pred, ev2) == steps
        c[same] += 1
    check(c[False] == 0, f'{c[True]} in-domain random sweeps: section 4 = observe-stage rules R1/R2/PUSH/O1/O2 and '
          f'step rule R3 (configuration, event and steps)' + (f'; {c[False]} differ' if c[False] else ''))


SECTIONS = dict(A=sec_A, B=sec_B, C=sec_C, D=sec_D, E=sec_E, F=sec_F, G=sec_G, H=sec_H)

if __name__ == '__main__':
    want = [a for a in sys.argv[1:] if a in SECTIONS] or list(SECTIONS)
    t0 = time.time()
    for s in want:
        SECTIONS[s]()
    print(f'\n{len(FAIL)} failure(s); {time.time() - t0:.1f} s')
    sys.exit(1 if FAIL else 0)