# Independent audit of A8-4/A8-5/A8-6 (E6/F4).  Bourbaki labels 1-3-4-5-6, node 2 on node 4.
# Roots: simple-root coords; coweights: coroot coords; <alpha,nu> = alpha^T A nu (simply laced).
from fractions import Fraction as Fr
import random
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
pair=lambda a,nu: sum(a[i]*A[i][j]*nu[j] for i in range(n) for j in range(n))
def refl(j,v):  # same formula for roots and coweights (simply laced)
    c=sum(A[j][i]*v[i] for i in range(n)); w=list(v); w[j]-=c; return tuple(w)
def orbit(v,gens=range(n)):
    seen={v}; fr=[v]
    while fr:
        nf=[]
        for x in fr:
            for j in gens:
                y=refl(j,x)
                if y not in seen: seen.add(y); nf.append(y)
        fr=nf
    return seen
E=[tuple(1 if i==j else 0 for i in range(n)) for j in range(n)]
Phi=set()
for e in E: Phi|=orbit(e)
assert len(Phi)==72
TH=[[-1,0,0,0,0,0],[-1,1,0,0,0,-1],[-2,0,1,0,0,-1],[-2,0,0,1,0,-2],[-1,0,0,0,1,-2],[0,0,0,0,0,-1]]
th=lambda v: tuple(sum(TH[i][j]*v[j] for j in range(n)) for i in range(n))
phi=lambda v: tuple(v[i]-th(v)[i] for i in range(n))
ab=lambda v: (pair(E[0],phi(v)), pair(E[5],phi(v)))   # (a,b) of phi(v)
# theta is an isometry and maps roots to roots with the same matrix
assert all(th(a) in Phi for a in Phi)
assert all(pair(th(a),th(v))==pair(a,v) for a in Phi for v in E)
black=[1,2,3,4]
# (1) no root vanishes on T' = (T^theta)^0 = torus of black coroots
van=[a for a in Phi if all(pair(a,E[j])==0 for j in black)]
print('(1) roots vanishing on the black torus T\':',len(van))
while True:   # generic cocharacter eta of T'
    eta=tuple(random.randint(-50,50) if j in black else 0 for j in range(n))
    if all(pair(a,eta)!=0 for a in Phi): break
print('    generic eta in X_*(T\') with <alpha,eta>!=0 for all 72 roots:',eta)
w1=(Fr(4,3),Fr(3,3),Fr(5,3),Fr(6,3),Fr(4,3),Fr(2,3)); w6=(Fr(2,3),Fr(3,3),Fr(4,3),Fr(6,3),Fr(5,3),Fr(4,3))
two_rho=(16,22,30,42,30,16)
PhiL=[a for a in Phi if a[0]==0 and a[5]==0]; assert len(PhiL)==24
assert all(th(a)==a for a in PhiL)
def audit(k,mirror=False):
    c=Fr(3*k-1,2)
    lam=tuple(int(c*w1[i]+w6[i]) if not mirror else int(w1[i]+c*w6[i]) for i in range(n))
    assert all((c*w1[i]+w6[i] if not mirror else w1[i]+c*w6[i])==lam[i] for i in range(n))
    mu=(3*k,0) if not mirror else (0,3*k)
    mup=(3*k-1,2) if not mirror else (2,3*k-1)
    assert ab(lam)==mup and all(pair(E[j],lam)>=0 for j in range(n))
    D=sum(two_rho[i]*A[i][j]*lam[j] for i in range(n) for j in range(n))//2  # <2rho,lam>
    O=orbit(lam)
    Fk=[v for v in O if ab(v)==mu]; Fkp=[v for v in O if ab(v)==mup]
    # Levi of mu: roots with zero alpha_1 (resp. alpha_6) coefficient -> D5
    idx=0 if not mirror else 5
    L5=[a for a in Phi if a[idx]==0]; assert len(L5)==40
    out5=[a for a in L5 if a not in PhiL]; assert len(out5)==16 and all(th(a)!=a and th(a) in out5 for a in out5)
    pairs={frozenset((a,th(a))) for a in out5}; assert len(pairs)==8
    dimZ=4+len(PhiL)+len(pairs)
    counts=[]
    for nu in Fk:
        nL=sum(1 for a in PhiL if pair(a,nu)<=0)
        npr=sum(1 for p in pairs if all(pair(a,nu)<=0 for a in p))
        counts.append((nL,npr,4+nL+npr))
    # rho_L pairing
    rhoL2=tuple(sum(a[i] for a in PhiL if all(x>=0 for x in a)) for i in range(n))
    rl=sorted(Fr(pair(rhoL2,nu),2) for nu in Fk)
    print('k=%d%s: lam=%s |W lam|=%d <2rho,lam>=%d delta(mu\')=8(a+b)=%d; #F^T\'=%d over mu, #F\'^T\'=%d over mu\'; dim Z=%d; dim(LieZ cap LieP) per fixed pt=%s; <rho_L,nu>=%s'%(
        k,' (mirror)' if mirror else '',lam,len(O),D,8*sum(mup),len(Fk),len(Fkp),dimZ,sorted(set(counts)),[str(x) for x in rl]))
    # A8-4 step 3: no dominant lam'' = lam - kappa, kappa in nonneg black cone, kappa != 0 (box check)
    bad=0
    for kap in __import__('itertools').product(range(7),repeat=4):
        if sum(kap)==0: continue
        v=list(lam)
        for j,cj in zip(black,kap): v[j]-=cj
        if all(pair(E[j],tuple(v))>=0 for j in range(n)): bad+=1
    print('    A8-4 step 3 box check (c_j<=6): dominant lam\'\' found =',bad)
for k in (1,3,5): audit(k)
for k in (1,3): audit(k,mirror=True)