bin2ket::usage = "bin2ket[ls] takes list ls of n 0's and 1's,
thought of as basis states of n qubits, and returns the corresponding ket in a
2^n dimensional Hilbert space. E.g. bin2ket[{0,1}] -> {0,1,0,0}. One can
multiply by coefficients, and add if the number of qubits is the
same. c*bin2ket[{0,1}] + d*bin2ket[{1,1}] -> {0,c,0,d}. Also see bket[]."
bin2ket[ls_]:= Module[{ket,ln=Length[ls],m},
ket=Table[0,{2^ln}];
m=1+Fold[2*#1+#2&,0,ls];
++ket[[m]];
ket](*END bin2ket*)

bket::usage = "bket[bin,n]. Returns stadard form of basis ket
corresponding to |bin>, where 'bin' is a string of n 0's or 1's. E.g.,
bket[01,2]={0,1,0,0}; bket[010,3]={0,0,1,0,0,0,0,0}. Also see bin2ket"
(*bket:comment. Due to difficulty in getting Mca to interpret
'000' as different from '0', this function adds 2*10^n to the n-bit 'bin'
and converts the result to a string, which is then converted back to a list
by Characters[] followed by ToExpression. Finally, the 2 is discarded using
Take[], and the result converted to the corresponding binary number using
Fold[]. Adding 1 to the result yields m, the position where the list
representing the starndard ket is changed from 0 to 1.*)
bket[bin_,n_]:=Module[{lst=Table[0,{2^n}],m},
m=1+Fold[2*#1+#2&,0,
Take[ToExpression[ Characters[ ToString[2*10^n+bin] ] ],-n]];
lst[[m]]=1;lst] (*END bket*)

blochket::usage = "blochket[{x,y,z}] takes the Cartesian
coordinates of a point on the Bloch sphere and returns the corresponding ket in
the form {cos(th/2),sin(th/2)e^i*phi}."
(*blochket: The If[] is intended to suppress an error message from
ArcTan[] if both x and y are 0.*)
blochket[ls_]:= Module[{theta,phi,x=ls[[1]],y=ls[[2]],z=ls[[3]]},
theta = ArcCos[z];
If[ ((0==x)&(0==y))||(0.==x)&&(0.==y), phi=0, phi = ArcTan[x,y],
phi = ArcTan[x,y]];
{Cos[theta/2], Sin[theta/2]*E^(I*phi)}]

cgate::usage = "cgate[W_] returns a controlled-W on A x B,
where A is the control qubit and W a unitary on B (any dimension), as a matrix"
cgate[w_]:= tenprod[{{1,0},{0,0}},IdentityMatrix[Length[w]]] +
tenprod[{{0,0},{0,1}},w]

diags::usage = "diags[M] takes a matrix M as a list of lists,
and extracts the diagonal elements as a single list."
(*diags:comment. This is the inverse to Mca DiagonalMatrix*)
diags[mat_] := Module[{j,v=Table[0,{l=Length[mat]}]},
v=Table[0,{l}]; For[j=1,j<=l,++j, v[[j]] = mat[[j,j]] ]; v]

expandout::usage = "expandout[op,ls,dl] takes an operator op as
a matrix defined on a list ls of Hilbert spaces in the tensor product of spaces
with dimensions given by dl, and returns it as a matrix on the full space.
E.g., expandout[cnot,{3,2},{4,2,2}] gives a controlled not with the last qubit
(3rd space) the control."
expandout[op_,ls_,dl_]:=
permmat[tenprod[op,IdentityMatrix[Fold[Times,1,dl]/Length[op]]],
Join[ls,invertlist[Length[dl],Sort[ls]]],dl];

expandout2::usage = "expandout2[op,ls,n] takes a matrix op
representing a gate or other operation, a list ls of the qubits which forms the
basis of the matrix, and the total number n of qubits in circuit, and forms the
2^n by 2^n matrix representing that operation";
expandout2[op_,ls_,n_]:=
permmat2[addidents[op,n-Log[2,Length[op]]],Join[ls,invertlist[n,Sort[ls]]]];

ketcofs::usage = "ketcofs[v_,b_,dl_] returns a list of kets
which are the expansion coefficients of ket v in the orthonormal basis b (list
of basis vectors) of the first factor in a tensor product BC.... Here dl is the
list of dimensions of the factors, e.g., {3,4}, in which case b is a 3x3
matrix."
(*Comment. The Map[Flatten...] is needed in order that the
final output is a list of kets and not a list of tensors, in the case in which
dl contains more than two elements*)
ketcofs[v_,b_,dl_]:= Map[Flatten,Conjugate[b] . ket2kten[v,dl]];

ket2bin::usage = "ket2bin[ket] assumes list of length 2^n
represents n-qubit k, and produces a list where each member of the ket list is
associated with a symbol of type, say |010>. E.g., ket = {al,0,bt,2} yields
{{al,|00>},{bt,|10>},{2,|11>}}."
(*ket2bin:comment. It could undoubtedly be made more readable, but
this crude form has the advantage that the ket list can be either numerical, or
symbols, or a combination. If one had just numbers, replacing 0==item with
0.==Abs[Chop[item]] would be advantageous*)
ket2bin[ket_]:=Module[
{it,item,jt,lng=Length[ket],nlist,nn,olist={},str},
nn=IntegerExponent[lng,2];
For[it=0,it<lng,++it,
nlist = IntegerDigits[it,2,nn];
str="|";
For[jt=1,jt<=nn,++jt, str = str<>ToString[nlist[[jt]]]; ];
str = str<>">";
item = ket[[it+1]];
If[ 0==item, Continue[]];
AppendTo[olist,{item,str}]
]; olist](*END ket2bin*)

matinq::usage = "matinq[amat,bmat]=sums amat[[j,k]]*bmat[[j,k]]
over j and k. Here amat must be a matrix, bmat could be a tensor of rank >2."
matinq[amat_,bmat_]:= Module[{ln=Length[amat]},
Sum[amat[[j]].bmat[[j]],{j,ln}] ](*END matinq*)

mat2oten2::usage = "mat2oten2[mt] assumes mt is a
2^n x 2^n matrix for some integer n, and converts it to an o-tensor."
mat2oten2[mt_]:=nten2oten[ mat2nten2[mt] ]

mat2paul::usage = "mat2paul[mat] is the tensor
c[[j1,j2,...jn]] of coefficients of the expansion of the 2^n x 2^n matrix mat
in the form Sum c[[j1,j2,...jn]] sigma^1_j1 ... sigma^n_jn"
mat2paul[mat_]:=oten2paul[mat2oten2[mat]]

paulten::usage="paulten[1,0,3] will generate the Pauli tensor
corresponding to sg_x ox I ox sg_z, and similarly for any number of arguments,
each of which must be an integer between 0 and 3."
paulten[args__] := Module[{ls = List[args],ln,lsf,lsp,ptn},
ln=Length[ls]; lsp=1+ls; lsf=Table[4,{ln}];
ptn=Array[0*#&,lsf]; ptn = ReplacePart[ptn,1,lsp]; ptn] (*END paulten*)

permket2::usage = "permket2[kt,pm] returns ket for a tensor
product of qubits in the permuted order corresponding to pm. E.g. for kt
defined on A ox B ox C, pm={2,3,1}, the new ket is defined on C ox A ox B."
permket2[kt_,pm_]:= Flatten[ transpose[ket2kten2[kt],pm] ];

permptrace::usage = "permptrace[n,q] returns a permutation of
the integers 1,2,3, ... 2n such that 2q-1 and 2q are moved to the beginning of
the list. Thus for n=3 and q = 2 it returns 3,4,1,2,5,6. Used to form a
partial trace"
permptrace[n_,q_]:=Array[If[#<(2*q-1),#+2,If[#>(2*q),#,If[OddQ[#],1,2]]]&,2*n];

quadn::usage = "quadn[ml] is the sum of the absolute squares
of the elements in ml, whether a scalar, vector, matrix or tensor."
(*Comment: Use Re[] to get rid of 0.I terms in output*)
quadn[m_]:= Module[{lng=Length[m]},
If[0==lng,Return[Conjugate[m]*m]];
Re[Conjugate[Flatten[m]].Flatten[m]] ];

quadr::usage = "quadr[ml] is the sum of the squares
of the elements in ml, assumed to be a REAL vector or matrix or tensor."
quadr[m_]:= Module[{lng=Length[m]},
If[0==lng,Return[m^2]]; Flatten[m].Flatten[m]]

tenprod::usage = "tenprod[mt1,mt2,...] returns the matrix of
the tensor product mt1 0x mt2 0x ... The matrices may be rectangular."
(*tenprod:comment. outer[] transforms the na matrices read in
to an o-tensor, which is transformed to an n-tensor by transpose[] using
the permutation pm. The product of numbers of columns of matrices
read in = dim = number of columns of output matrix.*)
tenprod[args__]:=Module[{dim,las=List[args],na,pm},
na=Length[las]; pm = permon[na];
dim=prodlist[ Map[Last[Dimensions[#]]&,las] ];
Partition[ Flatten[ transpose[outer[args],pm] ] , dim ]](*END tenprod*)