"""Independent exact verifier for "model g" RSC closure certificates (graded closure p2 <= p3 < 1/2, typed
exposures with credits omega2, omega3).  Written separately from rsc-verify.py; it imports no other project
code (not rsc-verify.py, rsc_verify_core.py or rsc_modelg.py) -- only the Python standard library.  The up-set
cone and the rational log bounds come from an earlier independent verifier for the model-e certificates.

usage:  python3 rsc-verify-independent.py cert.json[.gz] [procs]      (also accepts model-e certificates, see below)
        python3 rsc-verify-independent.py --selftest
MYV_SKIP_HEAVY=1 skips the pair and staged checks; it is for negative testing only, and its result is not a
verification.

Model (re-derived from the documented semantics in the docstrings of rsc_modelg.py, not copied):
 tree vertex states  A (sampled), R_c (closure-added, c = #sampled children, capped at 2),
 O00 (not in X, no X child), O01 (no sampled child, >= 1 R child), O1 (one sampled child), O2 (>= 2).
 Weights: sampled qx*y/2, closure-added qx*p_k/2 (needs >= 2 sampled tree-neighbours), outside 2^-m om_x
 (om_x = omega2^c2 omega3^c3, qx = q2^c2 q3^c3), x2 per edge with both ends in X = S u R, x1/2 for an
 undominated O vertex, x q_k for an O vertex with >= 2 sampled neighbours.  Graded rule of the 8-state
 recursion: k counts as ">= 3" (p3, q3) only when the parent is sampled and the vertex is in R2/O2; else p2, q2.
 The merge below is written from these weights; --selftest checks it EXACTLY against a brute-force
 enumeration of all S/R/O assignments of random small trees with exactly this rule, checks that it reduces
 to the uniform (model e) merge at p2 = p3, and checks that it is a lower bound on the "true" graded weight
 (p3, q3 whenever the actual count is >= 3), as the documentation claims.

Checks (all in Fractions; logs by rigorous rational series bounds, no floating point anywhere):
 (P)  0 < p2 <= p3 < 1/2, y > 0, e1 > 0, 0 < e2 < q3.
 (W)  1 <= omega_d <= 2 q_d, and log2(omega2) <= 1 - rho1, log2(omega3) <= (1 - rho1^2)(1 - eps3/ln 2),
      rho1 = 1/(1 + (y/2) q3), eps3 = -ln q3; upper bounds for ln(omega), eps3 and lower bounds for ln 2.
 (K)  K = {x : sum_{s in U} x_s >= 0 for every up-set U of the order R2>R1>R0, R2>O01>O00, O01>O1>O2} (the
      cone generated by covers e_b - e_w and all e_s; asserted); merge(K, R+^8) and merge(R+^8, K) lie in K
      (generators x unit vectors, both sides), final is K-monotone.
 (Z)  all z >= 0, c = min final(z) > 0.
 (a)  every ordered pair z_i o z_j  >=_K  sum w_k z_k with w >= 0, sum w >= 1 (parallel).
 (b)  Init(m)/(e1 e2^m), m = 0..m0, and (y/2, p2/2, 0, ...)/(e1 e2^(m0+1)) dominated likewise.
 (L)  for 1 <= m <= ME: d < dmin3(m) => 2^-m b^-d >= (1 + p3/y) eps3 y / ln(2 q3)   (test3),
      d < dmin2(m) => 2^-m b^-d >= eps2 y / ln 2 (test2), b = 2 + 2 max(p2, p3 - p2)/y; dmin3 <= dmin2.
 (E1) every leafexp entry present is verified (Init_g(m, c2, c3)/(e1 e2^m) dominated, 0 stages).
 (E2) every internal entry present is verified: its stage chain starting from MY Init_g(m, c2, c3)/(e1 e2^m);
      every stage point checked >= 0 (reported; see note), every x o z_b dominated by the next stage (or by Z
      at the last stage) with weight sum >= 1 (parallel); its stage count j = len(combos), len(stages) = j - 1.
 (H)  coverage, recomputed from scratch: for every m <= ME and every exposure the model allows, the point
      (qx, om_x) must lie in the exact convex hull of points that were verified with at most the allowed
      number of stages:  internal owner failing test3: all (0, c3), 0 <= c3 <= m, from points with
      j <= max(1, dmin3 - 1); failing test2: all c2 + c3 <= m from points with j <= max(1, dmin2 - 1);
      leaves (d = 1): if dmin3 <= 1 all (0, c3), c3 <= m - 1, if dmin2 <= 1 all c2 + c3 <= m - 1, from
      0-stage points.  (0, 0) is the plain init (inits for m <= m0, the tail point for m > m0).  The hull
      test is exact (monotone chain on Fractions; degenerate 1- and 2-point hulls handled), and every
      lattice point is tested, not only hull vertices.  The lattice hull vertices are also listed and
      compared with what the certificate provides.
 (T)  q3 > e2 and (q3/e2)^(ME+1) >= e2^-(m0+1): for m > ME every Init_g(m, c2, c3)/(e1 e2^m) dominates the
      tail point (qx >= q3^m since q3 <= q2; the O-part is >= 0).
Model-e certificates (model='e', keys p, omega, dmin, leafexp "m,c", internal 'c') are mapped to
p2 = p3 = p, omega2 = omega3 = omega, dmin2 = dmin3 = dmin, exposures (0, c).
Note on (E2) nonnegativity: with true messages M >= 0 the step merge(M,M') - merge(x,z) =
merge(M - x, z) + merge(M, M' - z) needs only z >= 0 and M >= 0, so stage-point signs are not logically
required; they are checked and reported anyway."""
import gzip, itertools, json, os, random, sys, time
from fractions import Fraction as Fr
from multiprocessing import Pool

