"""[Ancillary helper module, imported by rsc-verify.py.]"""
"""Model "g" of RSC: model e with graded closure and typed exposures.
A candidate with exactly two sampled F-neighbours joins R w.p. p2, one with >= 3 w.p. p3 (p2 <= p3 < 1/2).
The 8-state recursion knows the count >= 3 only for R2/O2 children of a sampled parent and uses
(p2, q2) elsewhere (a lower bound on W, since p_k is nondecreasing and f(in) >= f(out)).
Exposed outside vertices carry a type d = number of S-neighbours:
  type 2 costs q2 = 1 - p2 and gets the domination credit omega2 <= 2^(1 - kappa1), kappa1 = rho1;
  type 3 (d >= 3) costs q3 = 1 - p3 and gets omega3 <= 2^(1 - kappa2), kappa2 = rho1^2 + (eps3/ln 2)(1 - rho1^2),
  rho1 = (1 + (y/2) q3)^-1, eps3 = -ln q3; also omega_d <= 2 q_d.
Lightness tests (b = 2 + 2 max(p2, p3 - p2)/y, graded Lemma 5):
  test3: 2^-m b^-d >= (1 + p3/y) eps3 y / ln(2 q3)   (owners passing it expose nothing);
  test2: 2^-m b^-d >= eps2 y / ln 2                    (owners passing it expose only type-3 vertices).
A failing owner with m owned outside vertices exposes c2 type-2 and c3 type-3 vertices, c2 + c3 <= m
(leaves: <= m - 1), c2 >= 1 only if it fails test2.
init_g(m, c2, c3) = (qx y/2, qx p2/2, 0, 0, 2^-m omega2^c2 omega3^c3, 0, 0, 0), qx = q2^c2 q3^c3.
Since the final exposure point is linear in init_g and the antinorm region is convex, only the vertices of
the convex hull of {(qx, omx)} over the exposure lattice have to be checked (expo_vertices)."""
import math
from fractions import Fraction as Fr
import numpy as np

def kappas(p2, p3, y):
    q3 = 1 - p3; e3 = -math.log(q3)
    r1 = 1 / (1 + y / 2 * q3)
    return r1, r1 * r1 + e3 / math.log(2) * (1 - r1 * r1)

def omega_default(p2, p3, y):
    k1, k2 = kappas(p2, p3, y)
    w2 = min(2 ** (1 - k1), 2 * (1 - p2)); w3 = min(2 ** (1 - k2), 2 * (1 - p3))
    return math.floor(w2 * 1e6) / 1e6, math.floor(w3 * 1e6) / 1e6

def thresholds(p2, p3, y):
    """(thr3, thr2, lb): (m, d) fails test3 iff m + d lb > thr3, fails test2 iff m + d lb > thr2."""
    q2, q3 = 1 - p2, 1 - p3; e2, e3 = -math.log(q2), -math.log(q3)
    lb = math.log2(2 + 2 * max(p2, p3 - p2) / y)
    thr3 = math.log2(math.log(2 * q3) / (e3 * y)) - math.log2(1 + p3 / y)
    thr2 = math.log2(math.log(2) / (e2 * y))
    return thr3, thr2, lb

def dmin_of(m, thr, lb):
    d = 0
    while m + d * lb <= thr + 1e-12: d += 1
    return d

def jmin_of(m, thr, lb):
    return max(1, dmin_of(m, thr, lb) - 1)

def leaf_fails(m, thr, lb):
    return m >= 1 and m + lb > thr + 1e-12

def init_g(y, p2, p3, m, c2, c3, om2, om3):
    qx = (1 - p2) ** c2 * (1 - p3) ** c3
    return np.array([qx * y / 2, qx * p2 / 2, 0, 0, 2.0 ** (-m) * om2 ** c2 * om3 ** c3, 0, 0, 0])

def _hull(pts):
    """vertices of the convex hull of 2-d points (exact if the coordinates are Fractions); monotone chain."""
    P = sorted(set(pts))
    if len(P) <= 2: return P
    def cross(o, a, b): return (a[0] - o[0]) * (b[1] - o[1]) - (a[1] - o[1]) * (b[0] - o[0])
    lo = []
    for p in P:
        while len(lo) >= 2 and cross(lo[-2], lo[-1], p) <= 0: lo.pop()
        lo.append(p)
    up = []
    for p in reversed(P):
        while len(up) >= 2 and cross(up[-2], up[-1], p) <= 0: up.pop()
        up.append(p)
    return lo[:-1] + up[:-1]

def expo_vertices(cmax, q2, q3, om2, om3, exact=False):
    """(V3, Vmix): exposure counts to check.  V3 = all (0, c3), 1 <= c3 <= cmax (type-3-only exposures, allowed at
    owners failing test3).  Vmix = hull vertices with c2 >= 1 of {(q2^c2 q3^c3, om2^c2 om3^c3): c2 + c3 <= cmax}
    (owners failing test2); every other lattice point is a convex combination of Vmix and of (0, c3) points."""
    if cmax <= 0: return [], []
    V3 = [(0, c3) for c3 in range(1, cmax + 1)]
    if exact:
        q2, q3, om2, om3 = Fr(q2), Fr(q3), Fr(om2), Fr(om3)
    pts = {}
    for c2 in range(cmax + 1):
        for c3 in range(cmax + 1 - c2):
            pts[(q2 ** c2 * q3 ** c3, om2 ** c2 * om3 ** c3)] = (c2, c3)
    H = _hull(list(pts.keys()))
    Vmix = sorted(pts[h] for h in H if pts[h][0] >= 1)
    return V3, Vmix

if __name__ == '__main__':
    import sys
    y, p2, p3 = 0.31, 0.07, 0.077
    if len(sys.argv) > 3: y, p2, p3 = [float(a) for a in sys.argv[1:4]]
    om2, om3 = omega_default(p2, p3, y)
    thr3, thr2, lb = thresholds(p2, p3, y)
    print(f"y={y} p2={p2} p3={p3}: omega2={om2} omega3={om3} thr3={thr3:.4f} thr2={thr2:.4f} lb={lb:.4f}")
    print("jmin3:", [jmin_of(m, thr3, lb) for m in range(1, 9)], "jmin2:", [jmin_of(m, thr2, lb) for m in range(1, 9)])
    print("leaf fails3 from m =", min(m for m in range(1, 40) if leaf_fails(m, thr3, lb)),
          "fails2 from m =", min(m for m in range(1, 40) if leaf_fails(m, thr2, lb)))
    tot3 = totm = 0
    for m in range(1, 33):
        V3, Vm = expo_vertices(m, 1 - p2, 1 - p3, om2, om3)
        V3x, Vmx = expo_vertices(m, Fr(1 - p2).limit_denominator(10**6), Fr(1 - p3).limit_denominator(10**6),
                                 Fr(om2).limit_denominator(10**6), Fr(om3).limit_denominator(10**6), exact=True)
        tot3 += len(V3); totm += len(Vm)
        if m <= 8 or m == 32: print(f"m={m}: |V3|={len(V3)} |Vmix|={len(Vm)} (exact {len(Vmx)}) Vmix={Vm if m <= 6 else '...'}")
    print("totals over m<=32:", tot3, totm, "lattice size", sum((m + 1) * (m + 2) // 2 - (m + 1) for m in range(1, 33)))