"""Core machinery for the exact rational SOS-certificate task.

Problem: given an integer-coefficient polynomial p in n = 5 variables of
degree 6, certify a global lower bound
        p(x) >= gamma   for all x in R^5
by exhibiting an EXACT rational sum-of-squares certificate

        p(x) - gamma  =  z(x)^T Q z(x),      Q symmetric, Q positive definite

where z(x) is the vector of all 56 monomials of degree <= 3.  The submission
gives integers: a denominator `den` > 0, an integer `G`, and an integer
matrix `Qint`, meaning gamma = G/den and Q = Qint/den.

Everything is verified in EXACT INTEGER ARITHMETIC -- there are no numerical
tolerances anywhere in this task:
  (1) the polynomial identity is a finite set of integer equations,
  (2) positive definiteness is decided exactly by Sylvester's criterion,
      with the leading principal minors computed by the fraction-free
      Bareiss algorithm (integer arithmetic, no rounding).
Positive definiteness of Q is what makes the bound a *proof*: for every real
x, p(x) - gamma = z(x)^T Q z(x) >= 0.
"""

from itertools import product

N_VARS = 5
DEG_P = 6          # total degree of p
DEG_Z = 3          # z holds all monomials of degree <= DEG_Z
INT_CAP = 1 << 220  # every submitted integer must satisfy |v| < INT_CAP
DEN_CAP = 1 << 20   # the denominator may not exceed this

# Calibrated per-instance acceptance thresholds, as exact rationals
# (num, den): target = gamma*_k - 2e-5, partial = gamma*_k - 0.1, where the
# per-instance SOS optima gamma*_k were measured by bisection to 1e-9 and
# cross-checked against Newton-polished global minimisation (agreement to
# ~2e-7). With den <= 2^20 the rounding grid is coarser than the spectral
# margin at each optimum, so certificates require a deliberate
# margin-vs-tightness trade-off. See calibration.md.
GAMMA_TARGETS = {
    "p1": (-13816645, 1000000),   # gamma* = -13.8166248571
    "p2": (-86681433, 1000000),   # gamma* = -86.6814130632
    "p3": (-15831777, 1000000),   # gamma* = -15.8317566696
}
GAMMA_PARTIALS = {
    "p1": (-13917, 1000),
    "p2": (-86781, 1000),
    "p3": (-15932, 1000),
}
WITNESS_GAP_TOL = (1, 100000)          # p(x0) - gamma <= 1e-5, all instances
# backward-compat aliases (single-instance verifier paths)
GAMMA_TARGET = GAMMA_TARGETS["p1"]
GAMMA_PARTIAL = GAMMA_PARTIALS["p1"]


def monomials(n_vars, max_deg):
    """All exponent tuples with total degree <= max_deg, in a fixed order:
    sorted by (total degree, then the exponent tuple lexicographically)."""
    out = []
    for e in product(range(max_deg + 1), repeat=n_vars):
        if sum(e) <= max_deg:
            out.append(e)
    out.sort(key=lambda e: (sum(e), e))
    return out


Z_BASIS = monomials(N_VARS, DEG_Z)          # 56 monomials
P_MONOMIALS = monomials(N_VARS, DEG_P)      # 462 monomials
M = len(Z_BASIS)


def bareiss_leading_minors(A):
    """Leading principal minors of an integer symmetric matrix, exactly.

    Fraction-free (Bareiss) elimination: every intermediate value stays an
    integer, and after step k the pivot equals the (k+1)-st leading
    principal minor. Returns None if a zero pivot is hit (then the matrix is
    not positive definite, since PD forces all leading minors > 0).
    """
    n = len(A)
    a = [row[:] for row in A]
    minors = []
    prev = 1
    for k in range(n):
        piv = a[k][k]
        minors.append(piv)
        if piv <= 0:
            return minors          # not PD; caller sees a non-positive minor
        for i in range(k + 1, n):
            aik = a[i][k]
            for j in range(k + 1, n):
                a[i][j] = (a[i][j] * piv - aik * a[k][j]) // prev
            a[i][k] = 0
        prev = piv
    return minors