ST = ['A', 'R0', 'R1', 'R2', 'O00', 'O01', 'O1', 'O2']
IX = {s: i for i, s in enumerate(ST)}
COV = [('R2', 'R1'), ('R1', 'R0'), ('R2', 'O01'), ('O01', 'O00'), ('O01', 'O1'), ('O1', 'O2')]

def _closure(rel):
    rel = set(rel)
    while True:
        new = {(a, d) for (a, b) in rel for (c, d) in rel if b == c} - rel
        if not new: return rel
        rel |= new
GE = _closure(COV)
UPS = []
for r in range(1, 9):
    for S in itertools.combinations(range(8), r):
        Ss = {ST[i] for i in S}
        if all(b in Ss for (b, w) in GE if w in Ss): UPS.append(S)

def inK(D): return all(sum(D[i] for i in U) >= 0 for U in UPS)

# ---------------------------------------------------------------- graded merge (re-derived)
def merge(P, Q, pr):
    """P: partial message of parent v, Q: completed message of a new child u.  pr = (p2, p3)."""
    p2, p3 = pr; q2, q3 = 1 - p2, 1 - p3
    A, R0, R1, R2, O00, O01, O1, O2 = P
    a, r0, r1, r2, o00, o01, o1, o2 = Q
    # v sampled: u sampled (edge X-X: 2); u in R with 1 sampled child -> k=2 (p2, carried); u in R2 -> k>=3:
    # replace the carried p2 by p3 (factor p3/p2), edge X-X: 2; u in R0 invalid (k=1);
    # u outside: dominated by v; O1 -> k=2 not added q2, O2 -> k>=3 not added q3.
    toS = 2 * a + 2 * r1 + 2 * (p3 / p2) * r2 + o00 + o01 + q2 * o1 + q3 * o2
    # v in R: u sampled (2, increments v's count); u in R needs its own 2 sampled children (k=2: p2), edge 2;
    # u outside: dominated by v, O2 has k = 2 counted (parent not sampled): q2.
    toR_s = 2 * a
    toR_u = 2 * r2 + o00 + o01 + o1 + q2 * o2
    # v outside: u sampled (increments v's sampled count), u in R2 (gives v an R child), u outside:
    # O00 child is undominated (1/2), O2 -> q2.
    toO_s = a
    toO_r = r2
    toO_u = o00 / 2 + o01 + o1 + q2 * o2
    return [A * toS,
            R0 * toR_u,
            R1 * toR_u + R0 * toR_s,
            R2 * (toR_u + toR_s) + R1 * toR_s,
            O00 * toO_u,
            O01 * (toO_u + toO_r) + O00 * toO_r,
            O1 * (toO_u + toO_r) + (O00 + O01) * toO_s,
            O2 * (toO_u + toO_r + toO_s) + O1 * toO_s]

