polydiv_H(b,a,c)

function [q,r,qc,rc,c] = polydiv_H(b,a,c)
 %
 %   Polynomial division by convolution.
 %
 %   Given b(x) and a(x), we will find q(x) and r(x)  in
 %        b(x) = a(x)*q(x) + r(x)
 %   where
 %       b(x) = b(0)+ ...+b(k)*x^k + ... +b(n)*x^n
 %       a(x) = a(0)+ ...+a(k)*x^k + ... +a(m)*x^m
 %       q(x) = q(0)+ ...+q(k)*x^k + ... +q(n-m)*x^(n-m)
 %       r(x) = r(0)+ ...+r(k)*x^k + ... +r(m-1)*x^(m-1)
 %
 %   If coefficients of b(x) and a(x) are all integers, we may set  c = 0, 
 %   so that the process involve mostly integer arithmetric multiplications.
 %   The round-off errors may be eliminated.
 %
 %   This code is similar to the MATLAB built-in function: 'deconv.m'.
 %
 %   by F C Chang    04/12/08   Updated   03/30/09  07/07/11

  if nargin < 3,  c = 1;  end;
     n = length(b);  m = length(a);
  if m > n,  q = 0; r = b; qc = 0; rc = b;   return;  end;
  if c == 0,  w(1) = a(1)^(n-m+1);   else,  w(1) = 1; c = 1;  end;
 for k = 2:n-m+2,
     w(k) = [b(k-1),-a(min(k-1,m):-1:2)]*[w(1),w(max(2,k+1-m):k-1)]'/a(1);
 end
 for k = n-m+3:n+1,
     w(k) = [b(k-1),-a(min(m,k-1):-1:k-1+m-n)]*[w(1),w(max(k+1-m,2):n-m+2)]';
 end;
     qc = w(2:n-m+2);
     rc = w(n-m+3:n+1);
     q = qc;
     r = rc;
  if c == 0,  c = w(1);  q = qc/c; r = rc/c;   end;