# G7 independent check of the E6/F4 root-datum computations used by R-3, A7-19, A7-25, A7-32 and the verdict page.
# Standard library only. Bourbaki labels: chain 1-3-4-5-6, node 2 attached to 4.
from fractions import Fraction as Fr
from itertools import product

N = 6
edges = [(1,3),(3,4),(4,5),(5,6),(2,4)]
C = [[0]*N for _ in range(N)]
for i in range(N): C[i][i] = 2
for a,b in edges:
    C[a-1][b-1] = -1; C[b-1][a-1] = -1

# Simply laced: identify roots/coroots; weights/coweights in fundamental coordinates.
# Vectors in "simple-root coordinates" (coefficients on alpha_j).  Pairing <x,y> = x^T C y in root coords.
def pair_rr(x, y):
    return sum(x[i]*C[i][j]*y[j] for i in range(N) for j in range(N))

def root_to_fund(x):  # coordinates <x, alpha_j> for j
    return [sum(x[i]*C[i][j] for i in range(N)) for j in range(N)]

# inverse Cartan matrix (rational)
def inv(M):
    n = len(M); A = [[Fr(M[i][j]) for j in range(n)] + [Fr(int(i==j)) for j in range(n)] for i in range(n)]
    for c in range(n):
        p = next(r for r in range(c,n) if A[r][c] != 0)
        A[c],A[p] = A[p],A[c]
        pv = A[c][c]; A[c] = [v/pv for v in A[c]]
        for r in range(n):
            if r != c and A[r][c] != 0:
                f = A[r][c]; A[r] = [A[r][k]-f*A[c][k] for k in range(2*n)]
    return [row[n:] for row in A]
Cinv = inv(C)
def fund_to_root(c):  # weight with fundamental coords c -> root coords
    return [sum(Cinv[i][j]*c[j] for j in range(N)) for i in range(N)]

omega = [fund_to_root([int(i==j) for i in range(N)]) for j in range(N)]  # omega_j in root coords
print("omega_1 in root coords:", [str(v) for v in omega[0]])
print("omega_6 in root coords:", [str(v) for v in omega[5]])

# positive roots
def all_roots():
    simple = [tuple(int(i==j) for i in range(N)) for j in range(N)]
    roots = set(simple); frontier = list(simple)
    while frontier:
        new = []
        for r in frontier:
            for j in range(N):
                f = root_to_fund(list(r))
                # r + alpha_j is a root if <r,alpha_j> < 0 ... use string: for simply laced, r+a_j root iff <r,a_j> = -1
                if f[j] == -1:
                    s = list(r); s[j] += 1; s = tuple(s)
                    if s not in roots:
                        roots.add(s); new.append(s)
        frontier = new
    return sorted(roots)
pos = all_roots()
print("number of positive roots:", len(pos))
assert len(pos) == 36

# Weyl group action on fundamental-coordinate vectors: s_j(c) = c - c_j * (row j of C)
def refl(c, j):
    cj = c[j]
    return tuple(c[k] - cj*C[j][k] for k in range(N))

def dominant_rep(c):
    c = tuple(c); word = []
    while True:
        neg = [j for j in range(N) if c[j] < 0]
        if not neg: return c, word
        j = neg[0]; c = refl(c, j); word.append(j+1)

lam1 = (-8,1,1,1,1,-2); lam2 = (-8,1,1,1,1,1)
d1,w1 = dominant_rep(lam1); d2,w2 = dominant_rep(lam2)
print("dominant rep of rho-9w1-3w6:", d1, "word length", len(w1))
print("dominant rep of rho-9w1    :", d2, "word length", len(w2))
assert d1 == d2 == (1,1,1,0,1,1)
# norm check (in root coords)
r1 = fund_to_root(list(lam1)); r2 = fund_to_root(list(lam2))
print("norms:", pair_rr(r1,r1), pair_rr(r2,r2))

# theta = -w0(D4) on X_*(T); black nodes 2,3,4,5.  Compute via projection: theta(x) = 2 P_black x - x.
black = [1,2,3,4]  # zero-based indices of nodes 2,3,4,5
Bvec = [[int(i==j) for i in range(N)] for j in black]  # black simple roots in root coords
G = [[pair_rr(Bvec[a],Bvec[b]) for b in range(4)] for a in range(4)]
Ginv = inv(G)
def P_black(x):
    coeffs = [pair_rr(Bvec[a], x) for a in range(4)]
    sol = [sum(Ginv[a][b]*coeffs[b] for b in range(4)) for a in range(4)]
    out = [Fr(0)]*N
    for a in range(4):
        for i in range(N): out[i] += sol[a]*Bvec[a][i]
    return out
