#!/usr/bin/env python3
"""Standalone independent verifier for the C_3a lower-bound certificate (wins/gsd_3a).

Recomputes EVERYTHING from gsd_3a_cert.json with inline code (stdlib + mpmath only); imports
nothing from the search/harness. Exit 0 iff the certified large-deviation limit lower bound
L_lower is a VALID GHR lower bound on C_3a that STRICTLY beats the previous record.

THE MATH (self-contained).  Capped sparse-digit construction, for depth d = 1,2,3,...:
    U_d = { sum_{i<d} a_i B^i : a_i in A, sum_i a_i <= floor(rho*d) },   B = 2*max(A)+1.
No-carry (B >= 2 max A + 1) makes |U_d+U_d|, |U_d-U_d|, max(U_d) exactly countable. The GHR2007
digit lemma: for any finite U ⊂ Z>=0 with 0 in U and |U-U| < 2 max(U)+1,
    C_3a >= theta(U) := 1 + log(|U-U|/|U+U|) / log(2 max(U)+1).
So C_3a >= theta(U_d) for every d (a finite, checkable bound). As d -> inf:
  * D_d := |U_d - U_d| is SUPERMULTIPLICATIVE (concatenate two feasible difference-words: caps add),
    so (1/d) log D_d -> I_D := lim exists (Fekete), and I_D >= H(pi) for ANY law pi on A-A feasible
    at rho (a feasible pi yields >= exp(d H(pi)) feasible difference-words);
  * S_d := |U_d + U_d| <= #{sum-words y : sum_i y_i <= 2T} (feasible => exists split a with sum a<=T,
    sum(y-a)<=T => sum y <= 2T), so (1/d) log S_d <= I_S_up := 2*rho*s + log sum_y exp(-s*y) for any s;
  * 2 max(U_d)+1 ~ B^d (put max digit on top), so (1/d) log(2 max(U_d)+1) -> log B = log(2 max A+1).
Hence theta(U_d) -> L := 1 + (I_D - I_S)/log(2 max A+1) and, since each theta(U_d) is a valid bound,
    C_3a >= L >= L_lower := 1 + (H(pi) - I_S_up)/log(2 max A+1).
We take pi(delta) ∝ exp(-alpha (p_delta + q_delta)) with (p,q) the minimal digit pair for delta, and
DEFINE rho := E_pi[p] (we are free to pick the cap ratio), making pi feasible by construction. A dual
tilt gives I_D <= 2*rho*alpha + log sum_delta exp(-alpha(p+q)); requiring this < log(2 max A+1) secures
the lemma's d<q hypothesis for large d. All decisions are exact interval-arithmetic comparisons.
"""
import json
import math
import os
import sys
from fractions import Fraction

HERE = os.path.dirname(os.path.abspath(__file__))
# previous record: Mosaic Intelligence 2026 (PR #95), 40-digit certified floor
RECORD_PREV = Fraction(11835129324218615106564747894020784326947, 10**40)
PROVEN_UB = Fraction(4, 3)
DIGIT_LEMMA_CAP = Fraction(5, 4)


def delta_pairs(A):
    Aset = set(A)
    out = []
    for delta in sorted({a - b for a in A for b in A}):
        q = min(b for b in A if b + delta in Aset)
        out.append((q, q + delta))
    return out


# ---- exact finite-d DP counters (independent; brute-force checked at tiny d) ----
def max_u(A, B, d, T):
    A = sorted(set(A)); best = {0: 0}
    for i in range(d):
        w = B**i; nxt = {}
        for sd, v in best.items():
            for a in A:
                ns = sd + a
                if ns <= T and v + a * w > nxt.get(ns, -1):
                    nxt[ns] = v + a * w
        best = nxt
    return max(best.values())


def diff_count(A, B, d, T):
    pairs = delta_pairs(A); states = {(0, 0): 1}
    for _ in range(d):
        nxt = {}
        for (x, y), c in states.items():
            for p, q in pairs:
                nx, ny = x + p, y + q
                if nx <= T and ny <= T:
                    nxt[(nx, ny)] = nxt.get((nx, ny), 0) + c
        states = nxt
    return sum(states.values())


def sum_upper(A, B, d, T):
    A = sorted(set(A)); Aset = set(A)
    Y = sorted({a + b for a in A for b in A})
    minP = {y: min(a for a in A if (y - a) in Aset) for y in Y}
    states = {(0, 0): 1}
    for _ in range(d):
        nxt = {}
        for (u, v), c in states.items():
            for y in Y:
                nu, nv = u + minP[y], v + y
                if nu <= T and nv <= 2 * T:
                    nxt[(nu, nv)] = nxt.get((nu, nv), 0) + c
        states = nxt
    return sum(states.values())


def brute_tiny(A, B, d, T):
    A = sorted(set(A)); vals = [(0, 0)]
    for i in range(d):
        w = B**i
        vals = [(v + a * w, sd + a) for (v, sd) in vals for a in A if sd + a <= T]
    U = sorted({v for v, _ in vals})
    return max(U), len({u + v for u in U for v in U}), len({u - v for u in U for v in U})


