# E6/F4: Borel orbits via conjugates of w_L, types of (Y, alpha), stabilizer dimensions.
# Standard library only. Roots in simple-root coordinates, Bourbaki labels (index 0..5 = nodes 1..6).
from fractions import Fraction
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 pair(b, i):  # <beta, alpha_i^vee>
    return sum(b[j]*A[j][i] for j in range(n))
def refl(b, i):
    c = pair(b, i)
    return tuple(b[j] - (c if j == i else 0) for j in range(n))
simple = [tuple(1 if j == i else 0 for j in range(n)) for i in range(n)]
pos = set(simple); fr = list(simple)
while fr:
    new = []
    for b in fr:
        for i in range(n):
            r = refl(b, i)
            if all(x >= 0 for x in r) and r not in pos:
                pos.add(r); new.append(r)
    fr = new
pos = sorted(pos, key=lambda b: (sum(b), b))
assert len(pos) == 36
roots = pos + [tuple(-x for x in b) for b in pos]
N = len(roots); idx = {b: k for k, b in enumerate(roots)}
neg = lambda k: idx[tuple(-x for x in roots[k])]
ispos = lambda k: k < 36
S = [tuple(idx[refl(roots[k], i)] for k in range(N)) for i in range(n)]
comp = lambda p, q: tuple(p[q[k]] for k in range(N))       # p o q
ident = tuple(range(N))
def inv(p):
    r = [0]*N
    for k in range(N): r[p[k]] = k
    return tuple(r)
def length(p): return sum(1 for k in range(36) if not ispos(p[k]))
# w_L: longest element of W(D4), D4 = nodes 2,3,4,5 -> indices 1,2,3,4
black = [1, 2, 3, 4]
wL = ident; word = []
while True:
    js = [j for j in black if ispos(wL[idx[simple[j]]])]
    if not js: break
    wL = comp(wL, S[js[0]]); word.append(js[0] + 1)
assert length(wL) == 12
theta = tuple(neg(wL[k]) for k in range(N))           # theta = -w_L on roots
PhiL = [k for k in range(N) if all(roots[k][i] == 0 for i in (0, 5))]
assert len(PhiL) == 24
print("reduced word of w_L (Bourbaki nodes):", word)
print("theta(alpha_1) =", roots[theta[idx[simple[0]]]], " theta(alpha_6) =", roots[theta[idx[simple[5]]]])
# (C1) no real roots: theta(b) != -b ; (C2) b + theta(b) never a root (nor 0)
real = [k for k in range(N) if theta[k] == neg(k)]
print("real roots for theta:", len(real))
sums = [k for k in range(N) if tuple(a+b for a, b in zip(roots[k], roots[theta[k]])) in idx]
print("roots b with b+theta(b) a root:", len(sums))
fixed = [k for k in range(N) if theta[k] == k]
print("theta-fixed (imaginary) roots = Phi_L:", sorted(fixed) == sorted(PhiL))
# conjugacy class of w_L, with representatives w (x = w wL w^-1)
cls = {wL: ident}; fr = [wL]
while fr:
    new = []
    for x in fr:
        for i in range(n):
            y = comp(S[i], comp(x, S[i]))
            if y not in cls:
                cls[y] = comp(S[i], cls[x]); new.append(y)
    fr = new
print("size of conjugacy class of w_L (= number of B-orbits):", len(cls), "; |C_W(w_L)| = 51840/45 =", 51840 // len(cls))
w0 = ident
while True:
    js = [j for j in range(n) if ispos(w0[idx[simple[j]]])]
    if not js: break
    w0 = comp(w0, S[js[0]])
assert length(w0) == 36
print("w_0 in the class:", w0 in cls)
# per orbit data
data = []
for x, w in cls.items():
    assert comp(w, comp(wL, inv(w))) == x
    th = tuple(neg(x[k]) for k in range(N))            # theta_x = -x = w theta w^-1
    assert th == comp(w, comp(theta, inv(w)))
    imag = [k for k in range(N) if th[k] == k]
    assert len(imag) == 24 and all(th[k] != neg(k) for k in range(N))
    assert all(tuple(a+b for a, b in zip(roots[k], roots[th[k]])) not in idx for k in range(N))
    c = sum(1 for k in range(36) if th[k] != k and ispos(th[k])) // 2
    ell = length(x)
    assert ell == 12 + 2*c
    types = []
    for i in range(n):
        a = idx[simple[i]]
        if th[a] == a: types.append('G')
        elif ispos(th[a]): types.append('U+')   # Y is the smaller orbit in P_a Y
        else: types.append('U-')                # Y is open in P_a Y
        # type G iff w^-1 alpha in Phi_L
        assert (types[-1] == 'G') == (inv(w)[a] in PhiL)
    data.append((x, w, c, ell, types))
from collections import Counter
print("orbit dimensions 26-c:", sorted(Counter(26 - d[2] for d in data).items()))
print("closed orbits (theta_x Phi+ = Phi+):", sum(1 for d in data if all(t != 'U-' for t in d[4])), " x = w0 for it:", [d[0] == w0 for d in data if all(t != 'U-' for t in d[4])])
tc = Counter(t for d in data for t in d[4])
print("pairs by type:", dict(tc), " total", sum(tc.values()))
# raising consistency: U- at (x,a) <=> U+ at (s_a x s_a, a), dims differ by one
cdict = {d[0]: d for d in data}
for x, w, c, ell, types in data:
    for i in range(n):
        y = comp(S[i], comp(x, S[i]))
        if types[i] == 'G': assert y == x
        else:
            assert y != x
            t2 = cdict[y][4][i]; c2 = cdict[y][2]
            assert {types[i], t2} == {'U+', 'U-'}
            assert (c2 == c + 1) if types[i] == 'U-' else (c2 == c - 1)
print("raising/lowering consistency: ok")
# open orbit data
d0 = cdict[wL]; print("open orbit: c =", d0[2], "types", d0[4])
# every non-closed orbit has a U- simple root (base case of the induction is the closed orbit)
assert all(('U-' in d[4]) or d[0] == w0 for d in data)
print("every orbit other than the closed one has a simple root of type U- : ok")