def fin(P, pr):
    """root: R0, R1 invalid (k < 2), O00 undominated (1/2), O2 a candidate with k = 2 counted (q2)."""
    return P[0] + P[3] + P[4] / 2 + P[5] + P[6] + (1 - pr[0]) * P[7]

def init_g(y, pr, m, c2, c3, om2, om3):
    qx = (1 - pr[0]) ** c2 * (1 - pr[1]) ** c3
    return [qx * y / 2, qx * pr[0] / 2, Fr(0), Fr(0), Fr(1, 2 ** m) * om2 ** c2 * om3 ** c3, Fr(0), Fr(0), Fr(0)]

# ---------------------------------------------------------------- domination
def _w(x, mode):
    if isinstance(x, list): return Fr(x[0], x[1])
    return Fr(x).limit_denominator(10 ** 12) if mode == 0 else Fr(x)

def dom(P, lst, pts):
    """P >=_K sum w_k pts_k with w >= 0 (non-positive weights dropped), sum w >= 1.  Two exact decodings of the
    float weights are tried (limit_denominator(1e12), then the exact binary value); either suffices."""
    for mode in (0, 1):
        ws = [(k, _w(w, mode)) for k, w in lst]
        ws = [(k, w) for k, w in ws if w > 0]
        if sum(w for _, w in ws) < 1: continue
        D = list(P)
        for k, w in ws:
            z = pts[k]
            for d in range(8):
                if z[d]: D[d] -= w * z[d]
        if inK(D): return True
    return False

# ---------------------------------------------------------------- rigorous logs
def nlog1m(x, K):
    """(lo, hi) for -ln(1 - x), 0 <= x < 1, x rational."""
    x = Fr(x); s = Fr(0); t = Fr(1)
    for k in range(1, K + 1):
        t *= x; s += t / k
    return s, s + t * x / ((K + 1) * (1 - x))

LN2 = nlog1m(Fr(1, 2), 240)

def ln_hi(w):
    """upper bound for ln(w), w >= 1: ln w = -ln(1 - (1 - 1/w))."""
    return nlog1m(1 - 1 / Fr(w), 200)[1]

# ---------------------------------------------------------------- exact 2-d convex hull
def cross(o, a, b): return (a[0] - o[0]) * (b[1] - o[1]) - (a[1] - o[1]) * (b[0] - o[0])

def hull(pts):
    P = sorted(set(pts))
    if len(P) <= 2: return P
    lo, up = [], []
    for p in P:
        while len(lo) >= 2 and cross(lo[-2], lo[-1], p) <= 0: lo.pop()
        lo.append(p)
    for p in reversed(P):
        while len(up) >= 2 and cross(up[-2], up[-1], p) <= 0: up.pop()
        up.append(p)
    H = lo[:-1] + up[:-1]
    return H if len(H) >= 1 else P

def in_hull(H, x):
    """x in conv(H), H from hull() (CCW, no collinear vertices); exact, degenerate cases included."""
    if not H: return False
    if len(H) == 1: return x == H[0]
    if len(H) == 2:
        a, b = H
        return cross(a, b, x) == 0 and min(a[0], b[0]) <= x[0] <= max(a[0], b[0]) and min(a[1], b[1]) <= x[1] <= max(a[1], b[1])
    return all(cross(H[i], H[(i + 1) % len(H)], x) >= 0 for i in range(len(H)))

# ---------------------------------------------------------------- workers
G = {}
def _winit(Z, pr): G['Z'] = Z; G['pr'] = pr

def _pairs(args):
    keys, combos = args
    Z, pr = G['Z'], G['pr']
    for (i, j), cb in zip(keys, combos):
        if not dom(merge(Z[i], Z[j], pr), cb, Z): return (i, j)
    return None