def rigorous(A, alpha, s, prec=90):
    import mpmath
    iv = mpmath.iv; iv.prec = prec * 4
    A = sorted(set(A)); m = A[-1]
    pairs = delta_pairs(A)
    al = iv.mpf(str(alpha))
    wts = [iv.exp(-al * iv.mpf(p + q)) for (p, q) in pairs]
    Zd = sum(wts); ps = [w / Zd for w in wts]
    Ep = sum(pi * iv.mpf(p) for pi, (p, q) in zip(ps, pairs))
    Eq = sum(pi * iv.mpf(q) for pi, (p, q) in zip(ps, pairs))
    H = -sum(pi * iv.log(pi) for pi in ps)
    rho = Ep.b
    Y = sorted({a + b for a in A for b in A})
    sv = iv.mpf(str(s))
    Zs = sum(iv.exp(-sv * iv.mpf(y)) for y in Y)
    I_S_up = 2 * iv.mpf(rho) * sv + iv.log(Zs)
    I_D_up = 2 * iv.mpf(rho) * al + iv.log(Zd)
    logq = iv.log(iv.mpf(2 * m + 1))
    L_iv = 1 + (H - I_S_up) / logq
    tol = iv.mpf(10) ** (-prec // 2)
    return {
        "L_iv": L_iv, "L_lower": float(L_iv.a), "rho": float(rho),
        "feasible": bool(Ep.b <= rho + tol and Eq.b <= rho + tol),
        "d_lt_q": bool(I_D_up.b < logq.a),
        "I_D_up": float(I_D_up.b), "H_lo": float(H.a), "I_S_up": float(I_S_up.b),
        "logq_lo": float(logq.a),
    }


def main():
    import mpmath
    cert = json.load(open(os.path.join(HERE, "gsd_3a_cert.json")))
    A = sorted(set(int(x) for x in cert["alphabet"]))
    B = int(cert["base"]); alpha = cert["alpha_tilt"]; s = cert["s_tilt"]
    fail = 0
    print(f"alphabet A (m={A[-1]}, |A|={len(A)}): {A}")
    print(f"base B = {B}  (need >= 2 max A + 1 = {2*A[-1]+1})")
    if A[0] != 0 or len(A) < 2 or B < 2 * A[-1] + 1:
        print("FAIL: admissibility (0 in A, |A|>=2, no-carry base)"); fail += 1

    r = rigorous(A, alpha, s)
    print(f"PASS feasibility (E[p]=E[q]<=rho={r['rho']:.6f})" if r["feasible"]
          else "FAIL: tilted law not feasible"); fail += 0 if r["feasible"] else 1
    print(f"PASS lemma d<q (I_D <= {r['I_D_up']:.5f} < log(2m+1)={r['logq_lo']:.5f})" if r["d_lt_q"]
          else "FAIL: I_D upper bound >= log(2m+1)"); fail += 0 if r["d_lt_q"] else 1

    iv = mpmath.iv
    rec_iv = iv.mpf(str(RECORD_PREV.numerator)) / iv.mpf(str(RECORD_PREV.denominator))
    ub_iv = iv.mpf(4) / iv.mpf(3)
    cap_iv = iv.mpf(5) / iv.mpf(4)
    beats = bool(r["L_iv"] > rec_iv)
    under_ub = bool(r["L_iv"] < ub_iv)
    under_cap = bool(r["L_iv"] < cap_iv)
    print(f"L_lower           = {r['L_lower']:.18f}")
    print(f"previous record   = {float(RECORD_PREV):.18f}  (MI2026, PR #95)")
    print(f"beats record (rigorous interval): {beats}   margin ~ {r['L_lower']-float(RECORD_PREV):+.7f}")
    print(f"under proven UB 4/3: {under_ub}    under digit-lemma cap 5/4: {under_cap}")
    if not beats:
        print("FAIL: does not strictly beat the previous record"); fail += 1
    if not (under_ub and under_cap):
        print("FAIL: violates a known upper bound"); fail += 1

    # finite-d demonstration: theta_lower(U_d) is a genuine climbing GHR value
    bm, bs, bd = brute_tiny(A, B, 2, int(r["rho"] * 2))
    dp_ok = (max_u(A, B, 2, int(r["rho"]*2)) == bm and diff_count(A, B, 2, int(r["rho"]*2)) == bd)
    print("PASS finite-d DP counters == brute force at d=2" if dp_ok else "FAIL: DP != brute")
    fail += 0 if dp_ok else 1
    print("finite-d GHR values theta_lower(U_d) climbing toward L_lower (rigorous <= true theta):")
    prev = -1
    for d in (4, 8, 12, 16):
        T = int(r["rho"] * d)
        D = diff_count(A, B, d, T); S = sum_upper(A, B, d, T); m = max_u(A, B, d, T)
        th = 1 + (math.log(D) - math.log(S)) / math.log(2 * m + 1)
        print(f"   d={d:>2} T={T:>4}: theta_lower={th:.7f} {'up' if th>=prev else 'DOWN'}   [<= {r['L_lower']:.5f}]")
        if th > r["L_lower"] + 1e-6:
            print(f"   FAIL: theta_lower(U_{d}) exceeds L_lower"); fail += 1
        prev = th

    print("RESULT:", "ALL PASS — certified C_3a >= L_lower > previous record" if fail == 0
          else f"{fail} FAILURES")
    return 0 if fail == 0 else 1


if __name__ == "__main__":
    sys.exit(main())