> 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))) > >