(* :Context: Pauli` *) (* :Title: Pauli *) (* :Author: Heiko Feldmann *) (* :Version: Mathematica 3.0 *) (* :Package Version: 0.5 *) (* :Summary: This package implements the composition and decomposition of 2^n*2^n matrices in terms of direct products of Pauli("\[Sigma]") matrices *) BeginPackage["Pauli`"]; \[Sigma]::usage = "\[Sigma][0] returns the 2x2 unit matrix, \[Sigma][i] returns the ith Pauli \ matrix."; CircleTimes::usage = "a \[CircleTimes] b \[CircleTimes] c or CircleTimes[a,b,c] returns the \ outer product of the matrices a,b,c." ; \[Sigma]Extract::usage = "\[Sigma]Extract[m,{ind}] gives the coefficient of the Pauli matrix indexed \ by ind in the \[Sigma]Decomposition of m."; \[Sigma]Decomposition::usage = "\[Sigma]Decomposition[m], where m is a matrix of size 2^n*2^n, gives the \ coefficients of a Pauli decomposition of m as a tensor of rank n"; \[Sigma]Composition::usage = "\[Sigma]Composition[c], where c is a tensor of rank n, gives the Pauli composed matrix of size 2^n*2^n with coefficients c"; Begin["`Private`"]; \[Sigma][0]={{1,0},{0,1}}; \[Sigma][1]={{0,1},{1,0}}; \[Sigma][2]={{0,-I},{I,0}}; \[Sigma][3]={{1,0},{0,-1}}; matrixFlatten[matrix_]:= Flatten[(Flatten[#,1]& /@ Transpose[#])& /@matrix,1] ; CircleTimes[a__]:= Nest[matrixFlatten,Outer[Times,a],Length[List[a]]-1]; Tr[m_]:= Sum[m[[i,i]],{i,1,Length[m]}]; TrProd[m1_,m2_]:= Sum[m1[[i,j]]*m2[[j,i]],{i,1,Length[m1]},{j,1,Length[m1]}]; \[Sigma]Extract[m_,indices_]:= TrProd[m,CircleTimes[Sequence @@ \[Sigma] /@ indices]]/Length[m]; \[Sigma]Decomposition[m_/;IntegerQ[Log[2,Length[m]]]] := Module[{tiefe = Log[2,Length[m]],indextable,indexlist}, indextable = Table[{{0},{1},{2},{3}},{Log[2,Length[m]]}]; indexlist = Outer[Join,Sequence @@indextable,1]; Map[\[Sigma]Extract[m,#]&, indexlist,{tiefe}]]; \[Sigma]Composition[c_]:= Module[{tiefe = TensorRank[c],indextable,indexlist,result =0}, indextable = Table[{{0},{1},{2},{3}},{tiefe}]; indexlist = Outer[Join,Sequence @@indextable,1]; indexlist = Flatten[indexlist,TensorRank[indexlist]-2]; For[i=1,i<=Length[indexlist],i++, result +=c[[Sequence @@ (indexlist[[i]]+1)]] * CircleTimes[Sequence @@ \[Sigma] /@ indexlist[[i]]]];result]; (* Error messages not jet implemented *) End[]; (* Pauli`Private` *) EndPackage[]; (* Pauli`*)