# Checks the embedding of the Jacobi group in SO_5 used in the section on even orthogonal descent
# (basis e1,e2,u,f2,f1; B(e_i,f_j)=delta_ij, B(u,u)=1), and the
# commutator [x_{-e_m}(s), x_{e_r}(t)] = x_{e_r-e_m}(st).
import sympy as sp
x,y,z,x1,y1,z1,b,a,s,t=sp.symbols('x y z x1 y1 z1 b a s t')
idx={'e1':0,'e2':1,'u':2,'f2':3,'f1':4}
def E(p,q):
    M=sp.zeros(5,5); M[idx[p],idx[q]]=1; return M
I=sp.eye(5)
J=sp.zeros(5,5)
for p,q in (('e1','f1'),('e2','f2')): J[idx[p],idx[q]]=1; J[idx[q],idx[p]]=1
J[2,2]=1
xe=lambda i,tt: I+tt*E('e%d'%i,'u')-tt*E('u','f%d'%i)-sp.Rational(1,2)*tt**2*E('e%d'%i,'f%d'%i)
xplus=lambda tt: I+tt*(E('e1','f2')-E('e2','f1'))
xminus=lambda tt: I+tt*(E('e1','e2')-E('f2','f1'))
X=lambda v: xe(2,-2*v); Y=lambda v: xe(1,v); Z=lambda v: xplus(-2*v); n=lambda v: xminus(v/2)
iso=lambda M: sp.simplify(M.T*J*M-J)==sp.zeros(5,5)
g=lambda xx,yy,zz: X(xx)*Y(yy)*Z(zz-xx*yy/2)
disp=sp.Matrix([[1,0,y,x*y-2*z,-y**2/2],[0,1,-2*x,-2*x**2,x*y+2*z],[0,0,1,2*x,-y],[0,0,0,1,0],[0,0,0,0,1]])
print('isometries', all(iso(M) for M in (X(x),Y(y),Z(z),n(b))))
print('displayed matrix', sp.simplify(g(x,y,z)-disp)==sp.zeros(5,5))
print('Heisenberg law', sp.simplify(g(x,y,z)*g(x1,y1,z1)-g(x+x1,y+y1,z+z1+(x*y1-y*x1)/2))==sp.zeros(5,5))
print('n(b) action', sp.simplify(n(b)*g(x,y,z)*n(b).inv()-g(x,y-b*x,z))==sp.zeros(5,5))
low=I+2*a*(E('e2','e1')-E('f1','f2'))
print('lower root action', sp.simplify(low*g(x,y,z)*low.inv()-g(x-a*y,y,z))==sp.zeros(5,5))
# commutator in the ordered five-space (e_r,e_m,u,f_m,f_r): identify e1=e_r, e2=e_m
xneg=lambda ss: I+ss*(E('f2','u')-E('u','e2'))-sp.Rational(1,2)*ss**2*E('f2','e2')
A=xneg(s); B=xe(1,t)
C=sp.simplify(A*B*A.inv()*B.inv())
print('neg root isometry', iso(A), ' commutator =', [sp.simplify(C-xminus(c*s*t))==sp.zeros(5,5) for c in (1,-1)], '(for +st, -st)')