function [numPQ,denPQ]=combind(numPs,denPs,numQs,denQs)
% [numPQ,denPQ]=combind(numPs,denPs,numQs,denQs)
% Multiply P and Q and cancel the common roots (if posible).
% The function take the numerators and denominators of P and Q
% , and return the numerator and denominator of P*Q. 

[numPQ,temp1]=deconv(numPs,denQs);
[denPQ,temp2]=deconv(denPs,numQs);

% Make sure terms within 0.05 of each other cancel out

if sum(abs(temp1)) > 0.05
  if sum(abs(temp2)) > 0.05
   numPQ=conv(numPs,numQs);
   denPQ=conv(denPs,denQs);
  else
   numPQ=numPs;
   denPQ=conv(denPQ,denQs);
  end
 else
  if sum(abs(temp2)) > 0.05
   numPQ=conv(numPQ,numQs);
   denPQ=denPs;
  end
end

% Cancellation algorithm
rnum=[];
rden=[];
kn=length(roots(numPQ));
kd=length(roots(denPQ));

rnum=[roots(numPQ) zeros(1,kn)'];
rden=[roots(denPQ) zeros(1,kd)'];

k=0;
rcom=[];
for i=1:kn
 for j=1:kd
  if ( ( abs(rnum(i,1)-rden(j,1)) < 0.05*abs(rnum(i,1)) ) | ...
     ( abs(rnum(i,1)-rden(j,1)) < 0.05 ) ) & ...
     ( rnum(i,2) == 0 & rden(j,2) == 0)
     k=k+1;
     rcom(k)=(rnum(i,1)+rden(j,1))/2;
     rnum(i,2)=1;
     rden(j,2)=1;
   end
  end
 end

[numPQ,temp1]=deconv(numPQ,poly(rcom));
[denPQ,temp2]=deconv(denPQ,poly(rcom));

% Clean up any imaginary residues after "deconvolution"

numPQ=real(numPQ);
denPQ=real(denPQ);