function varargout = maxEfrac(rho)
% maxEfrac   Maximum entangled fraction
% requires: unpauli.m
% author: Toby Cubitt
% license: GPL2
%
%    S = maxEfrac(RHO) returns the maximum entangled fraction (also
%    known as maximum singlet fraction) of 2x2 density matrix RHO,
%    defined as the
%
%      max phi'*rho*phi
%
%    where phi is a two-qubit maximally entangled state
%
%    [O1,S,O2] = maxEfrac(RHO) in addition returns the orthogonal
%    matrices that achieve the maximum in the Hilbert-Schmidt
%    representation.


%% Copyright (C) 2004-2009 Toby Cubitt
%%
%% 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; either version 2
%% of the License, or (at your option) any later version.
%%
%% 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.


R = pauli(rho);
T = abs(R(2:4,2:4)).*sign(real(R(2:4,2:4)));

[O1,S,O2] = svd(T);

S = diag([1,1,det(O1)])*S*diag([1,1,det(O2)]);
O1 = O1*diag([1,1,det(O1)]);
O2 = O2*diag([1,1,det(O2)]);

Q = (O2*diag([1,1,-1])*O1.').';
P(1,1) = 1;
P(2:4,2:4) = Q;

phi = unpauli(P);

varargout{1} = trace(phi*rho);
if nargout >= 2
  varargout{2} = phi;
end
if nargout == 4
  varargout{3} = O1;
  varargout{4} = O2;
end