function [ d3, c3, d4, c4 ] = r8poly_division ( d1, c1, d2, c2 )

%*****************************************************************************80
%
%% r8poly_division() 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 MIT license. 
%
%  Modified:
%
%    24 October 2019
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    integer D1: the degree of the numerator.
%
%    real C1(1:D1+1): the coefficients of the numerator.
%
%    integer D2: the degree of the divisor.
%
%    real C2(1:D2+1): the coefficients of the divisor.
%
%  Output:
%
%    integer D3: the degree of Q.
%    If the divisor polynomial is zero, D3 is returned as -1.
%
%    real C3(1:D3+1): the quotient.
%
%    integer D4: the degree of R.
%    If the divisor polynomial is zero, D4 is returned as -1.
%
%    real C4(1:D4+1): the remainder.
%
  d1 = r8poly_degree ( d1, c1 );
  d2 = r8poly_degree ( d2, c2 );

  if ( c2(d2) == 0.0 )
    d3 = -1;
    d4 = -1;
    return
  end

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

  for i = d3 : -1 : 0
    c3(i+1) = c1(i+d2+1) / c2(d2+1);
    c1(i+d2+1) = 0.0;
    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
