"""G6 independent re-derivation of the E6/F4 numbers used in A8-2..A8-9, A7-26..A7-31.
Standard library only.  Conventions: Bourbaki labels 1..6 (chain 1-3-4-5-6, node 2 on node 4).
Coroots are written in the basis alpha_1^vee..alpha_6^vee (Lambda_T = coroot lattice, G simply connected).
Simply laced: identify roots and coroots, (x,y) = x^T A y in simple-(co)root coordinates.
"""
from fractions import Fraction as Fr
from itertools import product
from collections import Counter

A = [[2,0,-1,0,0,0],[0,2,0,-1,0,0],[-1,0,2,-1,0,0],[0,-1,-1,2,-1,0],[0,0,0,-1,2,-1],[0,0,0,0,-1,2]]
n = 6
def ip(x, y):  # (x,y) in simple-root coords
    return sum(x[i]*A[i][j]*y[j] for i in range(n) for j in range(n))
def pair(alpha, lam):  # <alpha, lam> = (alpha, lam) simply laced
    return ip(alpha, lam)
def add(x, y): return tuple(a+b for a, b in zip(x, y))
def sub(x, y): return tuple(a-b for a, b in zip(x, y))
def scal(c, x): return tuple(c*a for a in x)
def refl(i, x):  # s_i(x) = x - <alpha_i, x> alpha_i^vee
    c = sum(A[i][j]*x[j] for j in range(n))
    y = list(x); y[i] -= c; return tuple(y)

# --- roots by closure under reflections
simple = [tuple(1 if j == i else 0 for j in range(n)) for i in range(n)]
roots = set(simple); frontier = list(simple)
while frontier:
    new = []
    for r in frontier:
        for i in range(n):
            s = refl(i, r)
            if s not in roots:
                roots.add(s); new.append(s)
    frontier = new
roots = sorted(roots)
pos = [r for r in roots if all(c >= 0 for c in r)]
print('roots', len(roots), 'positive', len(pos))
two_rho = tuple(sum(r[i] for r in pos) for i in range(n))
print('2rho (simple-root coords):', two_rho)

