decompose.m

该算法是用于量子运算的matlab程序

% 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 "zyx" or "Z^{(1)}Y^{(2)}X^{(3)}".
%    Form decompose(rho,mode,threshold) makes it possible
%    to give a threshold. Below this value a correlation is
%    considered zero. Default is 1e-14.

% Copyright (C) 2005  Geza Toth    E.mail: [e-mail]
%
% This program is free software; you can redistribute it and/or
% modify it under the terms of the GNU General Public License
% as published by the Free Software Foundation; see gpl.txt
% of this subroutine package.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
% 
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc., 51 Franklin Street, Fifth Floor, 
% Boston, MA  02110-1301, USA.

function pstring=decompose(rho,varargin)

% Corraltions 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