def is_positive_definite(Qint):
    """Exact Sylvester test: symmetric integer Q is PD iff every leading
    principal minor is strictly positive."""
    minors = bareiss_leading_minors(Qint)
    return len(minors) == len(Qint) and all(mv > 0 for mv in minors)


def identity_defect(Qint, den, G, p_coeffs):
    """Exact check of  z^T Q z == p - gamma  (as integer equations).

    Multiplying by `den`, the identity says, for every monomial alpha,
        sum over ordered pairs (i, j) with z_i + z_j = alpha of Qint[i][j]
            == den * p_alpha - G * [alpha == 0].
    Returns the number of violated monomial equations (0 means exact).
    """
    acc = {}
    rng_v = range(N_VARS)
    for i, ei in enumerate(Z_BASIS):
        row = Qint[i]
        for j, ej in enumerate(Z_BASIS):
            a = tuple(ei[v] + ej[v] for v in rng_v)
            acc[a] = acc.get(a, 0) + row[j]
    zero = (0,) * N_VARS
    bad = 0
    seen = set()
    for a, v in acc.items():
        seen.add(a)
        want = den * p_coeffs.get(a, 0) - (G if a == zero else 0)
        if v != want:
            bad += 1
    for a, c in p_coeffs.items():          # monomials z^T Q z cannot produce
        if a not in seen:
            want = den * c - (G if a == zero else 0)
            if want != 0:
                bad += 1
    if zero not in seen and zero not in p_coeffs and G != 0:
        bad += 1
    return bad


def verify_submission(sub, p_coeffs, gamma_target, gamma_partial):
    """Full exact verification (certificate + witness). `sub` is parsed JSON."""
    out = {}
    try:
        den = int(sub["den"])
        G = int(sub["G"])
        Q = [[int(v) for v in row] for row in sub["Q"]]
    except Exception as exc:
        return {"structurally_valid": False, "reason": f"malformed fields: {exc}"}

    if den <= 0:
        return {"structurally_valid": False, "reason": "den must be positive"}
    if den > DEN_CAP:
        return {"structurally_valid": False,
                "reason": f"den = {den} exceeds the cap {DEN_CAP} (= 2^20)"}
    if len(Q) != M or any(len(r) != M for r in Q):
        return {"structurally_valid": False,
                "reason": f"Q must be {M}x{M}"}
    if abs(den) >= INT_CAP or abs(G) >= INT_CAP or \
            any(abs(v) >= INT_CAP for r in Q for v in r):
        return {"structurally_valid": False,
                "reason": "an integer exceeds the size cap"}
    if any(Q[i][j] != Q[j][i] for i in range(M) for j in range(i)):
        return {"structurally_valid": False, "reason": "Q is not symmetric"}
    out["structurally_valid"] = True

    out["identity_violations"] = identity_defect(Q, den, G, p_coeffs)
    out["positive_definite"] = bool(is_positive_definite(Q))
    out["gamma"] = G / den
    out["den_bits"] = int(den).bit_length()
    out["max_entry_bits"] = max(abs(v) for r in Q for v in r).bit_length()

    gt_num, gt_den = gamma_target
    gp_num, gp_den = gamma_partial
    # exact rational comparisons  G/den >= gt_num/gt_den
    out["meets_target"] = bool(G * gt_den >= gt_num * den)
    out["meets_partial"] = bool(G * gp_den >= gp_num * den)
    out.update(verify_witness(sub, p_coeffs, G, den))
    return out


