> with(linalg): with(group):

Warning, the protected names norm and trace have been redefined and unprotected

> kr:=proc(A,B)
# This procedure computes the Kronecker product of two matrices
local mm,nn,ans,ans1;
ans:=[];
for mm from 1 to rowdim(A) do
for nn from 1 to coldim(A) do
ans:=[op(ans), scalarmul(B, A[mm,nn] )];
od;
od;
ans1:=blockmatrix(rowdim(A),coldim(A), ans );
RETURN(evalm(ans1));
end:

> sigma||0:=diag(1,1):
sigma||1:=matrix(2,2,[0,1,1,0]):
sigma||2:=matrix(2,2,[0,-I,I,0]):
sigma||3:=matrix(2,2,[1,0,0,-1]):
#

> for i from 0 to 3 do for j from 0 to 3 do
SS||(i+4*j):=kr(sigma||i, sigma||j)
od; od;
#

> for i from 0 to 3 do
SSS||i:=evalm(sigma||i);
od:
#

> for i from 0 to 3 do for j from 0 to 3 do
SSS||i||j:=evalm(SS||(i+4*j));
od; od;

> for i from 0 to 3 do for j from 0 to 3 do
for k from 0 to 3 do
SSS||i||j||k:=kr(SS||(i+4*j), sigma||k):
od;
od; od;
#

> for i from 0 to 3 do for j from 0 to 3 do
ij:=i+4*j;
for k from 0 to 3 do for m from 0 to 3 do
SSS||i||j||k||m:=kr(SS||(ij), SS||(k+4*m)):
od; od;
od; od;

> L3x3||0:=matrix(3,3,[1,0,0, 0,1,0, 0,0,1]):
L3x3||1:=matrix(3,3, [0,0,0, 0,0,1, 0,-1,0]):
L3x3||2:=matrix(3,3, [0,0,-1, 0,0,0, 1,0,0]):
L3x3||3:=matrix(3,3, [0,1,0, -1,0,0, 0,0,0]):
for i from 0 to 3 do
LLL||i:=evalm(L3x3||i);
od:
for i from 0 to 3 do for j from 0 to 3 do
LLL||i||j:=kr(LLL||i, LLL||j);
od; od;
#

> LL||1:=matrix(4,4,[0,1,0,0, -1,0,0,0, 0,0,0,0, 0,0,0,0]):
LL||2:=matrix(4,4,[0,0,1,0, 0,0,0,0, -1,0,0,0, 0,0,0,0]):
LL||3:=matrix(4,4,[0,0,0,1, 0,0,0,0, 0,0,0,0, -1,0,0,0]):
LL||4:=matrix(4,4,[0,0,0,0, 0,0,1,0, 0,-1,0,0, 0,0,0,0]):
LL||5:=matrix(4,4,[0,0,0,0, 0,0,0,1, 0,0,0,0, 0,-1,0,0]):
LL||6:=matrix(4,4,[0,0,0,0, 0,0,0,0, 0,0,0,1, 0,0,-1,0]):
#

> ket||0:=matrix(2,1,[1,0]):
ket||1:=matrix(2,1,[0,1]):
#
for i from 0 to 1 do for j from 0 to 1 do
ket||i||j:=kr(ket||i,ket||j);
od; od;
#
for i from 0 to 1 do for j from 0 to 1 do
for k from 0 to 1 do
ket||i||j||k:=kr(ket||i||j, ket||k);
od;
od; od;
#
for i from 0 to 1 do for j from 0 to 1 do for k from 0 to 1 do
for m from 0 to 1 do
ket||i||j||k||m:=kr(ket||i||j||k, ket||m);
od;
od; od; od;

> Ket_To_Rho:=v->evalm(v &* map(conjugate, transpose(v))):

