#TP(A,B): the tensor product of A and B
TP:=proc(A,B) local i,j,i1,j1,a1,a2,b1,b2,AB,AB1:

#a1 is the number of rows of A
#a2 is the number of columns of A
a1:=nops(A):
a2:=nops(A[1]):


#b1 is the number of rows of B
#b2 is the number of columns of B
b1:=nops(B):
b2:=nops(B[1]):

#TP(A,B) has a1*b1 rows labeled [seq(seq([i,i1],i1=1..b1),i=1..a1)]]:
#it is has a2*b2 columns labeled [seq(seq([j,j1]ji1=1..b2),j=1..a2)]]:
b1:=nops(B):
AB:=[]:
for i from 1 to a1 do
for i1 from 1 to b1 do
AB1:=[]:
for j from 1 to a2 do
for j1 from 1 to b2 do
AB1:=[op(AB1),A[i][j]*B[i1][j1]]:
od:
od:
AB:=[op(AB),AB1]:
od:
od:
end:


#ES(): the four entangled mixed states Sing,T1,T2,T3
ES:=proc()
[[0,1,-1,0]/sqrt(2),
 [0,1,1,0]/sqrt(2),
 [1,0,0,1]/sqrt(2),
 [1,0,0,-1]/sqrt(2)
]:
end: