! Program Example from "Essential Fortran 90 & 95" by Loren P. Meissner
! Copyright 1996. Copying for sale requires permission from the author.
! Otherwise, distribution is permitted if these three lines are included.

! Example 3.10. Find all roots of a quadratic equation, given coefficients A, B, C.
  program D10
    implicit none
    real :: A, B, C, Part_1, Part_2, Discriminant, First_Root
    complex :: C_Root
! start program D10
    read (unit = *, fmt = *) A, B, C
    write (unit = *, fmt = *) " Coefficients are: ", A, B, C
    if (A /= 0) then
  ! A is not zero. Find discriminant.
      Discriminant = B ** 2 - 4.0 * A * C
      if (Discriminant >= 0.0) then   ! .. Roots are real.
        Part_1 = -B / (2.0 * A)
        ! .. For best results, force sign of Part_2 to agree with sign of Part_1.
        Part_2 = sign( sqrt( Discriminant ) / (2.0 * A), Part_1 )
        First_Root = Part_1 + Part_2
        ! .. For best results, unless First_Root is zero, compute second root as
        !    C / (A * First_Root) rather than as Part_1 - Part_2.
        if (First_Root /= 0.0) then
          write (unit = *, fmt = *) " Real roots: ", First_Root, C / (A * First_Root)
        else
          write (unit = *, fmt = *) " Real roots: ", First_Root, Part_1 - Part_2
        end if
      else   ! .. Discriminant < 0; roots are complex.
        C_Root = (-B + sqrt( cmplx( Discriminant ))) / (2.0 * A)
        write (unit = *, fmt = *) " Complex roots: ", C_Root, conjg( C_Root )
      end if
    else   ! .. A is zero. Equation is not a full quadratic.
      if (B /= 0) then
        write (unit = *, fmt = *) " Not a quadratic.  Root is: ", -C / B
      else
        write (unit = *, fmt = *) " Equation has the form C = 0; C is: ", C
      end if
    end if
    stop
  end program D10
