# Verify the characteristic-two matrix identities used in the even orthogonal descent.
from sympy import symbols, Matrix, eye, zeros, Poly, expand
s,a,v1,v2,v3,w1,w2,w3,t=symbols('s a v1 v2 v3 w1 w2 w3 t')
gens=(s,a,v1,v2,v3,w1,w2,w3,t)
def red(M):
    return M.applyfunc(lambda x: Poly(expand(x),*gens,modulus=2).as_expr())
def eq(A,B): return red(A-B)==zeros(*A.shape)
def inv(M): return red(M.inv())
# --- check 1: basis (e_r,e_m,u,f_m,f_r) = 0..4 ; E(a,b) sends basis b to a
def E(n,i,j):
    M=zeros(n,n); M[i,j]=1; return M
er,em,u,fm,fr=range(5)
x_er_em=lambda c: eye(5)+c*(E(5,er,fm)+E(5,em,fr))
x_mem=lambda c: eye(5)+c*E(5,u,em)+c**2*E(5,fm,em)
x_er=lambda c: eye(5)+c*E(5,u,fr)+c**2*E(5,er,fr)
x_er_mem=lambda c: eye(5)+c*(E(5,er,em)+E(5,fm,fr))
A=x_er_em(s); B=x_mem(a)
comm=red(A*B*inv(A)*inv(B))
print('commutator (21) :',eq(comm,x_er(s*a)*x_er_mem(s*a**2)))
# isometry check for Q=z^2+x_r y_r+x_m y_m
def Q(vec): return vec[u]**2+vec[er]*vec[fr]+vec[em]*vec[fm]
c=symbols('c0:5'); vec=Matrix(c)
for name,M in [('x_{e_r+e_m}',A),('x_{-e_m}',B),('x_{e_r}',x_er(s)),('x_{e_r-e_m}',x_er_mem(s))]:
    d=Poly(expand(Q(M*vec)-Q(vec)),*c,*gens,modulus=2)
    print('isometry',name,d.is_zero)
# --- checks 2,3: basis e1,e2,e,u,f,f2,f1 ; terminal space <e,u,f>, Q=z^2+xy
e1,e2,e,uu,f,f2,f1=range(7)
dual={e1:f1,e2:f2}
def Qt(v): return v[1]**2+v[0]*v[2]          # v=(x,z,y) on (e,u,f)
def Bt(v,w): return v[0]*w[2]+v[2]*w[0]       # polar form (char 2)
def r(i,v):
    M=eye(7); fi=dual[i]
    M[e,fi]+=v[0]; M[uu,fi]+=v[1]; M[f,fi]+=v[2]; M[i,fi]+=Qt(v)
    # r_i(v) w = w + B(w,v) e_i for terminal w
    M[i,e]+=Bt((1,0,0),v); M[i,uu]+=Bt((0,1,0),v); M[i,f]+=Bt((0,0,1),v)
    return M
def xplus(c): return eye(7)+c*(E(7,e1,f2)+E(7,e2,f1))
def xdiff(c): return eye(7)+c*(E(7,e1,e2)+E(7,f2,f1))   # x_{e1-e2}
V=(v1,v2,v3); W=(w1,w2,w3)
print('additive r_i     :',eq(r(e1,V)*r(e1,W),r(e1,tuple(p+q for p,q in zip(V,W)))))
C=red(r(e1,V)*r(e2,W)*inv(r(e1,V))*inv(r(e2,W)))
print('[r_i(v),r_j(w)]  :',eq(C,xplus(Bt(V,W))))
L=red(xdiff(s)*r(e2,V)*inv(xdiff(s)))
print('conjugation      :',eq(L,r(e2,V)*r(e1,tuple(s*p for p in V))*xplus(s*Qt(V))))
U=(0,a,0)
C3=red(xdiff(s)*r(e2,U)*inv(xdiff(s))*inv(r(e2,U)))
print('[X_j(s),Y_j(a)]  :',eq(C3,r(e1,(0,s*a,0))*xplus(s*a**2)))
# --- split basis of Corollary: c0=s^2+s
c0=s**2+s
def lam(v): return v[1]+v[0]+c0*v[2]
es=(s**2,s,1); fs=(s**2+1,s+1,1)
P=lambda x: Poly(expand(x),s,modulus=2).as_expr()
print('split basis      :',P(Qt(es)),P(Qt(fs)),P(lam(es)),P(lam(fs)),P(Bt(es,fs)))
# --- pair matrix
c=symbols('c')
g=Matrix([[c,1],[1,0]]); Ct=Matrix([[1,0,0],[0,1,1]])
print('g^t C            :',(g.T*Ct).tolist(),' g*(1,0)=',(g*Matrix([1,0])).T.tolist())