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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% MATLAB FUNCTION poladd FOR THE ADDITION                   %
% OF TWO POLYNOMIAL MATRICES                                %
%                                                           %
% Version with clearing by EPS                              %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% H. Kwakernaak, August, 1990                               %
% modified by M. Sebek, October 1990                        %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


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

EPS= 1e-8;

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

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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% one matrix has deg -Inf                                   %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if degA < 0
  C=B;
  degC=degB;
  return;
end

if degB < 0
  C=A;
  degC=degA;
  return;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% addition                                                  %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


degC = max(degA,degB);
if degC >= degA
   A = [A zeros(rA,(degC-degA)*cA)];
end
if degC >= degB
   B = [B zeros(rB,(degC-degB)*cB)];
end
C = A+B;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% clearing                                                  %                
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

k=degC*cA;

for i=1:rA
  for j=1:cA
    l=k+j;
      while (max(abs(A(i,l)),abs(B(i,l)))*EPS >= abs(C(i,l)))
        C(i,l)=0;
        l=l-cA;
	if l<1, break,end
      end
  end
end

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

if degC>0
   while C(:,cB*degC+1:cB*(degC+1))==zeros(rA,cB)
      C(:,cB*degC+1:cB*(degC+1))=[];
      degC=degC-1;
      if degC<1, return, end
   end
end