def _stage(args):
    tag, prev, rows, combos, target = args
    Z, pr = G['Z'], G['pr']
    tgt = Z if target is None else target
    for (a, b), cb in zip(rows, combos):
        if cb is None or not dom(merge(prev[a], Z[b], pr), cb, tgt): return (tag, a, b)
    return None

# ---------------------------------------------------------------- self tests of the merge
def _rand_tree(rng, s):
    par = {0: None}
    for v in range(1, s): par[v] = rng.randrange(v)
    return par

def _brute(par, m, c2, c3, y, pr, om2, om3, rule):
    """sum over all S/R/O assignments of the documented weights.  rule='8': graded rule of the 8 states,
    rule='true': p3/q3 whenever the number of sampled neighbours is >= 3."""
    p2, p3 = pr; q2, q3 = 1 - p2, 1 - p3
    V = list(par); ch = {v: [u for u in V if par[u] == v] for v in V}
    nb = {v: ch[v] + ([par[v]] if par[v] is not None else []) for v in V}
    tot = Fr(0)
    for asg in itertools.product('SRO', repeat=len(V)):
        st = dict(zip(V, asg)); w = Fr(1); ok = True
        for v in V:
            qx = q2 ** c2[v] * q3 ** c3[v]
            k = sum(st[u] == 'S' for u in nb[v])
            kc = sum(st[u] == 'S' for u in ch[v]); ps = par[v] is not None and st[par[v]] == 'S'
            big = (kc >= 2 and ps) if rule == '8' else (k >= 3)
            if st[v] == 'S': w *= qx * y / 2
            elif st[v] == 'R':
                if k < 2: ok = False; break
                w *= qx * (p3 if big else p2) / 2
            else:
                w *= Fr(1, 2 ** m[v]) * om2 ** c2[v] * om3 ** c3[v]
                if not any(st[u] in 'SR' for u in nb[v]): w /= 2
                if k >= 2: w *= (q3 if big else q2)
        if not ok: continue
        for v in V:
            if par[v] is not None and st[v] in 'SR' and st[par[v]] in 'SR': w *= 2
        tot += w
    return tot

def _msg(par, v, m, c2, c3, y, pr, om2, om3):
    P = init_g(y, pr, m[v], c2[v], c3[v], om2, om3)
    for u in [u for u in par if par[u] == v]:
        P = merge(P, _msg(par, u, m, c2, c3, y, pr, om2, om3), pr)
    return P

def _merge_uniform(P, Qm, q):   # the model-e merge of the earlier independent verifier
    A, R0, R1, R2, O00, O01, O1, O2 = P
    a, r0, r1, r2, o00, o01, o1, o2 = Qm
    toS = 2 * a + 2 * r1 + 2 * r2 + o00 + o01 + q * o1 + q * o2
    toR_s, toR_u = 2 * a, 2 * r2 + o00 + o01 + o1 + q * o2
    toO_s, toO_r, toO_u = a, r2, o00 / 2 + o01 + o1 + q * o2
    return [A * toS, R0 * toR_u, R1 * toR_u + R0 * toR_s, R2 * (toR_u + toR_s) + R1 * toR_s, O00 * toO_u,
            O00 * toO_r + O01 * (toO_u + toO_r), (O00 + O01) * toO_s + O1 * (toO_u + toO_r), O1 * toO_s + O2 * (toO_u + toO_r + toO_s)]

