# A6: T^theta over Z (cokernel of theta^*-1 on X^*(T)=weight lattice), root restrictions, C_W(w_L).
import pickle, sys, itertools
from math import gcd
D=pickle.load(open(sys.argv[1],'rb'))
roots=D['roots']; C=D['C']; S=D['S']; wL=D['wL']
n=6
edges=[(1,3),(3,4),(4,5),(5,6),(2,4)]
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
idx={v:k for k,v in enumerate(roots)}
simple=[idx[tuple(1 if j==i else 0 for j in range(n))] for i in range(n)]
# w_L as matrix on the weight lattice (fundamental-weight coordinates): s_i(c)_j = c_j - c_i*A[i][j]
def s_w(i,c): return tuple(c[j]-c[i]*A[i][j] for j in range(n))
# reduced word of wL: peel left descents using the permutation rep
def mul(a,b): return tuple(a[b[k]] for k in range(72))
def inv(a):
    r=[0]*72
    for k in range(72): r[a[k]]=k
    return tuple(r)
ispos=[sum(v)>0 for v in roots]
word=[]; w=wL
I=tuple(range(72))
while w!=I:
    wi=inv(w)
    for i in range(n):
        if not ispos[wi[simple[i]]]:
            word.append(i); w=mul(S[i],w); break
print('reduced word of w_L (1-based):',[i+1 for i in word])
def wL_w(c):
    for i in reversed(word): c=s_w(i,c)
    return c
E=[tuple(1 if j==i else 0 for j in range(n)) for i in range(n)]
M=[wL_w(e) for e in E]            # M[i] = image of omega_i
th=[tuple(-c for c in M[i]) for i in range(n)]   # theta^* = -w_L on weights
Dm=[[th[i][j]-(1 if i==j else 0) for i in range(n)] for j in range(n)]  # column i = (theta^*-1)(omega_i)
def minors(Mx,k):
    g=0
    rows=range(len(Mx)); cols=range(len(Mx[0]))
    for R in itertools.combinations(rows,k):
        for Cc in itertools.combinations(cols,k):
            sub=[[Mx[r][c] for c in Cc] for r in R]
            # det
            def det(m):
                if len(m)==1: return m[0][0]
                return sum((-1)**j*m[0][j]*det([row[:j]+row[j+1:] for row in m[1:]]) for j in range(len(m)))
            g=gcd(g,abs(det(sub)))
    return g
d1=minors(Dm,1); d12=minors(Dm,2); d3=minors(Dm,3)
print('gcd 1x1 minors',d1,' gcd 2x2 minors',d12,' gcd 3x3 minors',d3,'-> elementary divisors',d1,d12//d1 if d1 else None)
# root restrictions to T^theta: alpha in m X + im(theta^*-1) ?  (alpha in weight coords: alpha_i = sum_j A[i][j] omega_j)
def root_w(v): return tuple(sum(v[i]*A[i][j] for i in range(n)) for j in range(n))
def rank_mod(rows,p):
    m=[[x%p for x in r] for r in rows]; rk=0; cols=len(m[0]) if m else 0
    for c in range(cols):
        piv=None
        for r in range(rk,len(m)):
            if m[r][c]%p: piv=r;break
        if piv is None: continue
        m[rk],m[piv]=m[piv],m[rk]
        invp=pow(m[rk][c],p-2,p)
        m[rk]=[(x*invp)%p for x in m[rk]]
        for r in range(len(m)):
            if r!=rk and m[r][c]: m[r]=[(a-m[r][c]*b)%p for a,b in zip(m[r],m[rk])]
        rk+=1
    return rk
imgD=[[Dm[j][i] for j in range(n)] for i in range(n)]  # rows = images
div={}
for p in [q for q in range(2,98) if all(q%r for r in range(2,q))]:
    r0=rank_mod(imgD,p)
    div[p]=sum(1 for v in roots if rank_mod(imgD+[list(root_w(v))],p)==r0)
print('number of roots whose restriction to T^theta is divisible by p:',div)
# C_W(w_L): generators s2..s5, s1 s_beta1, s6 s_beta6, beta_i = -theta alpha_i
theta=tuple(roots.index(tuple(-c for c in roots[wL[k]])) for k in range(72))
def refl_perm(a):
    def ipv(u,v): return sum(u[i]*A[i][j]*v[j] for i in range(n) for j in range(n))
    return tuple(idx[tuple(roots[k][i]-ipv(roots[k],roots[a])*roots[a][i] for i in range(n))] for k in range(72))
negi=[idx[tuple(-c for c in v)] for v in roots]
gens=[S[i] for i in (1,2,3,4)]
for i in (0,5):
    b=negi[theta[simple[i]]]
    print('beta_%d ='%(i+1),roots[b],' (alpha,beta)=',sum(roots[simple[i]][a]*A[a][c]*roots[b][c] for a in range(n) for c in range(n)),
          ' alpha+beta root?',tuple(x+y for x,y in zip(roots[simple[i]],roots[b])) in idx,' alpha-beta root?',tuple(x-y for x,y in zip(roots[simple[i]],roots[b])) in idx)
    gens.append(mul(S[i],refl_perm(b)))
grp={I}; fr=[I]
while fr:
    nf=[]
    for w in fr:
        for g in gens:
            u=mul(g,w)
            if u not in grp: grp.add(u); nf.append(u)
    fr=nf
print('|<gens>|=',len(grp),' all commute with w_L:',all(mul(g,wL)==mul(wL,g) for g in grp), ' 51840/45 =',51840//45)