# Restriction map r(m_mu)(y) = sum_{nu in W mu} q^{<rho_L,nu>} y^{phi(nu)} for E6/F4; windows; separating functional.
import sys, pickle, os; sys.path.insert(0, '/private/tmp/claude-501/-Users-miltonlin-Projects-Trees/75c2c3d9-ef06-4644-b59b-d12583740b0c/scratchpad/anim/work/G5/scripts')
from fractions import Fraction as F
from collections import defaultdict, Counter
from e6model import EXP
A = EXP; n = 6
# inverse Cartan matrix (exact)
def inv(M):
    m = len(M); X = [[F(M[i][j]) for j in range(m)] + [F(int(i == j)) for j in range(m)] for i in range(m)]
    for c in range(m):
        p = next(r for r in range(c, m) if X[r][c] != 0); X[c], X[p] = X[p], X[c]
        pv = X[c][c]; X[c] = [x / pv for x in X[c]]
        for r in range(m):
            if r != c and X[r][c] != 0:
                f = X[r][c]; X[r] = [x - f * y for x, y in zip(X[r], X[c])]
    return [[X[i][m + j] for j in range(m)] for i in range(m)]
Ainv = inv(A)
rhoL_on = [-3, 1, 1, 1, 1, -3]          # <rho_L, alpha_j^vee>, verified in e6_theta.py
def rhoL(c): return sum(rhoL_on[j] * c[j] for j in range(n))
def ab_of_coroot_vec(c):  # phi(c) = c_1 e1 + c_6 e6 (phi kills black coroots); e1=(2,-1), e6=(-1,2)
    return (2 * c[0] - c[5], -c[0] + 2 * c[5])
def s(j, c):
    p = sum(A[j][i] * c[i] for i in range(n)); c = list(c); c[j] -= p; return tuple(c)
def orbit(c):
    seen = {c}; fr = [c]
    while fr:
        nf = []
        for x in fr:
            for j in range(n):
                y = s(j, x)
                if y not in seen: seen.add(y); nf.append(y)
        fr = nf
    return seen
def norm(c): return sum(A[i][j] * c[i] * c[j] for i in range(n) for j in range(n))
def dominant_coweights(maxnorm):
    # lam = sum_i f_i omega_i^vee with f_i >= 0 integers, in coroot coords lam = Ainv f ; keep integral ones
    out = []
    bound = 20
    def rec(i, f):
        if i == n:
            lam = [sum(Ainv[k][j] * f[j] for j in range(n)) for k in range(n)]
            if all(x.denominator == 1 for x in lam):
                lam = tuple(int(x) for x in lam); nm = norm(lam)
                if nm <= maxnorm: out.append((tuple(f), lam, nm))
            return
        for v in range(bound):
            f2 = f + [v]
            # norm of partial >= (4/3)*v^2 lower bound per coordinate is enough to prune crudely
            if sum(Ainv[k][k] * f2[k] ** 2 for k in range(len(f2))) > maxnorm: break
            rec(i + 1, f2)
    rec(0, [])
    return out
def s1(k): a, b = k; return (-a, a + b)
def s6(k): a, b = k; return (a + b, -b)
def wx_orbit(k):
    seen = {k}; fr = [k]
    while fr:
        nf = []
        for x in fr:
            for y in (s1(x), s6(x)):
                if y not in seen: seen.add(y); nf.append(y)
        fr = nf
    return seen
def is_dom(k): return k[0] >= 0 and k[1] >= 0
def r_of(lam):
    """dict: W_X-dominant kappa -> Counter{exponent of q: multiplicity}"""
    res = defaultdict(Counter)
    for nu in orbit(lam):
        k = ab_of_coroot_vec(nu)
        if is_dom(k): res[k][rhoL(nu)] += 1
    return dict(res)
def poly_str(cn):
    return ' + '.join(('%d*' % m if m != 1 else '') + 'q^%d' % e for e, m in sorted(cn.items()))
def ev(cn, q): return sum(m * F(q) ** e for e, m in cn.items())