def grade(m):
    if not m.get("structurally_valid"):
        return 0.0, m.get("reason", "invalid submission")
    if m["identity_violations"] != 0:
        return 0.0, (f"{m['identity_violations']} monomial equations of "
                     "z^T Q z = p - gamma are violated")
    if not m["positive_definite"]:
        return 0.0, "Q is not positive definite (exact Sylvester test failed)"
    if not m["meets_partial"]:
        return 0.0, "valid certificate, but the bound is far too loose"
    if not m["meets_target"]:
        return 0.5, "valid exact certificate, but the bound is loose"
    if not m.get("witness_valid"):
        return 0.5, ("lower bound certified, but the witness point is missing "
                     "or malformed: " + str(m.get("reason", "")))
    if not m.get("enclosure_ok"):
        return 0.5, (f"lower bound certified, but the witness leaves an "
                     f"enclosure of width {m.get('enclosure_width'):.3e}")
    return 1.0, "exact SOS certificate plus a witness: global minimum enclosed"


# ---------------------------------------------------------------------------
# Second deliverable: a rational witness point (two-sided enclosure)
# ---------------------------------------------------------------------------
# The SOS certificate proves  min p >= gamma.  A witness point x0 with
# p(x0) small proves  min p <= p(x0).  Requiring both turns the answer into a
# certified ENCLOSURE of the global minimum, and forces the solver to do
# global optimisation as well as semidefinite programming. p has integer
# coefficients and x0 is submitted as exact rationals, so p(x0) is evaluated
# in exact rational arithmetic -- still no floating point anywhere.
WITNESS_GAP = (1, 100000)      # require p(x0) - gamma <= 1e-5, exactly


def eval_poly_exact(p_coeffs, num, den):
    """p(x0) as an exact Fraction, with x0_i = num[i]/den[i]."""
    from fractions import Fraction
    xs = [Fraction(int(n), int(d)) for n, d in zip(num, den)]
    total = Fraction(0)
    for e, c in p_coeffs.items():
        term = Fraction(int(c))
        for xi, ei in zip(xs, e):
            if ei:
                term *= xi ** int(ei)
        total += term
    return total


def verify_witness(sub, p_coeffs, G, den):
    """Check the witness point and the resulting enclosure, exactly."""
    from fractions import Fraction
    out = {}
    try:
        num = [int(v) for v in sub["x0_num"]]
        wden = [int(v) for v in sub["x0_den"]]
    except Exception as exc:
        return {"witness_valid": False, "reason": f"malformed witness: {exc}"}
    if len(num) != N_VARS or len(wden) != N_VARS:
        return {"witness_valid": False,
                "reason": f"x0 must have {N_VARS} coordinates"}
    if any(d <= 0 for d in wden):
        return {"witness_valid": False, "reason": "x0 denominators must be positive"}
    if any(abs(v) >= INT_CAP for v in num + wden):
        return {"witness_valid": False, "reason": "witness integer exceeds the cap"}
    pv = eval_poly_exact(p_coeffs, num, wden)
    gamma = Fraction(int(G), int(den))
    gap = pv - gamma
    gn, gd = WITNESS_GAP
    out["witness_valid"] = True
    out["p_at_x0"] = float(pv)
    out["enclosure_width"] = float(gap)
    out["enclosure_ok"] = bool(gap >= 0 and gap * gd <= gn)
    return out


def grade_multi_submission(subs, polys):
    """subs: {key: submission dict}; polys: {key: coeff dict}.
    Full credit only if every instance scores 1.0; 0.5 if every instance
    reaches at least 0.5; else 0. Returns (score, per_instance_report)."""
    per = {}
    scores = []
    for key in sorted(polys):
        if key not in subs:
            per[key] = {"score": 0.0, "reason": "missing instance"}
            scores.append(0.0)
            continue
        m = verify_submission(subs[key], polys[key],
                              GAMMA_TARGETS[key], GAMMA_PARTIALS[key])
        sc, reason = grade(m)
        per[key] = {"score": sc, "reason": reason, "metrics": m}
        scores.append(sc)
    if all(s == 1.0 for s in scores):
        return 1.0, per
    if all(s >= 0.5 for s in scores):
        return 0.5, per
    return 0.0, per