PROGRAM QUBIC_ROOT_COMPUTATION
C=======================================================================
C                  SOLUTION OF X**3+A*X**2+B*X+C=0
C                           EIJI FUKUMORI
C                          29 DECEMBER 2012
C=======================================================================
      IMPLICIT REAL*8 ( A-H , O-Z )
      COMPLEX*16 BETA, T1, T2, ALPHA, X1, X2, X3
      WRITE (*,*) 'QUBIC ROOT COMPUTATION'
C------------ CASE OF TRIPLE ROOT --------------------
C    TRY  A =  -9.D0     B =  27.D0    C = -27.D0
C                  X1=3, X2=3, X3=3
C-----------------------------------------------------
C------------ CASE OF DOUBLE ROOT --------------------
C    TRY  A =  -7.D0     B =  16.D0    C = -12.D0
C                  X1=3, X2=2, X3=2
C-----------------------------------------------------
      A = -7.D0
      B = 16.D0
      C = -12.D0
C
      WRITE (*,*) 'A, B, C =', A, B, C
      IF ( A.EQ.0.D0 .AND. B.EQ.0.D0 .AND. C.EQ.0.D0 ) THEN
       X1 = DCMPLX(0.D0, 0.D0)
       X2 = DCMPLX(0.D0, 0.D0)
       X3 = DCMPLX(0.D0, 0.D0)
      ELSE
       P = -A**2/3.D0 + B
       Q = 2.D0*A**3/27.D0 -A*B/3.D0 + C
         IF ( P.EQ.0.D0 .AND. Q.EQ.0.D0 ) THEN
            X1 = DCMPLX(-A/3.D0, 0.D0)
            X2 = DCMPLX(-A/3.D0, 0.D0)
            X3 = DCMPLX(-A/3.D0, 0.D0)
         ELSE
            DISCR = Q**2 + (4.D0/27.D0)*P**3
            BETA = (-DCMPLX(Q) - CDSQRT(DCMPLX(DISCR)))/2.0D0
            WRITE (*,*) 'BETA =', BETA
            ALPHA = BETA**(1.D0/3.D0)
            T1 = DCMPLX(-0.5D0, +0.5D0*DSQRT(3.D0))
            T2 = DCMPLX(-0.5D0, -0.5D0*DSQRT(3.D0))
            WRITE (*,*) 'ALPHAS ARE ', ALPHA, ALPHA*T1, ALPHA*T2
            X1 = ALPHA    -DCMPLX(P)/(3.D0*ALPHA)    - DCMPLX(A)/3.D0
            X2 = ALPHA*T1 -DCMPLX(P)/(3.D0*ALPHA*T1) - DCMPLX(A)/3.D0
            X3 = ALPHA*T2 -DCMPLX(P)/(3.D0*ALPHA*T2) - DCMPLX(A)/3.D0
         END IF
      END IF
C
      WRITE (*,*) 'SOLUTIONS TO QUBIC EQ. ARE ',  X1, X2, X3
      WRITE (*,*) 'CONFIRMATION OF SOLUTIONS'
      WRITE (*,*) 'X1**3+A*X1**2+B*X1+C=', X1**3+A*X1**2+B*X1+C
      WRITE (*,*) 'X2**3+A*X2**2+B*X2+C=', X2**3+A*X2**2+B*X2+C
      WRITE (*,*) 'X3**3+A*X3**2+B*X3+C=', X3**3+A*X3**2+B*X3+C
      OPEN ( 1,FILE='QUBIC.SOL', STATUS='UNKNOWN' )
      WRITE (1,*) '0. 0.'
      WRITE (1,*) DREAL(X1), DIMAG(X1)
      WRITE (1,*)
      WRITE (1,*) '0. 0.'
      WRITE (1,*) DREAL(X2), DIMAG(X2)
      WRITE (1,*)
      WRITE (1,*) '0. 0.'
      WRITE (1,*) DREAL(X3), DIMAG(X3)
      R =            DSQRT(DREAL(X1)**2+DIMAG(X1)**2)
      R = DMAX1 ( R, DSQRT(DREAL(X2)**2+DIMAG(X2)**2) )
      R = DMAX1 ( R, DSQRT(DREAL(X3)**2+DIMAG(X3)**2) )
      PI = 4.D0*DATAN (1.D0)
      N = 40
      DANG = 2.D0*PI / N
      WRITE (1,*)
      WRITE (1,*) R, ' 0.'
      DO I = 1 , N
      ANG = I*DANG
      X = R*DCOS(ANG)
      Y = R*DSIN(ANG)
      WRITE (1,*) X , Y
      END DO
      CLOSE (1)
      STOP
      END