function [C,degC] = polmul(A,degA,B,degB,EPS)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% [C,degC]=polmul(A,degA,B,degB,EPS)                          %
%                                                             %
% MATLAB FUNCTION polmul FOR THE MULTIPLICATION               %
% OF TWO POLYNOMIAL MATRICES                                  %
%                                                             %
% Version with  clearing by  EPS                              %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% S. Pejchova, May  26st, 1993                                %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% setting of EPS                                              %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if nargin==4
   EPS=1e-8;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% size of C                                                   %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[rA,cA] = polsize(A,degA);
[rB,cB] = polsize(B,degB);
if cA ~= rB
   error('polmulp: Inconsistent dimensions of input matrices');
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% degC                                                        %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

degC = degA+degB;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% degC < 0                                                    %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if degC < 0
  C=zeros(rA,cB);
  return
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% preliminaries                                             %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

AZ=zeros(rA,cA);
BZ=zeros(rB,cB);
AL=[]; BL=[]; C=[]; N=[];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% finding the min. entries                                   %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

mina=1; minb=1; mr=1;

for i=1:rA
    MA=find(A(i,:));
    if isempty(MA)==0
       mr=min(abs(A(i,MA)));
    end
    if mr<mina
       mina=mr;
    end
end

mr=1;
for i=1:rB
    MB=find(B(i,:));
    if isempty(MB)==0
       mr=min(abs(B(i,MB)));
    end
    if mr<minb
       minb=mr;
    end
end

EP=mina*minb*EPS;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% multiplication                                             %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for i=1:degC+1

    if i<= degA+1
       AL=[AL,A(:,(i-1)*cA+1:i*cA)];
    else
       AL=[AL,AZ];
    end

    if i<= degB+1
       BL=[B(:,(i-1)*cB+1:i*cB); BL];
    else
       BL=[BZ;BL];
    end

    CL=AL*BL;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% clearing CL                                               %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    for j=1:rA
        MC=find(abs(CL(j,:))<=EP);

        if isempty(MC)==0
           [k,l]=size(MC);
           CL(j,MC)=zeros(k,l);
        end
    end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% creating  N                                               %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    if CL==zeros(rA,cB)
       N=[N,0];
    else
       N=[N,1];
    end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% creating C                                                %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    C=[C,CL];

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% setting degC                                               %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if isempty(find(N))==1
   C=zeros(rA,cB);
   degC=-inf;
else

degC=max(find(N))-1;
C=C(:,1:cB*(degC+1));
end