> Bell||0||0:=evalm( 1/sqrt(2) * evalm( ket||0||0 + ket||1||1 ) ):
Bell||0||1:=evalm( 1/sqrt(2) * evalm( ket||0||1 + ket||1||0 ) ):
Bell||1||0:=evalm( 1/sqrt(2) * evalm( ket||0||1 - ket||1||0 ) ):
Bell||1||1:=evalm( 1/sqrt(2) * evalm( ket||0||0 - ket||1||1 ) ):
#
Bell||0||0||0:=evalm( 1/sqrt(2) * evalm( ket||0||0||0 + ket||1||1||1 ) ):
Bell||0||0||1:=evalm( 1/sqrt(2) * evalm( ket||0||0||1 + ket||1||1||0 ) ):
Bell||0||1||0:=evalm( 1/sqrt(2) * evalm( ket||0||1||0 + ket||1||0||1 ) ):
Bell||0||1||1:=evalm( 1/sqrt(2) * evalm( ket||1||0||0 - ket||0||1||1 ) ):

Bell||1||0||0:=evalm( 1/sqrt(2) * evalm( ket||1||0||0 + ket||0||1||1 ) ):
Bell||1||0||1:=evalm( 1/sqrt(2) * evalm( ket||0||1||0 - ket||1||0||1 ) ):
Bell||1||1||0:=evalm( 1/sqrt(2) * evalm( ket||0||0||1 - ket||1||1||0 ) ):
Bell||1||1||1:=evalm( 1/sqrt(2) * evalm( ket||0||0||0 - ket||1||1||1 ) ):
#
Sam||0||0||0:=evalm( 1/2 * evalm( ket||1||1||0 + ket||1||0||1 + ket||0||1||1 - ket||0||0||0 ) ):
Sam||0||0||1:=evalm( 1/2 * evalm( ket||1||0||0 + ket||0||1||0 - ket||0||0||1 + ket||1||1||1 ) ):
Sam||0||1||0:=evalm( 1/2 * evalm( ket||1||0||0 - ket||0||1||0 + ket||0||0||1 + ket||1||1||1 ) ):
Sam||0||1||1:=evalm( 1/2 * evalm( ket||1||1||0 + ket||1||0||1 - ket||0||1||1 + ket||0||0||0 ) ):

Sam||1||0||0:=evalm( 1/2 * evalm( -ket||1||0||0 + ket||0||1||0 + ket||0||0||1 + ket||1||1||1 ) ):
Sam||1||0||1:=evalm( 1/2 * evalm( ket||1||1||0 - ket||1||0||1 + ket||0||1||1 + ket||0||0||0 ) ):
Sam||1||1||0:=evalm( 1/2 * evalm( -ket||1||1||0 + ket||1||0||1 + ket||0||1||1 + ket||0||0||0 ) ):
Sam||1||1||1:=evalm( 1/2 * evalm( ket||1||0||0 + ket||0||1||0 + ket||0||0||1 - ket||1||1||1 ) ):

> Bell4:=proc(a,b,c,d)
# Bell4(a,b,c,d) = ( ket.(1-a).(1-b).(1-c).(1-d) + ((-1)^a)*ket.a.b.c.d)/sqrt(2)
# where a,b,c,d are elts of {0,1}
#
local BellKet;
#
BellKet:=evalm( ket||(1-a)||(1-b)||(1-c)||(1-d) );
if (a=0) then
BellKet:=evalm( BellKet + ket||a||b||c||d );
else
BellKet:=evalm( BellKet - ket||a||b||c||d );
fi;
BellKet:=evalm( (1/sqrt(2)) * BellKet );
RETURN(evalm(BellKet));
end:
for a from 0 to 1 do for b from 0 to 1 do
for c from 0 to 1 do for d from 0 to 1 do
Bell||a||b||c||d:=Bell4(a,b,c,d);
od; od;
od; od;
a:='a':b:='b':c:='c':d:='d':

> QF:=proc(U)
# This procedure computes a unitary matrix Q whose rows are
# are the unit eigenvectors of a unitary matrix U
local Dim_U,List,n,nu,Q,j,temp_j,n_j,List_j,List_jj,i,Vec_i,Norm_Vec_i,k;

Dim_U:=rowdim(U):
List:=[eigenvects(U)];
n:=nops(List);
nu:=1;

Q:=matrix(Dim_U,Dim_U);
for j from 1 to n do
temp_j:=op(j,List);
n_j:=op(2,temp_j);
List_j:= convert(op(3,temp_j), list);
List_jj:=GramSchmidt(List_j);

