Code covered by the BSD License

# Polynomial division by convolution - quotient and reminder

### Feng Cheng Chang (view profile)

12 Apr 2008 (Updated )

Division of two polynomials to get quotient and reminder using convolution matrix.

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;
```