% This function takes as argument a string strA which describes a m\times m
% matrix A and an integer k which depends of a paramater a. The 
% coefficients of strA must be separated by
% a comma.
% The routine returns the following values:
% matrixA         : A string which describes a function of a which returns 
%                   A(a) when given a as argument.
% compound_matrix : A string which describes a function of a which returns 
%                   the n^{th}-compound matrix of A(a) when given a as 
%                   argument.
%
% When m=2k, the routine also returns the following matrices:
%
% sigma_matrix    : A matrix such that U\wedge V=U^T sigma_matrix V for U,V 
%                   two n-forms written in the standard basis.
% maslov_form     : A line matrix such that maslov_form*U gives the Maslov 
%                   angle associated to U.
% maslov_matrix   : A matrix such that maslov_matrix*U is the polynomial 
%                   defined by det((V-iW)-lambda(V+iW)), where V and W are
%                             ( V )
%                   such that (   ) has its minor determinants equal to U.
%                             ( W )
function[matrixA,compound_matrix,sigma_matrix,maslov_form,maslov_matrix]...
=exterior_algebra(strA,k)

    function[y,N]=Afun(i,j)
        sep=0;
        for rr=1:length(strA)
            if (strA(rr)==',' | strA(rr)==';' | strA(rr)==' ')
                sep=[sep,rr];
            end
        end
        sep=[sep,length(strA)+1];
        N=sqrt(length(sep)-1);
        alpha=(j-1)+N*(i-1)+1;
        y=char(strA(sep(alpha)+1:sep(alpha+1)-1)');
    end

    function[val]=isazero(str)
        if (length(str)==0 | (length(str)==1 & str(1)=='0'))
            val=(0==0);
        else
            val=(1==0);
        end
    end

    [a11,n]=Afun(1,1);

    subsets=nchoosek((1:n),k);
    N=size(subsets,1);


    %Computation of the compound matrix
    compound_matrix=sprintf('%s\n','function[y]=compound_matrix(a)');
    compound_matrix=sprintf('%sy=sparse(%d,%d);\n',compound_matrix,N,N);
    for a=1:N
        for b=1:N
            I=subsets(a,:);
            J=subsets(b,:);
            U=zeros(n,1);                                                                                                             
            U(I)=U(I)+1;
            U(J)=U(J)+1;
            R=find(U==1);
            if length(R)==2 
                if length(find(I==R(1)))==0
                R([2 1])=R; 
                end
                r=find(I==R(1));
                s=find(J==R(2));
                if ~(isazero(Afun(I(r),J(s))))
                     if mod(r+s,2)==0
                         compound_matrix=sprintf('%sy(%d,%d)=%s;\n', ...
                            compound_matrix,a,b,Afun(R(1),R(2)));
                     else
                         compound_matrix=sprintf('%sy(%d,%d)=-(%s);\n', ...
                            compound_matrix,a,b,Afun(R(1),R(2)));
                     end
                end
            end
            if a==b
                for r=1:k
                    if ~(isazero(Afun(I(r),J(r))))
                        compound_matrix=sprintf('%sy(%d,%d)=(%s)', ...
                            compound_matrix,a,b,Afun(I(r),J(r)));
                        for r=r+1:k
                            if ~(isazero(Afun(I(r),J(r))))
                                compound_matrix=sprintf('%s+(%s)', ...
                                    compound_matrix,Afun(I(r),J(r)));
                            end
                        end
                        compound_matrix=sprintf('%s;\n',compound_matrix);
                    end
                end
            end
        end
    end




    sigma_matrix=sparse(N,N);
    for a=1:N
        sigma_matrix(a,N+1-a)=(-1)^(sum(subsets(a,:)-(1:k)));
    end

    if (mod(n,2)==0 & k==n/2)
        maslov_form=sparse(1,N);
        maslov_matrix=sparse(k+1,N)
        for a=1:N
            I=subsets(a,find(subsets(a,:)<=k));
            J=subsets(a,find(subsets(a,:)>k));
            if prod(double(sort([I,(J-k)])==sort((1:k))))==1
                maslov_form(1,a)=...
                    i^(k-length(I))*(-1)^(sum(I-(1:length(I))));
                
                maslov_matrix(:,a)=...
                    (poly([-ones(k-length(I),1);ones(length(I),1)])')*...
                        (-1)^(length(I))*...
                        i^(k-length(I))*(-1)^(sum(I-(1:length(I))));                
            end
        end
    end

    matrixA=sprintf('%s\n y=[','function[y]=matrixA(a)');
    for r=1:n
        for s=1:n
            if isazero(Afun(r,s))
                matrixA=sprintf('%s %s',matrixA,'0');
            else
                matrixA=sprintf('%s (%s)',matrixA,Afun(r,s));
            end
        end
        if r==n
            matrixA=sprintf('%s];',matrixA);
        else
            matrixA=sprintf('%s; ... \n', matrixA);
        end
    end

end