for i from 1 to n_j do
Vec_i:=convert( op(i, List_jj), vector );
Norm_Vec_i:=norm(eval(Vec_i),2);

for k from 1 to Dim_U do
Q[nu,k]:=Vec_i[k]/Norm_Vec_i;
od;
nu:=nu+1;
od;
od;
Q:=map(radnormal, Q);
RETURN(evalm(Q))
end:

> Logz:=proc(z)
local zpolar;
if (z=0) then
RETURN(0);
else
zpolar:=convert(z,polar);
RETURN( (I*op(2,zpolar)) + ln(op(1,zpolar)) );
fi;
end:

> LogU:=proc(U)
# This procedure computes the natural log of a unitary matrix
local Q,IQ,DU,LDU,LU;

Q:=QF(U);
IQ:=map(conjugate,transpose(Q));
DU:=map(radnormal,evalm( Q &* U &* IQ));
LDU:=map(Logz, DU);
LU:=map(radnormal, evalm( IQ &* LDU &* Q));
RETURN(evalm(LU));
end:

> Singer2Log:=proc(W)
# Computes multiple valued log of matrix W
local LW,UQ,IUQ,DD,n1,n2,n3,n4;

LW:=LogU(W):
UQ:=QF(W):
IUQ:=map(conjugate,transpose(UQ)):
DD:=map(radnormal, evalm( IUQ &* diag(n1,n2,n3,n4) &* UQ)):
RETURN(evalm(LW + DD));
end:

> Pauli1:=proc(LU)
# This procedure expresses a 2x2 skew-Hermitian matrix as
# the sum of Pauli matrices.
# LU:= Sum(i=0..3, R[i+1]*sigma.i )
# Please note that the index of the matrix R runs from 1 to 4
local A, i;

A:=matrix(1,4);
for i from 0 to 3 do
A[1,i+1]:=trace( evalm( sigma||i &* LU) )/2;
od;
RETURN(evalm(A));
end:

> RPauli1:=proc(LU)
local R,RSUM,i;
#
# Expresses density operator rho as a formal sum of sigma.i's
R:=Pauli1(LU);
RSUM:=0;
for i from 0 to 3 do
RSUM:=RSUM + R[1,i+1]*S||i;
od;
RETURN(RSUM);
end:

> Pauli2:=proc(LU)
# This procedure expresses a 4x4 skew-Hermitian matrix as
# the sum of the tensor product of Pauli matrices.
# LU:= Sum( i=0..3, j=0..3, R[i+1,j+1]*SS.(i+4*j) )
# Please note that the indices of the matrix R run from 1 to 4
local A, i, j;

A:=matrix(rowdim(LU),coldim(LU));
for i from 0 to 3 do
for j from 0 to 3 do
A[i+1,j+1]:=trace( evalm( SS||(i+4*j) &* LU) )/4;
od;
od;
RETURN(evalm(A));
end:

> RPauli2:=proc(rho)
local R,RSUM,i,j;
#
# Expresses density operator rho as a formal sum of sigma.i(X)sigma.j's
#
R:=Pauli2(rho);
RSUM:=0;
for i from 0 to 3 do for j from 0 to 3 do
RSUM:=RSUM + R[i+1,j+1]*S||i||j;
od; od;
RETURN(RSUM);
end:

> Inv_Pauli2:=proc(R)
# Given a 4x4 matrix R over the reals, this procedure
# computes rho = Sum(p=0..3, q=0..3) R[p+1,q+q]*sigma.p(X)sigma.q
#
local p,q,rho;
#
rho:=diag(0,0,0,0):
for p from 0 to 3 do for q from 0 to 3 do
rho:= evalm( rho + evalm( R[p+1,q+1]*SS||(p+4*q) ) );
od; od;
RETURN(evalm(rho));
end:

> Pauli3:=proc(LU)
# This procedure expresses a 8x8 skew-Hermitian matrix as
# the sum of the tensor product of Pauli matrices.
# LU:= Sum( i=0..3, j=0..3, k=0..3, R[i,j,k]*SSS.i.j.k )
# Please note that the indices of the matrix R run from 0 to 3
local A, i, j, k;

