# Lemma on detecting non-nilpotence (label fl:lem:fixed-maximal-detectors), part (2): primes dividing the discriminant of the normalized invariant form
# on the coroot span of every standard Levi, compared with |W_L| and |W_G|.
from sympy import Matrix, factorint
from itertools import combinations
from math import factorial
def cartan(t,n):
    C=[[0]*n for _ in range(n)]
    for i in range(n): C[i][i]=2
    def e(i,j,a=-1,b=-1): C[i][j]=a; C[j][i]=b
    if t=='A':
        for i in range(n-1): e(i,i+1)
    # convention: C[i][j] = <a_i^v, a_j>
    if t=='B':   # alpha_n short: <a_{n-1}^v,a_n>=-1, <a_n^v,a_{n-1}>=-2
        for i in range(n-2): e(i,i+1)
        e(n-2,n-1,-1,-2)
    if t=='C':   # alpha_n long: <a_{n-1}^v,a_n>=-2, <a_n^v,a_{n-1}>=-1
        for i in range(n-2): e(i,i+1)
        e(n-2,n-1,-2,-1)
    if t=='D':
        for i in range(n-2): e(i,i+1)
        e(n-3,n-1)
    if t=='E':
        for a,b in [(0,2),(2,3),(3,4),(4,5),(5,6),(6,7)][:n-2]: e(a,b)
        e(1,3)
    if t=='F': e(0,1); e(1,2,-2,-1); e(2,3)
    if t=='G': e(0,1,-1,-3)
    return Matrix(C)
def symmetrizer(C):
    # positive integers e_j with (C*diag(e)) symmetric, minimal
    n=C.shape[0]; from sympy import Rational, ilcm
    e=[None]*n; e[0]=Rational(1); todo=[0]
    while todo:
        i=todo.pop()
        for j in range(n):
            if C[i,j]!=0 and e[j] is None:
                e[j]=e[i]*C[j,i]/C[i,j]; todo.append(j)   # C_ij e_j = C_ji e_i
    m=min(e); e=[x/m for x in e]
    assert all(x==int(x) for x in e); return [int(x) for x in e]
def weyl_order(t,n):
    return {'A':lambda:factorial(n+1),'B':lambda:2**n*factorial(n),'C':lambda:2**n*factorial(n),
            'D':lambda:2**(n-1)*factorial(n),'E':lambda:{6:51840,7:2903040,8:696729600}[n],
            'F':lambda:1152,'G':lambda:12}[t]()
def components(C,I):
    I=list(I); seen=set(); comps=[]
    for i in I:
        if i in seen: continue
        comp=[i]; seen.add(i); st=[i]
        while st:
            a=st.pop()
            for b in I:
                if b not in seen and C[a,b]!=0: seen.add(b); comp.append(b); st.append(b)
        comps.append(sorted(comp))
    return comps
def comp_weyl(C,comp):
    # order of the Weyl group of an irreducible subdiagram, by brute force on the reflection group
    import itertools
    n=len(comp); sub=C.extract(comp,comp)
    gens=[]
    for i in range(n):
        M=Matrix.eye(n)
        for j in range(n): M[i,j]-=sub[j,i]      # s_i(a_j^v)=a_j^v-<a_j^v,a_i>a_i^v, in coroot coordinates
        gens.append(M)
    seen={tuple(Matrix.eye(n))}; frontier=[Matrix.eye(n)]
    while frontier:
        new=[]
        for g in frontier:
            for s in gens:
                h=s*g; k=tuple(h)
                if k not in seen: seen.add(k); new.append(h)
        frontier=new
    return len(seen)
bad_local=[]; bad_global=[]
for t,ns in [('A',range(1,7)),('B',range(2,7)),('C',range(3,7)),('D',range(4,7)),('E',[6,7,8]),('F',[4]),('G',[2])]:
    for n in ns:
        C=cartan(t,n); e=symmetrizer(C); B=C*Matrix.diag(*e); assert B==B.T
        WG=weyl_order(t,n)
        cache={}
        for k in range(1,n+1):
            for I in combinations(range(n),k):
                d=abs(B.extract(list(I),list(I)).det())
                primes=set(factorint(d).keys())
                WL=1
                for comp in components(C,I):
                    key=tuple(comp)
                    if key not in cache: cache[key]=comp_weyl(C,comp) if len(comp)<=5 else None
                    WL = None if (WL is None or cache[key] is None) else WL*cache[key]
                for p in primes:
                    if WG%p: bad_global.append((t,n,I,p))
                    if WL is not None and WL%p: bad_local.append((t,n,I,p,WL))
print('primes of disc(kappa|Levi) not dividing |W_G| :',bad_global)
print('primes of disc(kappa|Levi) not dividing |W_L| :',bad_local[:10],len(bad_local))