% proj_sym  Gives the projector to the symmetric subspace
%    proj_sym(n)  Gives the projector to the symmetric subspace
%    of n qubits. proj_sym(n,d) does the same for qudits with dimension d.
%    Only the n=2 and the d=2, d=3 cases are implemented.

function P=proj_sym(N,varargin)

if isempty(varargin),
    d=2;
elseif length(varargin)==1,
    d=varargin{1};
else    
    error('Wrong number of input arguments.');
end %if

if N==2,
    P=zeros(d^2,d^2);
    for n=1:d
        for m=1:n-1
            v=zeros(d^2,1);
            v(n+(m-1)*d)=1;
            v(m+(n-1)*d)=-1;
            P=P+v*v'/2;   
        end %for
    end %for
    P=eye(d^2,d^2)-P;
elseif d==2,
    P=zeros(d^N,d^N);
    for n=0:N
        P=P+ketbra(dstate(n,N));
    end %for
elseif d==3
    su3;
    M1=coll(m1,N);
    M2=coll(m2,N);
    M3=coll(m3,N);
    M4=coll(m4,N);
    M5=coll(m5,N);
    M6=coll(m6,N);
    M7=coll(m7,N);
    M8=coll(m8,N);
    H=M1^2+M2^2+M3^2+M4^2+...
      M5^2+M6^2+M7^2+M8^2;
    me=max(real(eig(H)));
    [v,d]=eig(H);
    P=v*(me-0.05<d)*v';
elseif d==4
    paulixyz
    M1=coll(kron(x,e),N);
    M2=coll(kron(y,e),N);
    M3=coll(kron(z,e),N);
    M4=coll(kron(e,x),N);
    M5=coll(kron(e,y),N);
    M6=coll(kron(e,z),N);
    M7=coll(kron(x,x),N);
    M8=coll(kron(x,y),N);
    M9=coll(kron(x,z),N);
   M10=coll(kron(y,x),N);
   M11=coll(kron(y,y),N);
   M12=coll(kron(y,z),N);
   M13=coll(kron(z,x),N);
   M14=coll(kron(z,y),N);
   M15=coll(kron(z,z),N);
    H=M1^2+M2^2+M3^2+M4^2+...
      M5^2+M6^2+M7^2+M8^2+...
      M9^2+M10^2+M11^2+M12^2+...
      M13^2+M14^2+M15^2;
    me=max(real(eig(H)));
    [v,d]=eig(H);
    P=v*(me-0.05<d)*v';
elseif d==9
    su3;
    ma={eye(3)*sqrt(2/3),m1,m2,m3,m4,m5,m6,m7,m8};
    H=zeros(d^N,d^N);
    for n=1:9
        for m=1:9
            H=H+coll(kron(ma{n},ma{m}),N)^2;
        end %for
    end %for
    me=max(real(eig(H)));
    [v,d]=eig(H);
    P=v*(me-0.05<d)*v';
else 
    disp('Only d=2, d=3, d=4, d=9 or N=2 is realized.')
end %if