A:=array(0..3,0..3,0..3);
for i from 0 to 3 do
for j from 0 to 3 do
for k from 0 to 3 do
A[i,j,k]:=trace( evalm( SSS||i||j||k &* LU) )/8;
od;
od;
od;
RETURN(eval(A));
end:

> RPauli3:=proc(rho)
local R,RSUM,i,j,k;
#
# Expresses density operator rho as a formal sum of sigma.i(X)sigma.j(X)sigma.k's
#
R:=Pauli3(rho);
RSUM:=0;
for i from 0 to 3 do for j from 0 to 3 do for k from 0 to 3 do
RSUM:=RSUM + R[i,j,k]*S||i||j||k;
od; od; od;
RETURN(RSUM);
end:

> Pauli4:=proc(LU)
# This procedure expresses a 16x16 skew-Hermitian matrix as
# the sum of the tensor product of Pauli matrices.
# LU:= Sum( i=0..3, j=0..3, k=0..3, m=0..3, R[i,j,k,m]*SSS.i.j.k.m )
# Please note that the indices of the matrix R run from 0 to 3
local A, i, j, k, m;

A:=array(0..3,0..3,0..3,0..3);
for i from 0 to 3 do for j from 0 to 3 do
for k from 0 to 3 do for m from 0 to 3 do
A[i,j,k,m]:=trace( evalm( SSS||i||j||k||m &* LU) )/16;
od; od;
od; od;
RETURN(eval(A));
end:

> RPauli4:=proc(rho)
local R,RSUM,i,j,k,m;
#
# Expresses density operator rho as a formal sum of
# sigma.i(X)sigma.j(X)sigma.k'(X)sigma.m's
#
R:=Pauli4(rho);
RSUM:=0;
for i from 0 to 3 do for j from 0 to 3 do for k from 0 to 3 do for m from 0 to 3 do
RSUM:=RSUM + R[i,j,k,m]*S||i||j||k||m;
od; od; od; od;
RETURN(RSUM);
end:

> FLL4:=proc(LU)
# This procedure expresses a 4x4 skew-smmetric matrix as
# the sum a linear combination of LL.i's .
# LU:= Sum( i=1..6, A[1,i]*LL.i )
local A, i;

A:=matrix(1,6);
for i from 1 to 6 do
A[1,i]:=trace( evalm( LL||i &* LU) )/(-2);
od;
RETURN(evalm(A));
end:

> Matrep:=proc(P, P_Deg)
# This procedure converts a permutation represented as
# a product of (not necessarily disjoint) cycles to
# a permutation matrix
local f,Pmat,n,j,Cycle_j,n_j,a_j,i,b_j;

f:=(ii,jj)->if(ii=jj) then 1 else 0 fi;
Pmat:=matrix(P_Deg, P_Deg, f);
n:=nops(P);
if (n<1) then
RETURN(evalm(Pmat));
else
for j from 1 to n do
Cycle_j:=op(j,P);
n_j:=nops(Cycle_j);
if (n_j>1) then
a_j:=op(1,Cycle_j);
for i from 2 to n_j do
b_j:=op(i,Cycle_j);
Pmat:=swaprow(Pmat,a_j+1,b_j+1);
od;
fi;
od;
fi;
RETURN(evalm(Pmat));
end:

> RemoveIPi:=proc(A)
local ff, M;

ff:=(i,j)-> if (i=j) then 1/(Pi*I) else 0 fi;
M:=evalm( matrix(rowdim(A),coldim(A),ff) &* A);
RETURN(evalm(M));
end:

> f||0:=x->2*x:
f||1:=x->irem(x,2 ) + ( iquo(x,2 )*2^2):
f||2:=x->irem(x,2^2) + ( iquo(x,2^2)*2^3):
f||3:=x->irem(x,2^3) + ( iquo(x,2^3)*2^4):
f||4:=x->irem(x,2^4) + ( iquo(x,2^4)*2^5):
f||5:=x->irem(x,2^5) + ( iquo(x,2^5)*2^6):