def selftest(ntrees=250, seed=5, pr0=None, y0=None, om0=None):
    rng = random.Random(seed)
    worst_lb = None; nstrict = ndiff = 0
    for it in range(ntrees):
        s = rng.randint(1, 6)
        par = _rand_tree(rng, s)
        m = {v: rng.randint(0, 3) for v in par}
        c2 = {v: rng.randint(0, m[v]) for v in par}
        c3 = {v: rng.randint(0, m[v] - c2[v]) for v in par}
        if pr0 and it % 2 == 0: pr, y, (om2, om3) = pr0, y0, om0
        else:
            p2 = Fr(rng.randint(1, 30), 100); p3 = min(Fr(49, 100), p2 + Fr(rng.randint(0, 20), 100))
            pr, y = (p2, p3), Fr(rng.randint(10, 90), 100)
            om2 = 1 + Fr(rng.randint(0, 99), 100) * (1 - 2 * p2); om3 = 1 + Fr(rng.randint(0, 99), 100) * (1 - 2 * p3)
        a = fin(_msg(par, 0, m, c2, c3, y, pr, om2, om3), pr)
        b = _brute(par, m, c2, c3, y, pr, om2, om3, '8')
        assert a == b, ("merge != brute force", par, m, c2, c3, y, pr)
        t = _brute(par, m, c2, c3, y, pr, om2, om3, 'true')
        assert t >= b, ("8-state rule is not a lower bound", par, m, y, pr)
        if t > b: nstrict += 1
        if b != _brute(par, m, c2, c3, y, (pr[0], pr[0]), om2, om3, '8'): ndiff += 1
        r = t / b if b else None
        if r is not None and (worst_lb is None or r < worst_lb): worst_lb = r
    # reduction to the uniform merge at p2 = p3
    for it in range(300):
        p = Fr(rng.randint(1, 49), 100)
        P = [Fr(rng.randint(0, 9), rng.randint(1, 9)) for _ in range(8)]
        Q = [Fr(rng.randint(0, 9), rng.randint(1, 9)) for _ in range(8)]
        assert merge(P, Q, (p, p)) == _merge_uniform(P, Q, 1 - p)
    assert nstrict > 0 and ndiff > 0   # the graded rule is actually exercised
    return dict(trees=ntrees, min_true_over_8state=float(worst_lb), trees_true_gt_8state=nstrict, trees_graded_ne_uniform=ndiff)

# ---------------------------------------------------------------- main verification
def F(v): return Fr(v[0], v[1])

def load(path):
    op = gzip.open if path.endswith('.gz') else open
    with op(path, 'rt') as f: return json.load(f)

