function [ disc, p, r, d0, d1 ] = r8poly4_discriminant ( ...
  a, b, c, d, e )

%*****************************************************************************80
%
%% r8poly4_discriminant() returns the discriminant of a quartic polynomial.
%
%  Discussion:
%
%    The polynomial has the form
%
%      a * x^4 + b * x^3 + c * x^2 + d * x + e = 0
%
%    If the discriminant is negative, then the polynomial has two distinct
%    real roots and two distinct complex conjugate roots.
%
%    If the discriminant is positive, then all four roots are real, 
%    or all four roots are strictly complex.
%
%    If the discriminant is positive, let P=8*a*c-3*b^2
%    and D = 64*a^3*e - 16*a^2*c^2 + 16*a*b^2*c - 16*a^2*b*d - 3*b^4.
%    If P<0 and D<0, then all four roots are real and distinct.
%    If P>0 or D>0, there are two pairs of non-real complex conjugate roots
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    14 July 2025
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    John D Cook,
%    When a cubic or quartic has a double root,
%    https://www.johndcook.com/blog/2022/11/16/cubic-quartic-double-root/
%    Posted 16 November 2022.
%
%  Input:
%
%    real A, B, C, D, E: the coefficients of the polynomial.
%    A must not be zero.
%
%  Output:
%
%    real disc: the discriminant.
%
  disc = ...
      256.0 * a^3                   * e^3 ...
    - 192.0 * a^2 * b         * d   * e^2 ...
    - 128.0 * a^2       * c^2       * e^2 ...
    + 144.0 * a^2       * c   * d^2 * e   ...
    -  27.0 * a^2             * d^4       ...
    + 144.0 * a   * b^2 * c         * e^2 ...
    -   6.0 * a   * b^2       * d^2 * e   ...
    -  80.0 * a   * b   * c^2 * d   * e   ...
    +  18.0 * a   * b   * c   * d^3       ...
    +  16.0 * a         * c^4       * e   ...
    -   4.0 * a         * c^3 * d^2       ...
    -  27.0       * b^4             * e^2 ...
    +  18.0       * b^3 * c   * d   * e   ...
    -   4.0       * b^3       * d^3       ...
    -   4.0       * b^2 * c^3       * e   ...
    +               b^2 * c^2 * d^2;

  p = 8.0 * a * c - 3.0 * b^2;
  r = b^3 + 8.0 * d * a^2 - 4.0 * a * b * c;
  d0 = c^2 - 3.0 * b * d + 12.0 * a * e;
  d1 = ...
      64.0 * a^3 * e ...
    - 16.0 * a^2 * c^2 ...
    + 16.0 * a   * b^2 * c ...
    - 16.0 * a^2 * b * d ...
    - 3.0 * b^4;

  return
end