> Trace1:=proc(A)
local ans, ii, jj;
#
ans:=matrix(rowdim(A)/2, coldim(A)/2 );
for ii from 1 to rowdim(A)/2 do;
for jj from 1 to coldim(A)/2 do;
ans[ii,jj]:=A[1+f||0(ii-1),1+f||0(jj-1)] +A[1+f||0(ii-1)+1,1+f||0(jj-1)+1];
od;
od;
RETURN(evalm(ans))
end:
#
Trace2:=proc(A)
local ans, ii, jj;
#
ans:=matrix(rowdim(A)/2, coldim(A)/2 );
for ii from 1 to rowdim(A)/2 do;
for jj from 1 to coldim(A)/2 do;
ans[ii,jj]:=A[1+f||1(ii-1),1+f||1(jj-1)] +A[1+f||1(ii-1)+2,1+f||1(jj-1)+2];
od;
od;
RETURN(evalm(ans))
end:
#
Trace3:=proc(A)
local ans, ii, jj;
#
ans:=matrix(rowdim(A)/2, coldim(A)/2 );
for ii from 1 to rowdim(A)/2 do;
for jj from 1 to coldim(A)/2 do;
ans[ii,jj]:=A[1+f||2(ii-1),1+f||2(jj-1)] +A[1+f||2(ii-1)+4,1+f||2(jj-1)+4];
od;
od;
RETURN(evalm(ans))
end:
#
Trace4:=proc(A)
local ans, ii, jj;
#
ans:=matrix(rowdim(A)/2, coldim(A)/2 );
for ii from 1 to rowdim(A)/2 do;
for jj from 1 to coldim(A)/2 do;
ans[ii,jj]:=A[1+f||3(ii-1),1+f||3(jj-1)] +A[1+f||3(ii-1)+8,1+f||3(jj-1)+8];
od;
od;
RETURN(evalm(ans))
end:
#
Trace5:=proc(A)
local ans, ii, jj;
#
ans:=matrix(rowdim(A)/2, coldim(A)/2 );
for ii from 1 to rowdim(A)/2 do;
for jj from 1 to coldim(A)/2 do;
ans[ii,jj]:=A[1+f||4(ii-1),1+f||4(jj-1)] +A[1+f||4(ii-1)+16,1+f||4(jj-1)+16];
od;
od;
RETURN(evalm(ans))
end:
#
Trace6:=proc(A)
local ans, ii, jj;
#
ans:=matrix(rowdim(A)/2, coldim(A)/2 );
for ii from 1 to rowdim(A)/2 do;
for jj from 1 to coldim(A)/2 do;
ans[ii,jj]:=A[1+f||5(ii-1),1+f||5(jj-1)] +A[1+f||5(ii-1)+32,1+f||5(jj-1)+32];
od;
od;
RETURN(evalm(ans))
end:
#
PTrace:=(n,rho)->(Trace||(n+1))(rho):

> Fourier:=proc(n)
local omega,A,i,j;
# This program computes the nxn matrix ( omega^ij )/sqrt(n), where
# omega is the primitive n-th root of unity given by omega=exp(2PiI/n)
omega:=exp(2*Pi*I/n);
A:=matrix(n,n);
for i from 0 to n-1 do for j from 0 to n-1 do
A[i+1,j+1]:=evalc(omega^(i*j));
od; od;
A:=evalm( 1/sqrt(n) * A);
RETURN(evalm(A));
end:

> Vec:=proc(A)
local m,n,V_Dim,V,i,j;
# Vec transforms an mxn matrix A into an m*n column vector V, column by column
m:=rowdim(A); n:=coldim(A); V_Dim:=m*n;
V:=vector(V_Dim);
for j from 0 to n-1 do for i from 0 to m-1 do
V[i+m*j +1]:=A[i+1,j+1];
od; od;
RETURN(evalm(V));
end:
#
#Vec:=proc(M)
# Vec(M) = Matrix( Dim,1, [col1(M),col2(M), ... , col.CD(M)]),
# where Dim = rowdim(M)*coldim(M), and CD = coldim(M)
#
# local RD,CD,Dim,V,i,j;
#
# RD:=rowdim(M); CD:=coldim(M); Dim:=RD*CD;
# V:=matrix(Dim,1);
# for j from 1 to CD do for i from 1 to RD do
# V[i+(j-1)*RD,1]:=M[i,j];
# od; od;
# RETURN(evalm(V));
#end:

