# Exact lattice checks for E6/F4 (Bourbaki labels 1..6; chain 1-3-4-5-6, node 2 on node 4).
# Coroot lattice Q^v = X_*(T) for simply connected E6. theta = -w_L on X_*(T), w_L = longest of W(D4) on nodes 2,3,4,5.
from fractions import Fraction
from itertools import product
edges = {(1,3),(3,4),(4,5),(5,6),(2,4)}
def a(i,j):  # Cartan matrix <alpha_i^vee, alpha_j>  (simply laced, symmetric)
    if i==j: return 2
    return -1 if (min(i,j),max(i,j)) in edges else 0
N=[1,2,3,4,5,6]
C=[[a(i,j) for j in N] for i in N]
def refl(v,i):  # simple reflection s_i on coweights written in coroot coordinates v (list of 6)
    # <alpha_i, v> = sum_j v_j <alpha_j^vee, alpha_i> = sum_j v_j C[j][i]
    p=sum(v[j]*C[j][i] for j in range(6))
    w=list(v); w[i]-=p; return w
# build w_L as a word: longest element of D4 on nodes 2,3,4,5 (indices 1,2,3,4); find by BFS the element mapping D4-positive coroots to negative
D4=[1,2,3,4]
def apply_word(v,word):
    for i in reversed(word): v=refl(v,i)
    return v
# generate W(D4) elements as words via BFS on action on a generic vector
gen=[Fraction(1,1)*x for x in [1,10,100,1000,10000,100000]]
start=tuple(gen); seen={start:[]}; frontier=[start]
while frontier:
    nf=[]
    for v in frontier:
        for i in D4:
            w=tuple(refl(list(v),i))
            if w not in seen:
                seen[w]=seen[v]+[i]; nf.append(w)
    frontier=nf
print("|W(D4)| =",len(seen))
# w_L = the element sending each D4 simple coroot to its negative (w_L = -1 on D4 span)
def wL(v):
    return None
cands=[]
for img,word in seen.items():
    ok=True
    for i in D4:
        e=[0]*6; e[i]=1
        if apply_word(e,word)!=[-x for x in e]: ok=False;break
    if ok: cands.append(word)
print("w_L found:",len(cands)==1, "length",len(cands[0]))
WL=cands[0]
def theta(v): return [-x for x in apply_word(list(v),WL)]
# fundamental coweights omega_i^vee in coroot coordinates: solve C^T x = e_i  (<alpha_j, omega_i^vee> = delta_ij)
def solve(M,b):
    n=len(M); A=[[Fraction(M[r][c]) for c in range(n)]+[Fraction(b[r])] for r 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]
        for r in range(n):
            if r!=c and A[r][c]!=0:
                f=A[r][c]/A[c][c]; A[r]=[A[r][k]-f*A[c][k] for k in range(n+1)]
    return [A[r][n]/A[r][r] for r in range(n)]
CT=[[C[j][i] for j in range(6)] for i in range(6)]  # row i: <alpha_i, alpha_j^vee>... pairing <alpha_i, v> = sum_j v_j C[j][i]
om={i:solve(CT,[1 if k==i else 0 for k in range(6)]) for i in range(6)}
print("omega1^v =",om[0]); print("omega6^v =",om[5])
e=lambda i: [x-y for x,y in zip([1 if k==i else 0 for k in range(6)], theta([1 if k==i else 0 for k in range(6)]))]
e1,e6=e(0),e(5)
print("e1 =",e1,"  e6 =",e6)
# check e1 = 2 om1 - om6, e6 = -om1 + 2 om6
print("e1 == 2w1-w6:", e1==[2*x-y for x,y in zip(om[0],om[5])], "  e6 == -w1+2w6:", e6==[-x+2*y for x,y in zip(om[0],om[5])])
# phi(alpha_j^vee) for j in D4 is 0
print("phi(alpha_j^v)=0 for j=2..5:", all(e(j)==[0]*6 for j in D4))
# X_*(S) = {v in Q^v : theta v = -v}: integer solutions in a box, compare with Z e1 + Z e6
def in_span_e1e6(v):
    # solve v = x e1 + y e6 using pairings with alpha_1, alpha_6
    p1=sum(v[j]*C[j][0] for j in range(6)); p6=sum(v[j]*C[j][5] for j in range(6))
    # <alpha1,e1>=2? compute
    q11=sum(e1[j]*C[j][0] for j in range(6)); q16=sum(e6[j]*C[j][0] for j in range(6))
    q61=sum(e1[j]*C[j][5] for j in range(6)); q66=sum(e6[j]*C[j][5] for j in range(6))
    det=q11*q66-q16*q61
    x=Fraction(p1*q66-q16*p6,det); y=Fraction(q11*p6-p1*q61,det)
    return x.denominator==1 and y.denominator==1 and [x*a_+y*b_ for a_,b_ in zip(e1,e6)]==list(v)
R=range(-4,5); bad=0; cnt=0
for v in product(R,repeat=6):
    v=list(v)
    if theta(v)==[-x for x in v]:
        cnt+=1
        if not in_span_e1e6(v): bad+=1
print("anti-invariant coroot vectors in box:",cnt," not in Z e1+Z e6:",bad)
# rank of X_*(T)^theta (should be 4 = rank F4)
inv=[list(v) for v in product(range(-2,3),repeat=6) if theta(list(v))==list(v)]
def rank(vs):
    M=[[Fraction(x) for x in v] for v in vs]; r=0
    for c in range(6):
        p=next((i for i in range(r,len(M)) if M[i][c]!=0),None)
        if p is None: continue
        M[r],M[p]=M[p],M[r]
        for i in range(len(M)):
            if i!=r and M[i][c]!=0:
                f=M[i][c]/M[r][c]; M[i]=[M[i][k]-f*M[r][k] for k in range(6)]
        r+=1
    return r
print("rank X_*(T)^theta =",rank(inv))
# 3 omega1^v = 2 e1 + e6 ?
print("3w1 == 2e1+e6:", [3*x for x in om[0]]==[2*x+y for x,y in zip(e1,e6)])