function r = polynomial_resultant_roots ( p1, p2 )

%*****************************************************************************80
%
%% polynomial_resultant_roots() computes the resultant of polynomials p1 and p2.
%
%  Discussion:
%
%    This algorithm computes the resultant using knowledge of the roots
%    of the polynomials.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    30 January 2024
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Michael Pohst, Hans Zassenhaus,
%    Algorithmic Algebraic Number Theory,
%    Cambridge University Press, 1989
%    LC: QA247.P58
%    ISBN: 0-521-33060-2
%
%  Input:
%
%    real p1(np1), p2(np2): two polynomials.
%    On input, the first entry is the constant coefficient.
%
%  Output:
%
%    real r: the resultant.
%

%
%  Reverse order of polynomial coefficients to satisfy stupid MATLAB convention.
%
  p1 = p1(end:-1:1);
  p2 = p2(end:-1:1);
%
%  Get roots of polynomials.
%
  r1 = roots ( p1 );
  r2 = roots ( p2 );
%
%  Retrieve highest order coefficients.
%
  a1 = p1(1);
  a2 = p2(1);
%
%  Determine number of roots.
%
  n1 = length ( p1 ) - 1;
  n2 = length ( p2 ) - 1;
%
%  Compute resultant.
%
  r = 1.0;
  for i = 1 : n1
    for j = 1 : n2
      r = r * ( r1(i) - r2(j) );
    end
  end

  r = r * a1^n2 * a2^n1;

  return
end

