******************************************************************************
**      This program calculates all roots (real and/or complex) to an       **
**      N:th-degree polynomial with complex coefficients. (N <= 10)         **
**                                                                          **
**                 n        n-1                                             **
**       P(z) = a z  + a   z    + ... + a z + a                             **
**               n      n-1              1     0                            **
**                                                                          **
**      The execution terminates if                                         **
**                                                                          **
**                 1) Abs (Z1-Z0) < EPS       ==>   Root found = Z1         **
**                 2) ITER > MAX              ==>   Slow convergence        **
**                                                                          **
**      The program sets EPS to 1.0E-7 and MAX to 30                        **
**                                                                          **
**      The NEWTON-RAPHSON method is used:   z1 = z0 - P(z0)/P'(z0)         **
**      The values of P(z0) and P'(z0) are calculated with HORNER'S SCHEME. **
**                                                                          **
**      The array A(0:10) contains the complex coefficients of the          ** 
**      polynomial P(z).                                                    **
**      The array B(1:10) contains the complex coefficients of the          **
**      polynomial Q(z),                                                    **
**                  where P(Z) = (z-z0)*Q(z) + P(z0)                        **
**      The coefficients to Q(z) are calculated with HORNER'S SCHEME.       **
**                                                                          **
**      When the first root has been obtained with the NEWTON-RAPHSON 	    **
**      method, it is divided away (deflation).                             **
**      The quotient polynomial is Q(z).                                    **
**      The process is repeated, now using the coefficients of Q(z) as      **
**      input data.                                                         **
**      As a starting approximation the value STARTV = 1+i is used          **
**      in all cases.                                                       **
**      Z0 is the previous approximation to the root.                       **
**      Z1 is the latest approximation to the root.                         **
**      F0 = P(Z0)                                                          **
**      FPRIM0 = P'(Z0)                                                     **
******************************************************************************
      COMPLEX      A(0:10), B(1:10), Z0, Z1, STARTV
      INTEGER      N, I, ITER, MAX
      REAL         EPS
      DATA  EPS/1E-7/, MAX /30/, STARTV /(1,1)/
******************************************************************************
20    WRITE(6,*)  'Give the degree of the polynomial'
      READ (5,*) N
      IF (N .GT. 10) THEN
         WRITE(6,*) 'The polynomial degree is maximum 10'
      GOTO 20
      WRITE (6,*) 'Give the coefficients of the polynomial,',
     , ' as complex constants'
      WRITE (6,*) 'Highest degree coefficient first'
      DO 30 I = N, 0, -1
          WRITE (6,*) 'A(' , I, ') = '
          READ  (5,*)  A(I)
30    CONTINUE
      WRITE (5,*) '    The roots are','        Number of iterations'
******************************************************************************
40    IF (N GT 0) THEN
C     ******     Find the next root                                     ******
          Z0 = (0,0)
          ITER = 0
          Z1 = STARTV
50        IF (ABS(Z1-Z0) .GE. EPS) THEN 
C         ++++++ Continue the iteration until two estimates             ++++++
C         ++++++ are sufficiently close.                                ++++++
              ITER = ITER + 1
              IF (ITER .GT. MAX) THEN
C             ------ Too many iterations  ==> TERMINATE                 ------
                     WRITE (6,*) 'Too many iterations.'
                     WRITE (6,*) 'The latest approximation to the root is ',Z1
              GOTO 200
              ENDIF
              Z0 = Z1
              HORNER (N, A, B, Z0, F0, FPRIM0)
C             ++++++   NEWTON-RAPHSON'S METHOD                          ++++++
              Z1 = Z0 - F0/FPRIM0
          GOTO 50
          ENDIF
C         ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
100       WRITE (6,*) Z1, '         ',Iter
C         ****** The root is found. Divide it away and seek the next one *****
          N = N - 1
          FOR I = 0 TO N DO
              A(I) = B(I+1)
      GOTO 40
      ENDIF
200   END

      SUBROUTINE  HORNER(N, A, B, Z, F, FPRIM)
******    For the meaning of the parameters - see the main program      ******
******    BI and CI are temporary variable.                             ******

      INTEGER N, I
      COMPLEX A(1:10), B(0:10), Z, F, FPRIM, BI, CI

      BI = A(N)
      B(N) = BI
      CI = BI

      DO 60  I = N-1, 1, -1
             BI = A(I) + Z*BI
             CI = BI + Z*CI
             B(I) = BI
C            ++++++ BI = B(I) for the calculation of P(Z)               ++++++
C            ++++++ CI  for the calculation of P'(Z)                    ++++++
60    CONTINUE

      F = A(0) + Z*BI
      FPRIM = CI

      RETURN
      END
******   END HORNER'S SCHEME                                            ******

******   This program is composed by Ulla Ouchterlony 1984              ******
******   The program is translated by Bo Einarsson 1995                 ******