# --- fundamental coweights: inverse Cartan matrix (rational)
def inv(M):
    m = len(M); Mx = [[Fr(M[i][j]) for j in range(m)] + [Fr(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 Mx[r][c] != 0); Mx[c], Mx[p] = Mx[p], Mx[c]
        pv = Mx[c][c]; Mx[c] = [v/pv for v in Mx[c]]
        for r in range(m):
            if r != c and Mx[r][c] != 0:
                f = Mx[r][c]; Mx[r] = [a-f*b for a, b in zip(Mx[r], Mx[c])]
    return [row[m:] for row in Mx]
Ainv = inv(A)
omega = [tuple(Ainv[i][j] for j in range(n)) for i in range(n)]  # omega_i^vee in coroot coords (symmetric)
for i in range(n):
    assert all(pair(simple[j], omega[i]) == (1 if i == j else 0) for j in range(n))
print('omega1^vee =', [str(c) for c in omega[0]], ' omega6^vee =', [str(c) for c in omega[5]])
print('<rho,omega1^vee>, <rho,omega6^vee> =', Fr(pair(two_rho, omega[0]), 2), Fr(pair(two_rho, omega[5]), 2))

# --- theta = -w_L on the coroot lattice, w_L = longest element of W(D4 on nodes 2..5)
black = [1, 2, 3, 4]  # 0-based indices of nodes 2,3,4,5
def wL(x):
    # w_L = -1 on span(black coroots), identity on the orthogonal complement span(omega1, omega6)
    # decompose x = p + q with p in black span: solve <alpha_j, p> = <alpha_j, x> for j in black
    M = [[A[i][j] for j in black] for i in black]
    Mi = inv(M)
    rhs = [sum(A[i][k]*x[k] for k in range(n)) for i in black]
    coeff = [sum(Mi[a][b]*rhs[b] for b in range(4)) for a in range(4)]
    p = [Fr(0)]*n
    for a, j in enumerate(black): p[j] = coeff[a]
    p = tuple(p); q = sub(x, p)
    return sub(q, p)
def theta(x): return scal(-1, wL(x))
Th = [theta(simple[j]) for j in range(n)]  # columns
print('theta columns (theta alpha_j^vee):')
for j in range(n): print('  ', j+1, [int(c) for c in Th[j]])
for r in roots:
    assert theta(theta(r)) == r and ip(theta(r), theta(r)) == ip(r, r)
assert all(theta(simple[j]) == simple[j] for j in black)
assert theta(omega[0]) == scal(-1, omega[0]) and theta(omega[5]) == scal(-1, omega[5])
def phi(x): return sub(x, theta(x))
def ab(x):  # coordinates in omega1^vee, omega6^vee (for x in the -1 eigenspace)
    a = pair(simple[0], x); b = pair(simple[5], x)
    assert all(pair(simple[j], x) == 0 for j in black), x
    return (a, b)
e1 = phi(simple[0]); e6 = phi(simple[5])
print('e1 = phi(alpha1^vee) =', [int(c) for c in e1], '(a,b)=', ab(e1), '; e6 =', [int(c) for c in e6], ab(e6))
for j in black: assert phi(simple[j]) == (0,)*n
# theta permutes Phi+ \ Phi_L into Phi^-
PhiL = [r for r in roots if r[0] == 0 and r[5] == 0]
print('|Phi_L| =', len(PhiL))
for r in pos:
    if r not in PhiL: assert all(c <= 0 for c in theta(r))
    else: assert theta(r) == r

def lam_ab(a, b):  # a*omega1 + b*omega6 in coroot coords
    return add(scal(a, omega[0]), scal(b, omega[5]))
def in_lattice(x): return all(Fr(c).denominator == 1 for c in x)
def in_Qplus(x): return in_lattice(x) and all(c >= 0 for c in x)
def dominant(x): return all(pair(simple[i], x) >= 0 for i in range(n))

# --- A8-2 (i): Lambda_S cap Q^+ = Z>=0 e1 + Z>=0 e6, checked on a box
ok = True
for a in range(-12, 13):
    for b in range(-12, 13):
        x = lam_ab(a, b)
        if not in_lattice(x): continue
        inQ = in_Qplus(x)
        # coordinates in e1,e6: (a,b) = c1*(2,-1) + c6*(-1,2) -> c1=(2a+b)/3, c6=(a+2b)/3
        c1 = Fr(2*a+b, 3); c6 = Fr(a+2*b, 3)
        incone = c1.denominator == 1 and c6.denominator == 1 and c1 >= 0 and c6 >= 0
        if inQ != incone: ok = False; print('mismatch', a, b)
print('A8-2(i) Lambda_S cap Q+ = N e1 + N e6 on box |a|,|b|<=12:', ok)
# lattice Lambda_S = {a = b mod 3}
print('Lambda_S membership a=b mod 3 on box:', all(in_lattice(lam_ab(a, b)) == ((a-b) % 3 == 0) for a in range(-9, 10) for b in range(-9, 10)))

def delta(x): return Fr(pair(two_rho, x), 2)
print('delta(3,0),delta(2,2),delta(1,1),delta(0,3) =', [delta(lam_ab(*p)) for p in [(3,0),(2,2),(1,1),(0,3)]])
assert all(delta(lam_ab(a, b)) == 8*(a+b) for a in range(0, 10) for b in range(0, 10))

def leq(x, y): return in_Qplus(sub(y, x))   # x <= y in dominance order of Lambda_T

# --- Weyl orbit
def worbit(x):
    seen = {x}; fr = [x]
    while fr:
        new = []
        for y in fr:
            for i in range(n):
                z = refl(i, y)
                if z not in seen: seen.add(z); new.append(z)
        fr = new
    return seen

rhoL2 = tuple(sum(r[i] for r in PhiL if all(c >= 0 for c in r)) for i in range(n))  # 2 rho_L
print('2rho_L =', rhoL2)

# dominant elements below lambda (all dominant lambda'' with lambda - lambda'' in Q+)
def dominant_below(lam):
    # enumerate kappa in Q+ with lam - kappa dominant; bound coefficients via (lam,lam) >= (lam'',lam'')? use explicit box from lam coords
    res = []
    bound = [int(c) for c in lam]
    for kappa in product(*[range(b+1) for b in bound]):
        x = sub(lam, kappa)
        if dominant(x): res.append(x)
    return res

for k in [1, 3, 5]:
    ck = Fr(3*k-1, 2)
    lam = lam_ab(ck, 1)
    mu = lam_ab(3*k, 0); mup = lam_ab(3*k-1, 2)
    assert in_lattice(lam) and dominant(lam)
    assert phi(lam) == mup and sub(mup, mu) == e6
    D = pair(two_rho, lam)
    # A8-2(iii): interval [mu, mu'] in Lambda_S^+
    interval = [(a, b) for a in range(0, 3*k+3) for b in range(0, 6)
                if in_lattice(lam_ab(a, b)) and leq(mu, lam_ab(a, b)) and leq(lam_ab(a, b), mup)]
    # A8-4: dominant lambda'' <= lambda with phi(lambda'') >= mu (any such g-cell can meet the fibre over x_mu)
    below = dominant_below(lam)
    hits = [x for x in below if leq(mu, phi(x))]
    hits_p = [x for x in below if leq(mup, phi(x))]
    orb = worbit(tuple(int(c) for c in lam))
    Fk = [nu for nu in orb if phi(nu) == mu]
    Fkp = [nu for nu in orb if phi(nu) == mup]
    grading = sorted(Fr(pair(rhoL2, nu), 2) for nu in Fk)
    print(f'k={k}: lambda_k={tuple(int(c) for c in lam)}, <2rho,lambda>={D}, delta(mu\')={delta(mup)}, delta(mu)={delta(mu)},'
          f' <rho_L,lambda>={Fr(pair(rhoL2, lam),2)}; interval={interval}')
    print(f'   #dominant lambda\'\'<=lambda_k: {len(below)}; those with phi>=mu_k: {[tuple(int(c) for c in x) for x in hits]};'
          f' with phi>=mu_k\': {len(hits_p)}')
    print(f'   |W lambda_k|={len(orb)}, #nu with phi(nu)=mu_k: {len(Fk)}, #nu with phi(nu)=mu_k\': {len(Fkp)} (={Fkp==[tuple(int(c) for c in lam)]}); <rho_L,nu>: {[str(g) for g in grading]}')
    # A8-6(d) Lie algebra counts at each fixed point
    D5 = [r for r in roots if r[0] == 0]   # Levi of mu_k: roots with no alpha_1 component
    assert all(pair(r, mu) == 0 for r in D5) and len(D5) == 40
    outside = [r for r in D5 if r not in PhiL]
    pairs = set()
    for r in outside:
        t = tuple(int(c) for c in theta(r)); assert t in roots and t != r and t in outside
        pairs.add(frozenset([r, t]))
    dimZ = 4 + len(PhiL) + len(pairs)
    counts = set()
    for nu in Fk:
        longneg = sum(1 for r in PhiL if pair(r, nu) <= 0)
        pairneg = sum(1 for P in pairs if all(pair(r, nu) <= 0 for r in P))
        counts.add((len(PhiL)-longneg, len(pairs)-pairneg, dimZ - (4+longneg+pairneg)))
    print(f'   dim Z_k = 4+{len(PhiL)}+{len(pairs)} = {dimZ}; (#long>0, #pairs with a positive member, orbit-dim lower bound) at fixed points: {counts}')

# k=1 Schubert stratification of Gr^{<= lambda_1}
lam1 = lam_ab(1, 1)
print('dominant coweights <= lambda_1 and dims:', [(tuple(int(c) for c in x), pair(two_rho, x)) for x in dominant_below(lam1)])

# mirror: mu=(0,3k), lambda = omega1 + c_k omega6, mu' = mu + e1
for k in [1, 3, 5]:
    ck = Fr(3*k-1, 2); lam = lam_ab(1, ck); mu = lam_ab(0, 3*k); mup = lam_ab(2, 3*k-1)
    assert phi(lam) == mup and sub(mup, mu) == e1 and in_lattice(lam) and dominant(lam)
    below = dominant_below(lam)
    hits = [x for x in below if leq(mu, phi(x))]
    orb = worbit(tuple(int(c) for c in lam))
    Fk = [nu for nu in orb if phi(nu) == mu]; Fkp = [nu for nu in orb if phi(nu) == mup]
    D5m = [r for r in roots if r[5] == 0]
    outside = [r for r in D5m if r not in PhiL]
    pairs = {frozenset([r, tuple(int(c) for c in theta(r))]) for r in outside}
    counts = set()
    for nu in Fk:
        longneg = sum(1 for r in PhiL if pair(r, nu) <= 0)
        pairneg = sum(1 for P in pairs if all(pair(r, nu) <= 0 for r in P))
        counts.add(4+len(PhiL)+len(pairs) - (4+longneg+pairneg))
    print(f'mirror k={k}: hits={len(hits)} (only lambda itself: {hits==[lam]}), #F={len(Fk)}, #F\'={len(Fkp)}, orbit-dim bounds {counts}')

# A7-28: no root vanishes on T' = torus of the black coroots
vanish = [r for r in roots if all(pair(r, simple[j]) == 0 for j in black)]
print('roots vanishing on T\':', len(vanish))

# A8-9(c) / family page: uncovered labels in a box
cover = set()
B = 16
for c in product(range(0, 9), repeat=6):
    x = tuple(sum(c[i]*omega[i][j] for i in range(n)) for j in range(n))
    if not in_lattice(x): continue
    a, b = ab(phi(x))
    if a <= B and b <= B: cover.add((a, b))
labels = [(a, b) for a in range(B+1) for b in range(B+1) if (a-b) % 3 == 0]
unc = [p for p in labels if p not in cover]
print('uncovered labels with a,b<=16:', unc)

# A8-11.3: pairing of lambda_1 with the 27 weights
w27 = worbit(tuple(omega[0]))
print('|W omega_1| =', len(w27), ' <w,lambda_1> on 27:', sorted(Counter(pair(w, lam1) for w in w27).items()))