function [d,degd] = polred(b,degb,a,dega,c,k)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% performs forced polynomial reduction              % 
%                 d(s) = b(s) + a(s)*c*s^k          % 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% M. Sebek, July 1990                               %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


if k<0
  error('polred: negative degree');
end

n=max([degb dega+k]);
if n<0
 d=0;
 degd=-Inf;
 return
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% setting EPS                                       %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

EPS=1e-8;     

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% reduction                                         %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if dega<0
  d=b;
elseif degb<0
  d=c*[zeros(1,k) a];
else
  b=[b zeros(1,n-degb)];
  a=[zeros(1,k) a zeros(1,n-dega-k)];
  d=b+c*a;

  for i=n+1:-1:1
    if abs(d(i))<=max(abs([c*a(i) b(i)]))*EPS
     d(i)=0;
    end
  end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% computing degd                                    %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

degd=n;
while d(degd+1)==0 
  d(degd+1)=[];
  degd=degd-1;
  if degd<0
    d=0;
    degd=-Inf;
    return
  end
end