def theta(x):
    pb = P_black(x); return [2*pb[i]-x[i] for i in range(N)]
def phi(x):
    t = theta(x); return [x[i]-t[i] for i in range(N)]
# theta integral and an involution, permutes roots
for r in pos:
    t = theta(list(r))
    assert all(v.denominator == 1 for v in t)
    tt = tuple(int(v) for v in t)
    assert tt in set(pos) or tuple(-v for v in tt) in set(pos)
print("theta preserves the root system: ok")
fixed_black_roots = [r for r in pos if tuple(int(v) for v in theta(list(r))) == r]
print("positive roots fixed by theta:", len(fixed_black_roots))
e1 = phi([1,0,0,0,0,0]); e6 = phi([0,0,0,0,0,1])
print("e1 = phi(alpha1^v) in coroot coords:", [str(v) for v in e1])
print("e6 = phi(alpha6^v) in coroot coords:", [str(v) for v in e6])

# rho_L: half sum of positive black roots; <rho_L, alpha_j^v>
black_pos = [r for r in pos if all(r[i] == 0 for i in (0,5))]
print("positive black roots:", len(black_pos))
twoRhoL = [sum(r[i] for r in black_pos) for i in range(N)]
print("2 rho_L root coords:", twoRhoL, " <rho_L,alpha_j^v>:", [Fr(v,2) for v in root_to_fund(twoRhoL)])

# Lambda_T = coroot lattice Q (E6 simply connected).  Lambda_S = Q cap black-perp.
# Find integer vectors in Q (root coords integers) orthogonal to black: solve <x,alpha_b>=0
# brute force small box to get lattice basis
cand = []
for x in product(range(-6,7), repeat=N):
    if all(root_to_fund(list(x))[b] == 0 for b in black) and any(x):
        cand.append(x)
# check whether every such vector is an integer combination of e1, e6
def in_span_e(x):
    # solve x = c1 e1 + c6 e6 using coordinates 1 and 6 (e1,e6 have coordinate pattern)
    # e1, e6 as vectors
    a11, a12 = e1[0], e6[0]; a21, a22 = e1[5], e6[5]
    det = a11*a22 - a12*a21
    c1 = (x[0]*a22 - a12*x[5]) / det; c6 = (a11*x[5] - a21*x[0]) / det
    ok = all(c1*e1[i] + c6*e6[i] == x[i] for i in range(N))
    return ok, c1, c6
bad = [x for x in cand if not (in_span_e(x)[0] and in_span_e(x)[1].denominator == 1 and in_span_e(x)[2].denominator == 1)]
print("vectors of Q cap black-perp in box:", len(cand), " not in Z e1 + Z e6:", len(bad))

# restricted coordinates (a,b) with kappa = c1 e1 + c6 e6, (a,b) = (2c1 - c6, -c1 + 2c6)
def ab_of(x):
    ok, c1, c6 = in_span_e(x)
    assert ok
    return (2*c1 - c6, -c1 + 2*c6)

# phi of fundamental coweights
for j in range(N):
    print("phi(omega_%d^v) in (a,b):" % (j+1), tuple(str(v) for v in ab_of(phi(omega[j]))))

# dominance of kappa in E6: <alpha_j, kappa> >= 0
def is_dom_root(x):
    return all(v >= 0 for v in root_to_fund(x))

# Lambda_S^+ in (a,b): kappa = c1 e1 + c6 e6 with integer c, dominant
LS_plus = set()
for c1 in range(-20, 41):
    for c6 in range(-20, 41):
        x = [c1*e1[i] + c6*e6[i] for i in range(N)]
        if is_dom_root(x):
            LS_plus.add(ab_of(x))
LS_box = sorted(k for k in LS_plus if k[0] <= 18 and k[1] <= 18)
# phi(Lambda_T^+): lambda = sum c_j omega_j with lambda in Q (root coords integral)
img = set()
B = 10
for c in product(range(B), repeat=N):
    lamr = fund_to_root(list(c))
    if all(v.denominator == 1 for v in lamr):
        k = ab_of(phi(lamr))
        img.add(k)
unc = [k for k in LS_box if k not in img]
print("Lambda_S^+ labels with a,b<=18:", len(LS_box))
print("uncovered labels (a,b<=18):", unc)

