"""Deeper analysis for the main session.
 (C) classification by equality: at L = f(n), is the system feasible WITHOUT a hyperplane class?
 (B) phantom split distribution at n=48, L=253 (integer DP with backpointers), and the 4th-moment test.
 (A) uniform convexity bound: L <= K(n) + n(81-n)/162 * sqrt(Dmax/6), Dmax = max admissible sum (k_i-n/3)^2.
"""
import sys
from math import comb, sqrt
from fractions import Fraction as Fr
import numpy as np
import bound
from bound import g, partitions, U, e2, abc, lb_construction

NEG = -10**9

def dp_max(n, Lstar, allowed):
    """max sum abc over 40-multisets of allowed splits with sum e2 = 27C(n,2); returns (value, multiset)."""
    P = allowed
    E = 27 * comb(n, 2)
    dp = np.full((41, E + 1), NEG, dtype=np.int64); dp[0, 0] = 0
    choice = np.full((41, E + 1), -1, dtype=np.int16)
    for step in range(40):
        for idx, p in enumerate(P):
            e, v = e2(p), abc(p)
            cand = dp[step, :E + 1 - e] + v
            better = cand > dp[step + 1, e:]
            dp[step + 1, e:][better] = cand[better]
            choice[step + 1, e:][better] = idx
    if dp[40, E] <= NEG // 2: return None, None
    # backtrack
    ms = []; s = E
    for step in range(40, 0, -1):
        idx = choice[step, s]; ms.append(P[idx]); s -= e2(P[idx])
    return int(dp[40, E]), ms

def feasible(n, L, allowed):
    v, ms = dp_max(n, L, allowed)
    return v is not None and v >= 18 * L + 9 * comb(n, 3), ms

def partC():
    print("== (C) classification by equality ==")
    for n in range(3, 48):
        L = lb_construction(n)
        adm = [p for p in partitions(n) if U(p) >= L]
        if n <= 27:
            forbid = lambda p: p[0] == n           # the hyperplane split (n,0,0)
            label = "(n,0,0)"
        else:
            forbid = lambda p: p[0] == 27          # a full hyperplane inside S
            label = "a 27-part"
        ok_all, _ = feasible(n, L, adm)
        ok_without, ms = feasible(n, L, [p for p in adm if not forbid(p)])
        print(f"n={n:2d} L={L:4d}: feasible with all admissible splits: {ok_all}; "
              f"feasible with {label} forbidden: {ok_without}" + (f"  e.g. {sorted(set(ms))}" if ok_without else ""))

def moments_test(n, L, ms):
    """4th-moment identity: sum_W sum_i C(k_i,4) = C(n,4) + 3Q, Q = L(n-3) + sum_z C(m_z,2) >= L(n-3) + J(L)."""
    F = sum(comb(k, 4) for p in ms for k in p)
    # integer Jensen: sum over z in S of C(m_z,2) with sum m_z = 3L over n points; z not in S: sum = C(n,2)-3L over 81-n
    def jmin(total, cnt):
        q, r = divmod(total, cnt); return r * comb(q + 1, 2) + (cnt - r) * comb(q, 2)
    J = jmin(3 * L, n) + jmin(comb(n, 2) - 3 * L, 81 - n)
    need = comb(n, 4) + 3 * (L * (n - 3) + J)
    return F, need

def partB():
    print("\n== (B) phantom at n=48, L=253 ==")
    n, L = 48, 253
    adm = [p for p in partitions(n) if U(p) >= L]
    v, ms = dp_max(n, L, adm)
    print("max sum abc =", v, " needed >=", 18 * L + 9 * comb(n, 3), " implied L =", Fr(v - 9 * comb(n, 3), 18))
    from collections import Counter
    print("phantom multiset:", sorted(Counter(ms).items(), reverse=True))
    F, need = moments_test(n, L, ms)
    print(f"4th moment of phantom: sum C(k,4) = {F}; lower bound from B4/Jensen: {need}; violated: {F < need}")
    # and the real construction H + 11 + 10 for comparison (its split distribution via a concrete set)
    for nn in (48, 49, 50, 51, 52):
        LL = lb_construction(nn) + 1
        admn = [p for p in partitions(nn) if U(p) >= LL]
        vv, mss = dp_max(nn, LL, admn)
        if vv is None: print(f"n={nn}: infeasible at L={LL}"); continue
        FF, nd = moments_test(nn, LL, mss)
        print(f"n={nn} L={LL}: DP feasible ({Fr(vv - 9*comb(nn,3),18)} >= {LL}); 4th-moment test on argmax multiset: {FF} vs need {nd} -> {'violated' if FF < nd else 'satisfied'}")

def partA():
    print("\n== (A) uniform convexity bound ==")
    print("identity: L = K(n) + (1/18) sum_W (a-n/3)(b-n/3)(c-n/3),  K(n) = n(n^2-1)/6 - 40 n^3/243")
    print("pair identity: sum_W sum_i (k_i - n/3)^2 = n(81-n)/3 =: B(n)")
    for n in range(10, 48):
        t = lb_construction(n) + 1
        adm = [p for p in partitions(n) if U(p) >= t]
        K = Fr(n * (n * n - 1), 6) - Fr(40 * n ** 3, 243)
        B = Fr(n * (81 - n), 3)
        if not adm:
            print(f"n={n:2d} t={t}: no admissible split"); continue
        Dmax = max(sum((Fr(k) - Fr(n, 3)) ** 2 for k in p) for p in adm)
        Dmin = min(sum((Fr(k) - Fr(n, 3)) ** 2 for k in p) for p in adm)
        if 40 * Dmax < B:
            print(f"n={n:2d} t={t}: 40*Dmax = {float(40*Dmax):.1f} < B = {float(B):.1f}  (pair-count contradiction)")
            continue
        if 40 * Dmin > B:
            print(f"n={n:2d} t={t}: 40*Dmin = {float(40*Dmin):.1f} > B = {float(B):.1f}  (pair-count contradiction)")
            continue
        ub = float(K) + float(B) * 2 * (float(Dmax) / 6) ** 1.5 / (18 * float(Dmax))
        print(f"n={n:2d} t={t}: K={float(K):.2f} Dmax={float(Dmax):.2f} -> convexity bound L <= {ub:.2f} {'< t OK' if ub < t else '>= t FAILS'}")

if __name__ == "__main__":
    which = sys.argv[1] if len(sys.argv) > 1 else "ABC"
    if "C" in which: partC()
    if "B" in which: partB()
    if "A" in which: partA()