```% decompose   Creates a string with the Pauli decomposition of a
%             Hermitian matrix.
%    Form of use: decompose(rho) where rho is the density matrix.
%    One can print in LaTeX format with decompose(rho,1).
%    For the numbering of qubits we note that mkron(x,y,z) is
%    decomposed into "xyz" or "X^{(1)}Y^{(2)}Z^{(3)}".
%    For the LaTeX output the numbering is different from
%    the one used in other commands (e.g., reorder) in which
%    the numbering of qubits were 3 2 1 and not 1 2 3. This is
%    to help to produce a LaTeX formula that can directly be
%    used in a paper. Form decompose(rho,mode,threshold)
%    makes it possible to give a threshold. Below this
%    value a correlation is considered zero. Default is 1e-14.

function pstring=decompose(rho,varargin)

% Correlations smaller than that are taken to be zero
mincorr=1e-14;

LaTeXON=0;
if ~isempty(varargin),
if length(varargin)>=1,
LaTeXON=varargin{1};
end %if
if length(varargin)==2,
mincorr=varargin{2};
end %if
if length(varargin)>2,
error('Wrong number of input arguments.')
end %if
end %if
x=[0 1;1 0];
z=[1 0;0 -1];
y=i*x*z;
e=eye(2);

% Use sparse matrices
e=sparse(e);
x=sparse(x);
y=sparse(y);
z=sparse(z);

% rho=ketbra2(rho);

[sy,sx]=size(rho);

N=log2(sx);

pstring='';

OPstr='1xyz';
OPstr2='EXYZ';

opindex=zeros(N);

while 1
switch opindex(1)
case 0, op=e;
case 1, op=x;
case 2, op=y;
case 3, op=z;
end %switch
for m=2:N
switch opindex(m)
case 0, op=kron(op,e);
case 1, op=kron(op,x);
case 2, op=kron(op,y);
case 3, op=kron(op,z);
end %switch
end %for
if trace(rho*op)~=0,
% Correlations: Use "real" to get rid off small imaginary parts
corr=real(trace(rho*op))/2^N;
termstr='';
% Construct string for the term
if LaTeXON,
if opindex==zeros(N),
termstr='E';
else
for m=1:N
if opindex (m) >0,
termstr=[termstr,OPstr2(opindex(m)+1) '^{(' num2str(m) ')}' ];
end %if
end %for

end %if
else
for m=1:N
termstr=[termstr,OPstr(opindex(m)+1)];
end %for
end %if
if isempty(pstring),
if abs(corr-1)<mincorr,
pstring=[termstr];
elseif abs(corr+1)<mincorr,
pstring=['-' termstr];
elseif abs(corr)>mincorr,
pstring=[num2str(corr) '*' termstr];
end %if
else
if corr>mincorr,
if abs(corr-1)<mincorr,
pstring=[pstring '+' termstr];
else
pstring=[pstring '+' num2str(corr) '*' termstr];
end %if
elseif corr<-mincorr,
if abs(corr+1)<mincorr,
pstring=[pstring '-' termstr];
else
pstring=[pstring '-' num2str(abs(corr)) '*' termstr];
end %if
end %if
end %if
end %if

% Increase counter
k=N;
opindex(k)=opindex(k)+1;
while opindex(k)==4
opindex(k)=0;
k=k-1;
if k==0,
if isempty(pstring),
pstring='0';
end %if
return
end %if
opindex(k)=opindex(k)+1;
end %while

end %while

```