> Ad:=proc(Q)
local M,p,q,Apq,R_pq,Vpq,i;
# Computes Ad_Q as a 16x16 matrix; Q is assumed to be a 4x4 unitary
#
M:=matrix(16,16);
for q from 0 to 3 do print([`q=`||q]); for p from 0 to 3 do
#print([p,q]);
Apq:=evalm( Q &* SS||(p+4*q) &* My_Adjoint(Q) );
R_pq:=Pauli2(Apq);
Vpq:=Vec(R_pq);
for i from 0 to 15 do
#print([i]);
M[i+1,p+4*q +1]:=Vpq[i+1];
#print([i]);
od;
od; od;
RETURN(evalm(M));
end:

> Small_ad:=proc(Q)
local M,p,q,Apq,R_pq,Vpq,i;
# Computes Small ad_Q as a 16x16 matrix; Q is assumed to be a 4x4 skew Hermitian
# Small_ad(Q) (M) = ad_iQ(M), where M is assumed to be skew Hermitian
#
M:=matrix(16,16);
for q from 0 to 3 do print([`q=`||q]); for p from 0 to 3 do
#print([p,q]);
Apq:=evalm( evalm( evalm( Q &* SS||(p+4*q) ) - evalm( SS||(p+4*q) &* Q ) ) );
Apq:=evalm( I * Apq);
R_pq:=evalm(-I*Pauli2(Apq));
Vpq:=Vec(R_pq);
for i from 0 to 15 do
#print([i]);
M[i+1,p+4*q +1]:=Vpq[i+1];
od;
od; od;
RETURN(evalm(M));
end:

> FCom:=proc(A,B)
local X;
# FCom(A,B) is the commutator [A,B] of matrices A and B
#
X:=evalm( evalm( A &* B ) - evalm( B &* A ) );
RETURN(evalm(X));
end:

> TraceRL:=proc(A)
local B;
# B[i1,j0] = Sum(a=0..1, A[a+2*i1,j0+2*a] )

B:=matrix(2,2);
B[1,1]:=A[1,1]+A[2,3]; B[1,2]:=A[1,2]+A[2,4];
B[2,1]:=A[3,1]+A[4,3]; B[2,2]:=A[3,2]+A[4,4];
RETURN(evalm(B));
end:
#
TraceLR:=proc(A)
local B;
# B[i0,j1] = Sum(a=0..1, A[i0+2*a,a+2*j1] )

B:=matrix(2,2);
B[1,1]:=A[1,1]+A[3,2]; B[1,2]:=A[1,3]+A[3,4];
B[2,1]:=A[2,1]+A[4,2]; B[2,2]:=A[2,3]+A[4,4];
RETURN(evalm(B));
end:

> My_Conj:=x->eval(subs( I=-I,x)):
My_Adjoint:=A->map( My_Conj, transpose(A)):
Ket:=(a,b)->matrix(2,1,[cos(a),expand((exp(I*b)*sin(a)))]):

> Id:=proc(n)
# Id(n) = nxn identity matrix
#
local f,i,j, A;
f:=(i,j)->if (i=j) then 1 else 0 fi:
A:=matrix(n,n,f);
RETURN(evalm(A));
end:

> CNOT12:=matrix(4,4,[1,0,0,0, 0,1,0,0, 0,0,0,1, 0,0,1,0]):
CNOT21:=matrix(4,4,[1,0,0,0, 0,0,0,1, 0,0,1,0, 0,1,0,0]):

> Hadamard:=proc(n)
# Hadamard(n) = (X)_(i=1..n) H, where H is the 2x2 Hadarmard matrix
local IR2,H,Ans,i;

IR2:=1/sqrt(2);
H:=matrix(2,2,[IR2,IR2,IR2,-IR2]);
Ans:=evalm(H);
if (n<>1) then
for i from 1 to n-1 do
Ans:=evalm( kr(Ans,H));
od;
fi;
RETURN(evalm(Ans));
end:

> Int_To_Bin:=proc(n,lgth)
# Int_To_Bin(n) = postive integer n written as a binary number of length lgth
local nn,s,i,nx;

