MODULE My_Module
  CONTAINS
    SUBROUTINE Cubic(a, b, c, d, root) !https://en.m.wikipedia.org/wiki/Cubic_equation#General_cubic_formula
      USE, INTRINSIC :: ISO_Fortran_env, dp=>REAL64 !modern DOUBLE PRECISION
      IMPLICIT NONE
      INTEGER :: k, m
      REAL(dp), INTENT(IN) :: a, b, c, d
      COMPLEX(dp), INTENT(IN OUT), DIMENSION(:), CONTIGUOUS :: root !dummy array for Cardano
      REAL(dp), PARAMETER :: ONE_THIRD = 1.0_dp / 3.0_dp
      COMPLEX(dp) :: Delta_0, Delta_1, radical, CC, x
      COMPLEX(dp), PARAMETER :: xi = CMPLX(-1.0_dp / 2.0_dp, SQRT(3.0_dp) / 2.0_dp, dp) !primitive cube root of unity
      Delta_0 = b**2 - 3.0_dp * a * c
      Delta_1 = 2.0_dp * b**3 - 9.0_dp * a * b * c + 27.0_dp * a**2 * d
      radical = SQRT(Delta_1**2 - 4.0_dp * Delta_0**3)
      CC = ((Delta_1 + radical) / 2.0_dp)**ONE_THIRD
      m = 1
      roots: DO k = 0, 2 !three roots
        x = -(b + xi**k * CC + Delta_0 / (xi**k * CC)) / (3.0_dp * a)
        root(m) = x
        !print *, m, root(m)
        m = m + 1
      END DO roots
    END SUBROUTINE Cubic
END MODULE My_Module

PROGRAM Quartic
  USE My_Module
  USE, INTRINSIC :: ISO_Fortran_env, dp=>REAL64 !modern DOUBLE PRECISION
  IMPLICIT NONE
  INTEGER :: i
  REAL(dp) :: alpha, beta, gamma, delta, epsilon !general quartic
  REAL(dp) :: p, q, r !depressed quartic
  REAL(dp) :: aa, bb, cc, dd !resolvent cubic
  REAL(dp) :: s_1 !depressed cubic
  COMPLEX(dp), DIMENSION(3) :: Cardano
  COMPLEX(dp) :: radical, x_1, x_2, x_3, x_4
  PRINT *, " "
  WRITE(6, '(a6)', advance='no') "  a = "
  READ *, alpha !https://math.stackexchange.com/a/57688/536410 (Equation 1)
  IF (alpha > 0.0_dp .OR. alpha < 0.0_dp) THEN !a /= 0
    GO TO 1
  ELSE
    PRINT *, " "
    STOP
  END IF
1 WRITE(6, '(a6)', advance='no') "  b = "
  READ *, beta !https://math.stackexchange.com/a/57688/536410 (Equation 1)
  WRITE(6, '(a6)', advance='no') "  c = "
  READ *, gamma !https://math.stackexchange.com/a/57688/536410 (Equation 1)
  WRITE(6, '(a6)', advance='no') "  d = "
  READ *, delta !https://math.stackexchange.com/a/57688/536410 (Equation 1)
  WRITE(6, '(a6)', advance='no') "  e = "
  READ *, epsilon !https://math.stackexchange.com/a/57688/536410 (Equation 1)
  PRINT *, " "
  !https://en.wikipedia.org/wiki/Quartic_function#Converting_to_a_depressed_quartic
  p = (8.0_dp * gamma * alpha - 3.0_dp * beta**2) / (8.0_dp * alpha**2)
  q = (beta**3 - 4.0_dp * gamma * beta * alpha + 8.0_dp * delta * alpha**2) / (8.0_dp * alpha**3)
  r = (-3.0_dp * beta**4 + 256.0_dp * epsilon * alpha**3 - 64.0_dp * delta * beta * alpha**2 + 16.0_dp * gamma * beta**2 * alpha)  &
    / (256.0_dp * alpha**4)
  !https://math.stackexchange.com/a/57688/536410 (Equation 4)
  aa = 8.0_dp
  bb = -4.0_dp * p
  cc = -8.0_dp * r
  dd = 4.0_dp * p * r - q**2
  Cardano = CMPLX(0.0_dp, 0.0_dp, dp) !initialize
  CALL Cubic(aa, bb, cc, dd, Cardano)
  i = 1 !initialize
  test: DO !There Can Be Only One
2   IF (INT(Cardano(i) / CONJG(Cardano(i))) == 1) THEN !are you for real?
      s_1 = REALPART(Cardano(i))
    ELSE
      i = i + 1 !no good, try again
      GO TO 2
    END IF
    radical = SQRT(2.0_dp * s_1 - p)
    IF (ISNAN(REALPART(radical))) THEN !when all three roots are real
      i = i + 1 !no good, try again
      GO TO 2
    END IF
    EXIT
  END DO test
  !print *, i, Cardano(i)
  !https://math.stackexchange.com/a/57688/536410
  x_1 = -0.5_dp * (radical - SQRT(-2.0_dp * s_1 - p + (2.0_dp * q) / radical)) - beta / (4.0_dp * alpha)
  IF (INT(x_1 / CONJG(x_1)) == 1) THEN !are you for real?
    WRITE(6, '(f20.15)') REALPART(x_1)
  ELSE
    WRITE(6, '(2f20.15, a1)') REALPART(x_1), IMAGPART(x_1), "i"
  END IF
  x_2 = -0.5_dp * (radical + SQRT(-2.0_dp * s_1 - p + (2.0_dp * q) / radical)) - beta / (4.0_dp * alpha)
  IF (INT(x_2 / CONJG(x_2)) == 1) THEN !are you for real?
    WRITE(6, '(f20.15)') REALPART(x_2)
  ELSE
    WRITE(6, '(2f20.15, a1)') REALPART(x_2), IMAGPART(x_2), "i"
  END IF
  x_3 =  0.5_dp * (radical + SQRT(-2.0_dp * s_1 - p - (2.0_dp * q) / radical)) - beta / (4.0_dp * alpha)
  IF (INT(x_3 / CONJG(x_3)) == 1) THEN !are you for real?
    WRITE(6, '(f20.15)') REALPART(x_3)
  ELSE
    WRITE(6, '(2f20.15, a1)') REALPART(x_3), IMAGPART(x_3), "i"
  END IF
  x_4 =  0.5_dp * (radical - SQRT(-2.0_dp * s_1 - p - (2.0_dp * q) / radical)) - beta / (4.0_dp * alpha)
  IF (INT(x_4 / CONJG(x_4)) == 1) THEN !are you for real?
    WRITE(6, '(f20.15)') REALPART(x_4)
  ELSE
    WRITE(6, '(2f20.15, a1)') REALPART(x_4), IMAGPART(x_4), "i"
  END IF
  PRINT *, " "
END PROGRAM Quartic