"""Lemma 23.4 (from the Siegel-Whittaker coefficient to frequency zero), steps 2 and 3, in type B_m, m=2n:
simulate, for r=m-1,...,n+1: row separation (Lemma 22.3 with b=0), removal of e_r+e_m, exchange of U_{e_r} for U_{-e_m},
and conjugation by w_r; check the hypotheses of the exchange lemma at the level of roots (including a+2b),
the root sets, and, with explicit matrices in SO_{2m+1}, that w_r is an isometry of determinant one carrying
the character of step 2 to lambda_{r-1} with all coefficients +1."""
import itertools, sys
import numpy as np
def comb(a,b):
    out=[]
    for i,j in ((1,1),(1,2),(2,1)):
        out.append(tuple(i*x+j*y for x,y in zip(a,b)))
    return out
def is_root(x):
    nz=[t for t in x if t]
    return (len(nz)==2 and all(abs(t)==1 for t in nz)) or (len(nz)==1 and abs(nz[0])==1)
class St:
    def __init__(s,S,chi): s.S=set(S); s.chi=set(chi); s.log=[]; s.grp(s.S,'init')
    def grp(s,S,what):
        for a in S:
            for b in S:
                if a<b:
                    for c in comb(a,b):
                        if is_root(c):
                            if c not in S: s.log.append(f'{what}: not closed {a} {b}')
                            elif c in s.chi: s.log.append(f'{what}: character on commutator {a} {b}')
    def norm(s,U,v,what):
        for u in U:
            for c in comb(u,v):
                if is_root(c):
                    if c not in U: s.log.append(f'{what}: does not normalize ({u})')
                    elif c in s.chi: s.log.append(f'{what}: changes character ({u})')
    def exchange(s,x,y,what):
        assert x in s.S and y not in s.S,(what,x,y)
        U=s.S-{x}; s.grp(U,what+' common')
        z=tuple(p+q for p,q in zip(x,y))
        if not(z in U and z in s.chi): s.log.append(what+': pairing')
        s.norm(U,x,what+' X'); s.norm(U,y,what+' Y'); s.S=U|{y}
    def remove(s,x,what):
        assert x in s.S
        if x in s.chi: s.log.append(what+': removing character root')
        U=s.S-{x}; s.grp(U,what); s.norm(U,x,what+' removed'); s.S=U
def run(n):
    m=2*n
    def v(*t):
        x=[0]*m
        for sg,i in t: x[i-1]+=sg
        return tuple(x)
    pos=set()
    for i in range(1,m+1):
        pos.add(v((1,i)))
        for j in range(i+1,m+1): pos.add(v((1,i),(-1,j))); pos.add(v((1,i),(1,j)))
    def chain(r):
        if r==m-1:
            return pos-{v((1,m))}, {v((1,i),(-1,i+1)) for i in range(1,m)}
        S=pos-{v((1,m))}-{v((1,i),(-1,m)) for i in range(r+1,m)}
        chi={v((1,i),(-1,i+1)) for i in range(1,r)}|{v((1,r),(-1,m))}|{v((1,i),(-1,i+1)) for i in range(r+1,m-1)}|{v((1,m-1),(1,m))}
        return S,chi
    # matrices: basis e_1..e_m,u,f_m..f_1 not needed; use e_1..e_m, f_1..f_m, u
    N=2*m+1; E_=lambda i:i-1; F_=lambda i:m+i-1; Uu=2*m
    def w_matrix(r):
        W=np.zeros((N,N),dtype=int)
        for i in range(1,r): W[E_(i),E_(i)]=1; W[F_(i),F_(i)]=1
        for i in range(r,m-1): W[E_(i+1),E_(i)]=-1; W[F_(i+1),F_(i)]=-1
        W[F_(m),E_(m-1)]=-1; W[E_(m),F_(m-1)]=-1
        W[E_(r),E_(m)]=1; W[F_(r),F_(m)]=1
        W[Uu,Uu]=-1
        return W
    Qm=np.zeros((N,N),dtype=int)   # Gram matrix of B: B(e_i,f_i)=1, B(u,u)=2 (Q(u)=1)
    for i in range(1,m+1): Qm[E_(i),F_(i)]=Qm[F_(i),E_(i)]=1
    Qm[Uu,Uu]=2
    def wroot(r,x):   # action of w_r on weights
        out=[0]*m
        for i in range(1,m+1):
            c=x[i-1]
            if not c: continue
            if i<r: out[i-1]+=c
            elif i<=m-2: out[i]+=c
            elif i==m-1: out[m-1]-=c
            else: out[r-1]+=c
        return tuple(out)
    # root vectors as matrices (Lie algebra of SO(Q), Q = sum x_i y_i + z^2)
    def Emat(a,b):
        M=np.zeros((N,N),dtype=int); M[a,b]=1; return M
    def rootvec(x):
        idx=[(i+1,t) for i,t in enumerate(x) if t]
        if len(idx)==2:
            (i,a),(j,b)=idx
            if a==1 and b==-1: return Emat(E_(i),E_(j))-Emat(F_(j),F_(i))
            if a==1 and b==1: return Emat(E_(i),F_(j))-Emat(E_(j),F_(i))
            if a==-1 and b==1: return (Emat(E_(j),E_(i))-Emat(F_(i),F_(j)))
            if a==-1 and b==-1: return Emat(F_(i),E_(j))-Emat(F_(j),E_(i))
        (i,a),=idx
        if a==1: return 2*Emat(E_(i),Uu)-Emat(Uu,F_(i))   # e_i: u -> 2 e_i, f_i -> -u  (preserves Q)
        return 2*Emat(F_(i),Uu)-Emat(Uu,E_(i))
    for x in list(pos)[:]:
        X=rootvec(x); assert not (X.T.dot(Qm)+Qm.dot(X)).any()
    problems=[]; ok_sets=True; ok_signs=True; ok_iso=True
    S,chi=chain(m-1)
    for r in range(m-1,n,-1):
        st=St(S,chi)
        for j in range(r+1,m):
            st.exchange(v((1,r),(-1,j)),v((1,j),(-1,m)),f'r={r} rowsep j={j}')
        st.remove(v((1,r),(1,m)),f'r={r} remove')
        com=pos-{v((1,m)),v((1,r)),v((1,r),(1,m))}-{v((1,r),(-1,j)) for j in range(r+1,m)}
        if st.S!=com|{v((1,r))}: problems.append(f'r={r}: U_com Z_r mismatch')
        st.exchange(v((1,r)),v((-1,m)),f'r={r} Z<->Y')
        problems+=st.log
        Sn,chin=chain(r-1)
        if {wroot(r,x) for x in Sn}!=st.S: ok_sets=False
        if {wroot(r,x) for x in chin}!=st.chi: ok_sets=False
        W=w_matrix(r)
        if (W.T.dot(Qm).dot(W)-Qm).any() or round(np.linalg.det(W))!=1: ok_iso=False
        Winv=np.linalg.inv(W).round().astype(int)
        # character of step 2 = sum over st.chi of coordinate on root vector; pull back by w_r and compare with lambda_{r-1}
        for x in chin:
            Y=W.dot(rootvec(x)).dot(Winv); tgt=rootvec(wroot(r,x))
            if (Y-tgt).any(): ok_signs=False
        S,chi=Sn,chin
    Sf,chif=chain(n)
    print(f'n={n} (B_{m}): exchange problems={len(problems)}; root sets under w_r match={ok_sets}; '
          f'w_r isometry of det 1={ok_iso}; character coefficients all +1={ok_signs}',problems[:3])
for n in (2,3,4): run(n)