nn:=n;
s:=``;
for i from 1 to lgth do
nx:=irem(nn,2);
s:=cat(nx,s);
nn:=iquo(nn,2);
od;
RETURN(s);
end:

> RKet:=proc(ket)
# RKet(ket) = column vector in symbolic ket form
local r,S,i,logr,bbb;
#
r:=rowdim(ket);
logr:=ceil( evalf( log[2](r) ) );
S:=0;
for i from 1 to r do
bbb:=Int_To_Bin(i-1,logr);
S:=S + ket[i,1]*(`|`||bbb||`>`);
od;
RETURN(S);
end:

> CNOT:=proc(T,C,W)
# CNOT(T,C,W) = permutation (written as a product of disjoint cyles) which represents
# a CNOT wit target bit T, control bit C, and with W wires labeled 0..(W-1).
# Convention: First=Bottom =Right and Last=Top=Left
local ET,EC,SETEC,EW,L,k,BITC,BITT;

ET:=2^T; EC:=2^C; SETEC:=ET+EC; EW:=2^W;
L:=[];
for k from 0 to (EW-1) do
BITC:=irem(floor(k/EC),2); BITT:=irem(floor(k/ET),2);
if (BITC+BITT=0) then
L:=[op(L), [EC+k, SETEC+k] ];
fi;
od;
RETURN(L);
end:

> Entropy:=proc(A)
local EV,nEV,Ans,n;
# Entropy of a Hermitian operator
# This procedure has not been completed
# The problem is that the eigenvalues, although
# real, are not always non-negative. A second
# problem is how to normalize the Hermitian matrix
# before computing the entropy. Should one make it
# of 2-norm one? Should one make it of trace one?
#
EV:=vector([eigenvals(evalm((-I)*A))]);
nEV:=1/norm(EV,2);
EV:=evalm( nEV * EV);
print(`EV=`,EV);
Ans:=0;
for n from 1 to nops(EV) do
Ans:=Ans - EV[n]*Logz(EV[n]);
od;
RETURN(Ans):
end:

> # norm2(X) = the square root of the spectral radius
# of evalm( X &* map(conjugate, transpose(X)))
norm2:=X-> sqrt(max(op(map(abs,[eigenvals(
evalm( X &* map(conjugate,transpose(X)) )
)])))):