# Hecke labels lambda_k = (3k-1)/2 omega1 + omega6, k odd
for k in (1,3,5):
    c = [Fr(3*k-1,2),0,0,0,0,1]
    lamr = fund_to_root(c)
    assert all(v.denominator == 1 for v in lamr), "lambda_k not in coroot lattice"
    ab = ab_of(phi(lamr))
    two_rho = 2*sum(lamr)  # <2rho, lambda> = 2 * sum of root coords (since <rho,alpha_i^v>=1)
    mu = (3*k, 0); mup = (3*k-1, 2)
    # delta(a,b) = <rho, kappa>
    x_mup = [ (Fr(2*mup[0]+mup[1],3))*e1[i] + (Fr(mup[0]+2*mup[1],3))*e6[i] for i in range(N)]
    x_mu  = [ (Fr(2*mu[0]+mu[1],3))*e1[i] + (Fr(mu[0]+2*mu[1],3))*e6[i] for i in range(N)]
    print("k=%d: phi(lambda_k)=%s, mu'_k=%s, <2rho,lambda_k>=%s, delta(mu'_k)=%s, codim=%s" % (
        k, tuple(str(v) for v in ab), mup, two_rho, sum(x_mup), sum(x_mup)-sum(x_mu)))
    assert ab == mup and two_rho == sum(x_mup) and sum(x_mup)-sum(x_mu) == 8

# W-orbit of lambda_1 and fixed points phi(nu) = mu_1, mu'_1
def orbit(c):
    c = tuple(c); seen = {c}; fr = [c]
    while fr:
        new = []
        for v in fr:
            for j in range(N):
                w = refl(v, j)
                if w not in seen: seen.add(w); new.append(w)
        fr = new
    return seen
for k in (1,3):
    lam_f = (Fr(3*k-1,2),0,0,0,0,1)
    orb = orbit(lam_f)
    cnt_mu = 0; cnt_mup = 0; rhoL_vals = []
    for v in orb:
        vr = fund_to_root(list(v))
        ab = ab_of(phi(vr))
        if ab == (3*k,0):
            cnt_mu += 1; rhoL_vals.append(Fr(pair_rr(twoRhoL, vr),2))
        if ab == (3*k-1,2): cnt_mup += 1
    print("k=%d: |W lambda_k|=%d, #phi(nu)=mu_k: %d, #phi(nu)=mu'_k: %d, <rho_L,nu> on the first set: %s" % (
        k, len(orb), cnt_mu, cnt_mup, sorted(str(x) for x in rhoL_vals)))

# W_X = S_3 on log coordinates (L1, L6) of y with u=y(e1), v=y(e6): s1:(L1,L6)->(-L1,L1+L6), s6:(L1,L6)->(L1+L6,-L6)
def wx_orbit(p):
    seen = {p}; fr = [p]
    while fr:
        new = []
        for (a,b) in fr:
            for q in ((-a, a+b), (a+b, -b)):
                if q not in seen: seen.add(q); new.append(q)
        fr = new
    return seen
O1 = wx_orbit((-5,1)); O2 = wx_orbit((-5,4))
print("W_X y1:", sorted(O1)); print("W_X y2:", sorted(O2)); print("disjoint:", O1.isdisjoint(O2))
print("sigma y1 = (1,-5) in W_X y2:", (1,-5) in O2)
# check that x = q^{rho_L} y gives the weights lam1, lam2: <lam, alpha_j^v> = <rho_L,alpha_j^v> + L(phi(alpha_j^v))
rhoL_f = [Fr(v,2) for v in root_to_fund(twoRhoL)]
for (L1,L6), lam in (((-5,1), lam1), ((-5,4), lam2)):
    vals = [rhoL_f[j] for j in range(N)]
    vals[0] += L1; vals[5] += L6
    print("weight of q^{rho_L} y for y=(%d,%d):" % (L1,L6), tuple(str(v) for v in vals), "expected", lam)
    assert tuple(vals) == tuple(Fr(v) for v in lam)
# orbit sum M_(3,0) exponents at y1,y2: y^{(a,b)} = u^{c1} v^{c6}
def M_exponents(label, L):
    exps = []
    for (a,b) in sorted(wx_orbit(label)):
        c1 = Fr(2*a+b,3); c6 = Fr(a+2*b,3)
        exps.append(c1*L[0] + c6*L[1])
    return sorted(exps)
# W_X acting on labels (a,b) of Lambda_S: s1(a,b)=(-a,a+b), s6(a,b)=(a+b,-b) (same formulas in (a,b) coordinates)
print("M_(3,0) exponents at y1:", [str(e) for e in M_exponents((3,0),(-5,1))])
print("M_(3,0) exponents at y2:", [str(e) for e in M_exponents((3,0),(-5,4))])
print("M_(1,1) at y1,y2:", [str(e) for e in M_exponents((1,1),(-5,1))], [str(e) for e in M_exponents((1,1),(-5,4))])