c ----------------------------------------------------------------------
c Solve a cubic equation where a, b, c, and d are real.
c   a*x**3 + b*x**2 + c*x + d = 0
c
c Variables used:
c   a, b, c, d  ... coefficients (input)
c   x(i)        ... three (generally) complex solutions (output)
c   nroot       ... number of roots
c
c Progamming Note: same as CUBICEQ1.FOR, except with subroutines
c ----------------------------------------------------------------------
c Instructor: Nam Sun Wang
c ----------------------------------------------------------------------

c Declare variables
      complex x(3)

c Print the program header ---------------------------------------------
      write(*,*) '-----------------------------------------------------'
      write(*,*) 'Solve a cubic equation:'
      write(*,*) '   a*x**3 + b*x**2 + c*x + d = 0'
      write(*,*) 'where a, b, c, and d are real.'
      write(*,*) '-----------------------------------------------------'

c Input the coefficients -----------------------------------------------
      write(*,*) 'Enter a: '
      read(*,*) a
      write(*,*) 'Enter b: '
      read(*,*) b
      write(*,*) 'Enter c: '
      read(*,*) c
      write(*,*) 'Enter d: '
      read(*,*) d

c Call a routine to solve the cubic equation ---------------------------
      call cubic(a, b, c, d, x, nroot)

c Output the results. --------------------------------------------------
      write(*,*)
      write(*,'(a,i2,a)') ' There are', nroot, ' root(s):'
      do i=1, nroot
        write(*,*) x(i)
      enddo

      end



c ----------------------------------------------------------------------
      subroutine quad(a, b, c, x, nroot)
c ----------------------------------------------------------------------
c Solve a quadratic equation where a, b, and c are real.
c   a*x*x + b*x + c = 0
c
c Public Variables
c   a, b, c     ... coefficients (input)
c   x(i)        ... two complex solutions (output)
c   nroot       ... number of roots (output)
c ----------------------------------------------------------------------

      complex x(1)

      if(a .eq. 0.)then
        if(b .eq. 0.)then
c         We have a non-equation; therefore, we have no valid solution
          nroot = 0
        else
c         We have a linear equation with 1 root.
          nroot = 1
          x(1) = cmplx(-c/b, 0.)
        endif
      else
c     We have a true quadratic equation.  Apply the quadratic formula to find two roots.
      nroot = 2
        DD = b*b-4.*a*c
        if(DD .ge. 0.)then
          x(1) = cmplx((-b+sqrt(DD))/2./a, 0.)
          x(2) = cmplx((-b-sqrt(DD))/2./a, 0.)
        else
          x(1) = cmplx(-b/2./a, +sqrt(-DD)/2./a)
          x(2) = cmplx(-b/2./a, -sqrt(-DD)/2./a)
        endif
      endif

      end



c ----------------------------------------------------------------------
      subroutine cubic(a, b, c, d, x, nroot)
c ----------------------------------------------------------------------
c Solve a cubic equation where a, b, c, and d are real.
c   a*x**3 + b*x**2 + c*x + d = 0
c
c Public Variables
c   a, b, c, d  ... coefficients (input)
c   x(i)        ... three (generally) complex solutions (output)
c   nroot       ... number of roots (output)
c Local Variables:
c   y1, y2, y3  ... three transformed solutions
c
c Formula used are given in Tuma, "Engineering Mathematics Handbook", p7
c   (McGraw Hill, 1978).
c   Step 0: If a is 0. use the quadratic formula to avoid dividing by 0.
c   Step 1: Calculate p and q
c           p = ( 3*c/a - (b/a)**2 ) / 3
c           q = ( 2*(b/a)**3 - 9*b*c/a/a + 27*d/a ) / 27
c   Step 2: Calculate discriminant D
c           D = (p/3)**3 + (q/2)**2
c   Step 3: Depending on the sign of D, we follow different strategy.
c           If D<0, three distinct real roots.
c           If D=0, three real roots of which at least two are equal.
c           If D>0, one real and two complex roots.
c   Step 3a: For D>0 and D=0,
c           Calculate u and v
c           u = cubic_root(-q/2 + sqrt(D))
c           v = cubic_root(-q/2 - sqrt(D))
c           Find the three transformed roots
c           y1 = u + v
c           y2 = -(u+v)/2 + i (u-v)*sqrt(3)/2
c           y3 = -(u+v)/2 - i (u-v)*sqrt(3)/2
c   Step 3b Alternately, for D<0, a trigonometric formulation is more convenient
c           y1 =  2 * sqrt(|p|/3) * cos(phi/3)
c           y2 = -2 * sqrt(|p|/3) * cos((phi+pi)/3)
c           y3 = -2 * sqrt(|p|/3) * cos((phi-pi)/3)
c           where phi = acos(-q/2/sqrt(|p|**3/27))
c                 pi  = 3.141592654...
c   Step 4  Finally, find the three roots
c           x = y - b/a/3
c ----------------------------------------------------------------------

c Declare variables
      complex x(1)
      data pi/3.141592654/

c Step 0: If a is 0 use the quadratic formula. -------------------------
      if(a .eq. 0.)then
        call quad(b, c, d, x, nroot)
        return
      endif

c Cubic equation with 3 roots
      nroot = 3

c Step 1: Calculate p and q --------------------------------------------
      p  = c/a - b*b/a/a/3.
      q  = (2.*b*b*b/a/a/a - 9.*b*c/a/a + 27.*d/a) / 27.

c Step 2: Calculate DD (discriminant) ----------------------------------
      DD = p*p*p/27. + q*q/4.

c Step 3: Branch to different algorithms based on DD -------------------
      if(DD .lt. 0.)then
c       Step 3b:
c       3 real unequal roots -- use the trigonometric formulation
        phi = acos(-q/2./sqrt(abs(p*p*p)/27.))
        temp1=2.*sqrt(abs(p)/3.)
        y1 =  temp1*cos(phi/3.)
        y2 = -temp1*cos((phi+pi)/3.)
        y3 = -temp1*cos((phi-pi)/3.)
      else
c       Step 3a:
c       1 real root & 2 conjugate complex roots OR 3 real roots (some are equal)
        temp1 = -q/2. + sqrt(DD)
        temp2 = -q/2. - sqrt(DD)
        u = abs(temp1)**(1./3.)
        v = abs(temp2)**(1./3.)
        if(temp1 .lt. 0.) u=-u
        if(temp2 .lt. 0.) v=-v
        y1  = u + v
        y2r = -(u+v)/2.
        y2i =  (u-v)*sqrt(3.)/2.
      endif

c Step 4: Final transformation -----------------------------------------
      temp1 = b/a/3.
      y1 = y1-temp1
      y2 = y2-temp1
      y3 = y3-temp1
      y2r=y2r-temp1

c Assign answers -------------------------------------------------------
      if(DD .lt. 0.)then
        x(1) = cmplx( y1,  0.)
        x(2) = cmplx( y2,  0.)
        x(3) = cmplx( y3,  0.)
      elseif(DD .eq. 0.)then
        x(1) = cmplx( y1,  0.)
        x(2) = cmplx(y2r,  0.)
        x(3) = cmplx(y2r,  0.)
      else
        x(1) = cmplx( y1,  0.)
        x(2) = cmplx(y2r, y2i)
        x(3) = cmplx(y2r,-y2i)
      endif

      end