> # List of functions included above:
# QF:=proc(U) ... Computes a unitary matrix QF such that QF &* U &* is a diagonal matrix
# Logz(Z) computes the natural log of z
# LogU:=proc(U) ... Computes natural log of the square matrix U
# kr:=proc(A,B) ... Computes the tensor product A(X)B of matrices A and B
# Pauli1:=proc(U)... Expresses a 2x2 Hermitian matrix as a linear sum of sigma.i's
# RPauli1:=proc(LU)... Expresses a 2x2 Hermitian matrix as a linear formal sum of sigma.i's
# Pauli2:=proc(LU)... Expresses a 4x4 Hermitian matrix as a linear sum of sigma.i(X)sigma.j's
# RPauli2:=proc(rho)... Expresses density operator rho as a formal sum of sigma.i(X)sigma.j's
# Inv_Pauli2:=proc(R)...For a given a 4x4 matrix R over the reals, this procedure
# computes rho = Sum(p=0..3, q=0..3) R[p+1,q+q]*sigma.p(X)sigma.q
# Pauli3:=proc(LU)... Expresses a 8x8 Hermitian matrix as a linear sum of
# sigma.i (X) sigma.j (X) sigma.k's
# RPauli3(rho) ... Expresses density operator rho as a formal sum of sigma.i(X)sigma.j(X)sigma.k's
# Pauli4:=proc(LU)... Expresses a 16x16 Hermitian matrix as a linear sum of
# sigma.i (X) sigma.j (X) sigma.k (X) sigma.m's
# RPauli4(rho) ... Expresses density operator rho as a formal sum of
# sigma.i(X)sigma.j(X)sigma.k(X)sigma.m's
# FLL4:=proc(LU) ... Expresses a skew symmetric 4x4 matrix as a linear combination of LL.i's
# Matrep:=proc(P, P_Deg) ... Converts the degree P_Deg permutation P (represented as
# a list of lists) into a P_Deg x P_Deg permutation matrix
# RemoveIPi:=proc(A) = A/(I*Pi), where A is a square matrix
# f.0, f.1, f.2, f.3, f.4, f.5 are functions called by Trace1, Trace2, Trace3
# Trace1, Trace2, Trace3, Trace4, Trace5, Trace6 are partial trace functions
# PTrace(n,rho) = partial trace of the n-th qubit = (Trace.(n+1))(rho)
# Qubits are listed right to left, and labeled 0,1, ...
# n = 0,1,2,3,4,5
# sigma.0, sigma.1, sigma.2, sigma.3 are predefined as the Pauli matrices
# SSS.i=sigma.i
# SS.(i+4*j) = kr(sigma.i,sigma.j) are predefined
# SSS.i.j=SS.i.j = kr(sigma.i,sigma.j) are predefined
# SSS.i.j.k = kr( SS.(i+4*j), sigma.k ) are predefined
# SSS.i.j.k.m = kr( SS.(i+4*j), SS(k+4*m) ) are predefined
# L3x3.j ... j=1..3 ... are predefined infinitesimal generators of SO(3)
# LL.j ... j=1..6 ... are predefined infinitesimal generators of SO(4)
# ket.0 = |0> = (1.0)^Transpose ket.1 =|1> = (0,1)^Transpose
# ket.i.j = ket.i (X) ket.j
# ket.i.j.k = ket.i (X) ket.j (X) ket.k
# ket.i.j.k.m = ket.i (X) ket.j (X) ket.k (X) ket.m
# Ket_To_Rho(ket) = density operator corresponding to ket
# Bell.i.j ..... Bell basis kets, i = 0..1, j = 0..1
# Bell.i.j.k ... Bell basis kets, i = 0..1, j = 0..1, k = 0..1
# Bell.i.j.k.m ... Bell basis kets, i = 0..1, j = 0..1, k = 0..1, m = 0..1
# Bell4(a,b,c,d) = ( ket.(1-a).(1-b).(1-c).(1-d) + ((-1)^a)*ket.a.b.c.d)/sqrt(2)
# where a,b,c,d are elts of {0,1}
# Sam.i.j.k .... Sam basis kets, i = 0..1, j = 0..1, k = 0..1
# Fourier:=proc(n) computes the nxn matrix ( omega^ij )/sqrt(n), where omega=exp(2PiI/n)
# Vec:=proc(A) ... transforms an mxn matrix A into an m*n column vector V, column by column
# Ad:=proc(Q) ... Computes Ad_Q as a 16x16 matrix; Q is assumed to be a 4x4 unitary matrix
# FCom(A,B) is the commutator [A,B] of matrices A and B
# (Small_ad:=proc(Q))(M) = Small_Ad_iQ(M), where Q is Hermitian and M is skew-Hermitian
# TraceRL(A) = Sum(a=0..1, A[a+2*i1,j0+2*a] )
# TraceLR(A) = Sum(a=0..1, A[i0+2*a,a+2*j1] )
# My_Conj(z) = the conjugate of the comple number z
# My_Adjoint(A) = the adjoint of the complex square matrix A
# Ket:=(a,b)->matrix(2,1,[cos(a),(exp(I*b)*sin(a))]):
# Vec:=proc(M), Vec(M) = Matrix( Dim,1, [col1(M),col2(M), ... , col.CD(M)]),
# where Dim = rowdim(M)*coldim(M), and CD = coldim(M)
# Id(n) = nxn identity matrix
# CNOT12:= (i,j)->(i,i+j) .... CNOT(21):=(i,j)->(i+j,j)
# Hadamard(n) = (X)_(i=1..n) H, where H is the 2x2 Hadarmard matrix
# Int_To_Bin(n) = postive integer n written as a binary number of length lgth
# RKet(ket) = column vector in symbolic ket form
# CNOT(T,C,W) = permutation (written as a product of disjoint cyles) which represents
# a CNOT wit target bit T, control bit C, and with W wires labeled 0..(W-1).
# Convention: First=Bottom =Right and Last=Top=Left
# norm2(X) = the square root of the spectral radius
# of evalm( X &* map(conjugate, transpose(X)))

>

>