function [AI,degAI,DE,degDE] = polinv(A,degA,mes)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% [AI,degAI,M,degM]=polinv(A,degA,mes)                  %
%                                                       %
% MATLAB FUNCTION polinv COMPUTE THE INVERSE MATRIX  AI %
% (A*AI=E) OF A NONSINGULAR SQUARE POLYNOMIAL MATRIX  A.%
%                                                       %
% PROCEDURE CALCULATES INVERSE FOR BOTH CASES  -        %
% IF A IS UNIMODULAR MATRIX OR NOT.                     %
%                                                       %
% IF MATRIX "A" IS UNIMODULAR, THEN MATRIX "DE"=[],     %
% IF NOT -  MATRIX "DE" IS ROW POLYNOM. MATRIX AND      % 
% DE(k) IS THE COMMON  DENOMINATORS OF THE RESPECTIVE   %
% COLUMNS OF "AI"                                       %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% S. Pejchova, March 7th,  1994                         %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% preliminary steps and initialization                  %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[rA,cA] = polsize(A,degA);

if rA ~= cA
   error('polinv: Input matrix is not square');
end

if nargin==2
   mes=0;
end

M=A; degM=degA;
Ds=[]; 
DE=ones(1,rA);
degDE=0;

V=eye(cA); degV=0;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Lower triangular form of A                           %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for k=1:rA

    [m,n]=size(M);
    Au(1,1:n)=M(k,:);
    [An,degAn,DEGM,LCM]=polmatcl(Au,degM);

    [Mi,Pi]=min(abs(DEGM(1,k:cA)));
    Pi=Pi+k-1;

    if k~=Pi
       EL=eye(cA);
       EL(k,k)=0;   EL(Pi,Pi)=0;
       EL(k,Pi)=1;  EL(Pi,k)=1;

       [M,degM]=polmul(M,degM,EL,0);
       a=DEGM(1,Pi);
       DEGM(1,Pi)=DEGM(1,k);
       DEGM(1,k)=a;
       [V,degV]=polmul(V,degV,EL,0);
    end

%%%%%%%%%%%%
    if mes>0
       disp('zamena sloupcu');
       k, M, degM,DEGM, keyboard
    end
%%%%%%%%%%%%

    [ak,degak]=polen(M,degM,k,k);
    [rV,cV]=polsize(V,degV);

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %  unimodular input matrix               %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    Ds =[Ds degak];

    if degak==0
       if ak~=1
          EL=eye(cA);
          EL(k,k)=1/(ak);
          ak=1/(ak);
          [M,degM]=polmul(M,degM,EL,0);
          [V,degV]=polmul(V,degV,EL,0);
       end

%%%%%%%%%%%%
       if mes>0
          disp('Jednicka na hl. diagonale');
          k, M,keyboard;
       end
%%%%%%%%%%%%

       for i=k+1 : cA
           if DEGM(1,i)>=0
              EM=[];
              [an,degan]=polen(M,degM,k,i);
              EL=eye(cA);
              [EM,degEM]=polseten(EL,0,-an,degan,k,i);
              [M,degM]=polmul(M,degM,EM,degEM);
              [V,degV]=polmul(V,degV,EM,degEM);

%%%%%%%%%%%%
              if mes>0
                 disp('Nulovani prvku vpravo - ak=1 ');
                 k, i, M, keyboard;
              end
%%%%%%%%%%%%

           end
       end

    else
       if degak<0
          error('polinv: Matrix is singular ');
       end
       for i=k+1 : cA
           if DEGM(1,i)>=0
              EM=[];
              [an,degan]=polen(M,degM,k,i);
              [u,degu,v,degv]=poldiv(an,degan,ak,degak);

              if degv ==-inf
                 EL=eye(cA);
                 [EM,degEM]=polseten(EL,0,-u,degu,k,i);
              else
                 EL=eye(cA);
                 [EM,degEM]=polseten(EL,0,ak,degak,i,i);
                 [EM,degEM]=polseten(EM,degEM,-an,degan,k,i);
              end

              [M,degM]=polmul(M,degM,EM,degEM);
              [V,degV]=polmul(V,degV,EM,degEM);
           end

%%%%%%%%%%%%
           if mes>0
              disp('Nulovani prvku vpravo - ak polynom');
              k, i, M, keyboard;
           end
%%%%%%%%%%%%

       end
    end

end

%%%%%%%%%%%%
if mes>0
   disp('Dolni trojuhelnikova matice');
   M,V, keyboard;
end
%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  Diagonal form                                %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[rM,cM]=polsize(M,degM);

for k=rA:-1:1
    Au(1,1:cM*(degM+1))=M(k,1:cM*(degM+1));
    [An,degAn,DEGM,LCM]=polmatcl(Au,degM);

    for i=1:k-1
       if DEGM(1,i)>=0
          [an,degan]=polen(M,degM,k,i);
          if Ds(k)==0
             EM=[];
             EL=eye(cA);
             [EM,degEM]=polseten(EL,0,-an,degan,k,i);
          else
             [ak,degak]=polen(M,degM,k,k);
             [u,degu,v,degv]=poldiv(an,degan,ak,degak);

             if degv == -inf
                EM=[];
                EL=eye(cA);
                [EM,degEM]=polseten(EL,0,-u,degu,k,i);
             else
                EM=[];
                EL=eye(cA);
                [EM,degEM]=polseten(EL,0,ak,degak,i,i);
                [EM,degEM]=polseten(EM,degEM,-an,degan,k,i);
             end
          end
          [M,degM]=polmul(M,degM,EM,degEM);
          [V,degV]=polmul(V,degV,EM,degEM);
       end

%%%%%%%%%%%%
          if mes>0
             disp('Nulovani levych prvku');
             k,i,M, keyboard;
          end
%%%%%%%%%%%%

    end

    [ap,degap]=polen(M,degM,k,k);
    [DE,degDE]=polseten(DE,degDE,ap,degap,1,k);


end

if degDE==0
   DE=[];
   degDE=-inf;
end


AI=V; degAI=degV;