function [ c3, c4 ] = poly_divide ( c1, c2 )

%*****************************************************************************80
%
%% poly_divide computes the quotient and remainder of two polynomials.
%
%  Discussion:
%
%    The power sum form of a polynomial is:
%
%      p(x) = c(0) + c(1) * x + ... + c(n-1) * x^(n-1) + c(n) * x^(n)
%
%    Because MATLAB doesn't allow 0 indexing, we use an adjusted formula:
%
%      p(x) = c(1) + c(2) * x + ... + c(n) * x^(n-1) + c(n+1) * x^(n)
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license. 
%
%  Modified:
%
%    09 November 2019
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    real C1(1:D1+1): the coefficients of the numerator.
%
%    real C2(1:D2+1): the coefficients of the divisor.
%
%  Output:
%
%    real C3(1:D3+1): the quotient.
%
%    real C4(1:D4+1): the remainder.
%
  d1 = length ( c1 ) - 1;
  d2 = length ( c2 ) - 1;

  d3 = d1 - d2;
  d4 = d2 - 1;

  c3 = zeros ( d3 + 1, 1 );
  for i = d3 : -1 : 0
    c3(i+1) = c1(i+d2+1) / c2(d2+1);
    c1(i+1:i+d2) = c1(i+1:i+d2) - c3(i+1) * c2(1:d2);
  end

  c4(1:d4+1) = c1(1:d4+1);

  return
end
