# A6 independent recomputation: E6/F4 Borel orbits via C = class of w_L, theta_x = -x.
import itertools, sys, json
n=6
edges=[(1,3),(3,4),(4,5),(5,6),(2,4)]          # Bourbaki E6
A=[[2 if i==j else 0 for j in range(n)] for i in range(n)]
for i,j in edges: A[i-1][j-1]=A[j-1][i-1]=-1
def ip(u,v): return sum(u[i]*A[i][j]*v[j] for i in range(n) for j in range(n))
roots=[v for v in itertools.product(range(-3,4),repeat=n) if ip(v,v)==2]
assert len(roots)==72 and all(all(c>=0 for c in v) or all(c<=0 for c in v) for v in roots)
idx={v:k for k,v in enumerate(roots)}
pos=[k for k,v in enumerate(roots) if sum(v)>0]; assert len(pos)==36
neg_of=[idx[tuple(-c for c in v)] for v in roots]
ispos=[sum(v)>0 for v in roots]
simple=[idx[tuple(1 if j==i else 0 for j in range(n))] for i in range(n)]
def refl(a,v): c=ip(v,roots[a]); return tuple(v[i]-c*roots[a][i] for i in range(n))
S=[tuple(idx[refl(simple[i],v)] for v in roots) for i in range(n)]   # s_i as permutation
I=tuple(range(72))
def mul(a,b): return tuple(a[b[k]] for k in range(72))   # (a*b)(g)=a(b(g))
def inv(a):
    r=[0]*72
    for k in range(72): r[a[k]]=k
    return tuple(r)
def length(w): return sum(1 for k in pos if not ispos[w[k]])
def neg(w): return tuple(neg_of[w[k]] for k in range(72))   # (-w)
# w_L: longest element of W(D4) on nodes 2,3,4,5 (indices 1,2,3,4)
L=[1,2,3,4]
grp={I:()}; fr=[I]
while fr:
    nf=[]
    for w in fr:
        for i in L:
            u=mul(S[i],w)
            if u not in grp: grp[u]=(i,)+grp[w]; nf.append(u)
    fr=nf
assert len(grp)==192
wL=max(grp,key=length); assert length(wL)==12
for i in L: assert wL[simple[i]]==neg_of[simple[i]]
theta=neg(wL)
a1=roots[theta[simple[0]]]; print('theta(alpha1)=',a1)
PhiL=[k for k in range(72) if theta[k]==k]; print('|fixed roots of -w_L|=',len(PhiL), 'all in D4 span:', all(roots[k][0]==0 and roots[k][5]==0 for k in PhiL))
# conjugacy class C with representatives w (x = w wL w^-1)
C={wL:I}; fr=[wL]
while fr:
    nf=[]
    for x in fr:
        for i in range(n):
            y=mul(mul(S[i],x),S[i])
            if y not in C: C[y]=mul(S[i],C[x]); nf.append(y)
    fr=nf
print('|C|=',len(C))
w0=None
for x in C:
    if all(not ispos[x[k]] for k in pos): w0=x
print('w0 in C:', w0 is not None)
def sum_root(a,b):
    v=tuple(roots[a][i]+roots[b][i] for i in range(n)); return v
table={}
bad=[]
for x,w in C.items():
    th=neg(x)
    assert mul(mul(w,wL),inv(w))==x
    Phx=[k for k in range(72) if th[k]==k]
    Phx_w=sorted(w[k] for k in PhiL)
    if sorted(Phx)!=Phx_w: bad.append('PhiX')
    real=[k for k in range(72) if th[k]==neg_of[k]]
    if real: bad.append('real')
    for k in range(72):
        v=sum_root(k,th[k])
        if v in idx or all(c==0 for c in v): bad.append('sum')
    pairs=set()
    for k in pos:
        if th[k]!=k and ispos[th[k]]: pairs.add(frozenset((k,th[k])))
    c=len(pairs); l=length(x)
    if l!=12+2*c: bad.append('len')
    types=[]; partners=[]
    for i in range(n):
        a=simple[i]
        xa=x[a]
        if a in Phx: t='G'
        elif ispos[xa]: t='Uo'
        else: t='Us'
        types.append(t)
        partners.append(mul(mul(S[i],x),S[i]))
        if (t=='G')!=(partners[-1]==x): bad.append('partner')
    table[x]=dict(c=c,dim=26-c,len=l,types=types,partners=partners,Phx=set(Phx))
print('bad checks:',sorted(set(bad)))
from collections import Counter
print('dim multiplicities:',sorted(Counter(t['dim'] for t in table.values()).items()))
print('type counts:',Counter(tt for t in table.values() for tt in t['types']))
print('open orbit x=wL types:',table[wL]['types'],'c=',table[wL]['c'])
print('w0 dim:',table[w0]['dim'])
# partner dimension changes
ok=True
for x,t in table.items():
    for i in range(n):
        p=t['partners'][i]; tp=table[p]
        if t['types'][i]=='Uo' and not (tp['dim']==t['dim']-1 and tp['types'][i]=='Us'): ok=False
        if t['types'][i]=='Us' and not (tp['dim']==t['dim']+1 and tp['types'][i]=='Uo'): ok=False
print('partner dims consistent:',ok)
import pickle
pickle.dump(dict(roots=roots,C=C,wL=wL,w0=w0,table={x:dict(c=t['c'],dim=t['dim'],types=t['types'],partners=t['partners'],Phx=t['Phx']) for x,t in table.items()},S=S),open(sys.argv[1],'wb'))