def verify(cert, procs=4, log=print):
    t0 = time.time()
    def L(*a): log(f"[{time.time() - t0:8.1f}s]", *a)
    model = cert.get('model')
    if model == 'g':
        y, p2, p3 = F(cert['y']), F(cert['p2']), F(cert['p3'])
        om2, om3 = F(cert['omega2']), F(cert['omega3'])
        dmin3 = {int(k): v for k, v in cert['dmin3'].items()}; dmin2 = {int(k): v for k, v in cert['dmin2'].items()}
        leafent = {}
        for k, v in cert['leafexp'].items():
            m, a, b = map(int, k.split(',')); leafent[(m, a, b)] = v
        internal = [(ex['m'], ex['c2'], ex['c3'], ex) for ex in cert['internal']]
    elif model == 'e':
        y, p2 = F(cert['y']), F(cert['p']); p3 = p2
        om2 = om3 = F(cert['omega'])
        dmin3 = {int(k): v for k, v in cert['dmin'].items()}; dmin2 = dict(dmin3)
        leafent = {}
        for k, v in cert['leafexp'].items():
            m, a = map(int, k.split(',')); leafent[(m, 0, a)] = v
        internal = [(ex['m'], 0, ex['c'], ex) for ex in cert['internal']]
    else:
        return False, "unknown model", {}
    e1, e2 = F(cert['e1']), F(cert['e2']); pr = (p2, p3); q2, q3 = 1 - p2, 1 - p3
    m0, ME = cert['m0'], cert['ME']
    info = dict(model=model, y=y, p2=p2, p3=p3, omega2=om2, omega3=om3, e1=e1, e2=e2)
    # (P)
    if not (0 < p2 <= p3 < Fr(1, 2) and y > 0 and e1 > 0 and 0 < e2 < q3 < 1): return False, "(P) parameter ranges", info
    # (W)
    eps2, eps3 = nlog1m(p2, 200), nlog1m(p3, 200)
    rho1 = 1 / (1 + y / 2 * q3)
    if not (1 <= om2 <= 2 * q2 and 1 <= om3 <= 2 * q3): return False, "(W) 1 <= omega_d <= 2 q_d", info
    w2 = ln_hi(om2) <= (1 - rho1) * LN2[0]
    w3 = ln_hi(om3) <= (1 - rho1 ** 2) * (LN2[0] - eps3[1])
    if not w2: return False, "(W) omega2: log2(omega2) <= 1 - rho1 fails (rigorous bounds)", info
    if not w3: return False, "(W) omega3: log2(omega3) <= (1-rho1^2)(1-eps3/ln2) fails (rigorous bounds)", info
    L("(W) omega admissible (rigorous rational log bounds)")
    # (K)
    E = [[Fr(int(i == j)) for j in range(8)] for i in range(8)]
    gens = [[E[IX[b]][d] - E[IX[w]][d] for d in range(8)] for b, w in COV] + E
    assert all(inK(g) for g in gens)
    # K = cone(gens): its dual = nonnegative monotone functions = cone of up-set indicators (layer cake)
    for g in gens:
        if fin(g, pr) < 0: return False, "(K) final not K-monotone", info
        for e in E:
            if not inK(merge(g, e, pr)) or not inK(merge(e, g, pr)): return False, "(K) cone invariance", info
    L(f"(K) cone invariant under graded merge, final monotone ({len(UPS)} up-sets)")
    # (Z)
    Z = [[F(v) for v in z] for z in cert['Z']]; n = len(Z)
    if not all(len(z) == 8 and all(v >= 0 for v in z) for z in Z): return False, "(Z) negative point", info
    c = min(fin(z, pr) for z in Z)
    if not c > 0: return False, "(Z) c <= 0", info
    info.update(c=c, npts=n)
    L(f"(Z) {n} points >= 0, c = {float(c):.8f}")
    def chk(P, key):
        return key in cert['inits'] and dom(P, cert['inits'][key], Z)
    # (b)
    for m in range(m0 + 1):
        if not chk([x / (e1 * e2 ** m) for x in init_g(y, pr, m, 0, 0, om2, om3)], str(m)): return False, f"(b) init m={m}", info
    uinf = [x / (e1 * e2 ** (m0 + 1)) for x in [y / 2, p2 / 2, 0, 0, 0, 0, 0, 0]]
    if not chk(uinf, 'inf'): return False, "(b) tail point", info
    L(f"(b) inits m=0..{m0} and tail point OK")
    # (T)
    if not (q3 > e2 and (q3 / e2) ** (ME + 1) >= 1 / e2 ** (m0 + 1)): return False, "(T) tail", info
    L(f"(T) tail OK (ME={ME}, m0={m0})")
    # (L)
    if sorted(dmin3) != list(range(1, ME + 1)) or sorted(dmin2) != list(range(1, ME + 1)): return False, "(L) dmin keys", info
    b = 2 + 2 * max(p2, p3 - p2) / y
    rhs3 = (1 + p3 / y) * eps3[1] * y / (LN2[0] - eps3[1])
    rhs2 = eps2[1] * y / LN2[0]
    nl = 0
    for m in range(1, ME + 1):
        if not (0 <= dmin3[m] <= dmin2[m]): return False, f"(L) dmin3 <= dmin2 at m={m}", info
        for d in range(1, dmin3[m]):
            if not Fr(1, 2 ** m) / b ** d >= rhs3: return False, f"(L) test3 m={m} d={d}", info
            nl += 1
        for d in range(1, dmin2[m]):
            if not Fr(1, 2 ** m) / b ** d >= rhs2: return False, f"(L) test2 m={m} d={d}", info
            nl += 1
    L(f"(L) classification OK ({nl} light (m,d) tests); dmin3={[dmin3[m] for m in range(1, 9)]}.. dmin2={[dmin2[m] for m in range(1, 9)]}..")
    # (E1) leaf entries
    leafok = set()
    for (m, a2, a3), cb in leafent.items():
        if not (m >= 1 and a2 >= 0 and a3 >= 0): continue
        if dom([x / (e1 * e2 ** m) for x in init_g(y, pr, m, a2, a3, om2, om3)], cb, Z): leafok.add((m, a2, a3))
        else: return False, f"(E1) leaf exposure m={m} c=({a2},{a3}) not dominated", info
    L(f"(E1) {len(leafok)} leaf exposure entries verified")
    # internal entries: parse (stage chains start from MY Init_g); verified below, after the coverage test
    ents = []
    nneg = 0
    for idx, (m, a2, a3, ex) in enumerate(internal):
        cmbs = ex['combos']; j = len(cmbs)
        if not (m >= 1 and a2 >= 0 and a3 >= 0 and j >= 1 and len(ex['stages']) == j - 1 and ex.get('j', j) == j):
            return False, f"(E2) malformed internal entry #{idx}", info
        Y0 = [[x / (e1 * e2 ** m) for x in init_g(y, pr, m, a2, a3, om2, om3)]]
        stages = [Y0] + [[[F(v) for v in z] for z in Yl] for Yl in ex['stages']]
        for Yl in stages:
            for z in Yl:
                if len(z) != 8: return False, f"(E2) malformed stage point entry #{idx}", info
                if any(v < 0 for v in z): nneg += 1
        ents.append((m, a2, a3, j, stages, cmbs))
    L(f"(E2) {len(ents)} internal entries parsed; stage points with a negative coordinate: {nneg}")
    info['neg_stage_points'] = nneg
    # (H) coverage (uses the entries' stage counts; every entry is verified below, any failure -> FAILED)
    intok = {}
    for (m, a2, a3, j, _, _) in ents:
        intok[(m, a2, a3)] = min(j, intok.get((m, a2, a3), j))
    def pt(a2, a3): return (q2 ** a2 * q3 ** a3, om2 ** a2 * om3 ** a3)
    nreq = 0; nvert = 0; missing_vertices = []
    for m in range(1, ME + 1):
        base = [pt(0, 0)]                      # plain init: inits (m <= m0) or tail point (m > m0)
        leafpts = [pt(a2, a3) for (mm, a2, a3) in leafok if mm == m]
        def avail(J): return base + leafpts + [pt(a2, a3) for (mm, a2, a3), jj in intok.items() if mm == m and jj <= J]
        reqs = []
        j3, j2 = max(1, dmin3[m] - 1), max(1, dmin2[m] - 1)
        reqs.append(('internal/test3', [(0, a3) for a3 in range(m + 1)], avail(j3)))
        reqs.append(('internal/test2', [(a2, a3) for a2 in range(m + 1) for a3 in range(m + 1 - a2)], avail(j2)))
        if dmin3[m] <= 1: reqs.append(('leaf/test3', [(0, a3) for a3 in range(m)], base + leafpts))
        if dmin2[m] <= 1: reqs.append(('leaf/test2', [(a2, a3) for a2 in range(m) for a3 in range(m - a2)], base + leafpts))
        for name, lat, av in reqs:
            H = hull(av)
            for (a2, a3) in lat:
                nreq += 1
                if not in_hull(H, pt(a2, a3)):
                    return False, f"(H) coverage {name} m={m} exposure ({a2},{a3}) not in hull of verified points", info
            # informational: hull vertices of the lattice vs points provided
            LV = hull([pt(a2, a3) for a2, a3 in lat])
            nvert += len(LV)
            avs = set(av)
            missing_vertices += [(name, m) for v in LV if v not in avs]
    L(f"(H) coverage OK: {nreq} (lattice point, owner class) pairs inside exact hulls; {nvert} lattice-hull vertices, "
      f"{len(missing_vertices)} not provided verbatim (covered anyway)")
    for m in range(1, min(ME, 8) + 1):
        lat = [(a2, a3) for a2 in range(m + 1) for a3 in range(m + 1 - a2)]
        LV = [v for v in hull([pt(a2, a3) for a2, a3 in lat])]
        nmix = sum(1 for (a2, a3) in lat if a2 >= 1 and pt(a2, a3) in LV)
        L(f"    m={m}: lattice c2+c3<=m has {len(LV)} hull vertices ({nmix} with c2>=1); cert internal entries at m: "
          f"{sorted((a2, a3, jj) for (mm, a2, a3), jj in intok.items() if mm == m)}")
    if os.environ.get('MYV_SKIP_HEAVY'):
        L("MYV_SKIP_HEAVY set: pair and staged checks SKIPPED (negative-testing mode, not a verification)")
        return None, "heavy checks skipped", info
    # (a) + (E2) in parallel
    pool = Pool(procs, initializer=_winit, initargs=(Z, pr))
    try:
        keys = [(i, j) for i in range(n) for j in range(n)]
        tasks = []
        B = 1500
        for s in range(0, len(keys), B):
            ks = keys[s:s + B]
            tasks.append((ks, [cert['combos'].get(f"{i},{j}", []) for i, j in ks]))
        done = 0
        for r in pool.imap_unordered(_pairs, tasks):
            if r is not None: return False, f"(a) pair {r}", info
            done += 1
            if done % max(1, len(tasks) // 10) == 0: L(f"(a) pairs {done}/{len(tasks)} chunks")
        L(f"(a) all {len(keys)} pairs OK")
        etasks = []
        for idx, (m, a2, a3, j, stages, cmbs) in enumerate(ents):
            for i in range(1, j + 1):
                prev = stages[i - 1]; target = stages[i] if i < j else None
                cb = cmbs[i - 1]
                rows = [(a, bb) for a in range(len(prev)) for bb in range(n)]
                for s in range(0, len(rows), 1500):
                    rr = rows[s:s + 1500]
                    etasks.append(((idx, m, a2, a3, i), prev, rr, [cb.get(f"{a},{bb}") for a, bb in rr], target))
        L(f"(E2) {sum(len(t[2]) for t in etasks)} staged products to check")
        done = 0
        for r in pool.imap_unordered(_stage, etasks, chunksize=1):
            if r is not None:
                (idx, m, a2, a3, i), a, bb = r
                return False, f"(E2) internal exposure m={m} c=({a2},{a3}) stage {i} row ({a},{bb})", info
            done += 1
            if done % max(1, len(etasks) // 10) == 0: L(f"(E2) {done}/{len(etasks)} chunks")
        L("(E2) all internal staged products OK")
    finally:
        pool.terminate(); pool.join()
    # isolated-root remark: an F consisting of one vertex with dmin3(m) = 0 would fail at d = 0
    zero = [m for m in range(1, ME + 1) if dmin3[m] == 0]
    if zero: L(f"remark: dmin3(m) = 0 for m in {zero[0]}..{zero[-1]} (single-vertex F would fail at d=0; leaf rule c <= m-1 used, model-level)")
    return True, "all conditions", info

if __name__ == '__main__':
    if sys.argv[1] == '--selftest':
        print("selftest: merge == brute force (exact), lower bound vs true graded, p2=p3 reduction:", selftest())
        sys.exit(0)
    path = sys.argv[1]; procs = int(sys.argv[2]) if len(sys.argv) > 2 else 4
    t = time.time()
    cert = load(path)
    print(f"loaded {path} in {time.time() - t:.1f}s", flush=True)
    pr0 = None
    try:
        if cert.get('model') == 'g': pr0, y0, om0 = (F(cert['p2']), F(cert['p3'])), F(cert['y']), (F(cert['omega2']), F(cert['omega3']))
        else: pr0, y0, om0 = (F(cert['p']), F(cert['p'])), F(cert['y']), (F(cert['omega']), F(cert['omega']))
        print("selftest (trees, min true/8-state ratio):", selftest(200, 3, pr0, y0, om0), flush=True)
    except AssertionError as e:
        print("FAILED selftest", e); sys.exit(1)
    ok, why, info = verify(cert, procs, log=lambda *a: print(*a, flush=True))
    out = {k: (float(v) if isinstance(v, Fr) else v) for k, v in info.items()}
    print(path, "VERIFIED" if ok else ("PARTIAL (not verified): " + why if ok is None else "FAILED: " + why), out, flush=True)
    if ok: print("exact values:", {k: f"{v.numerator}/{v.denominator}" for k, v in info.items() if isinstance(v, Fr)})
    sys.exit(0 if ok else 1)