C
C      ________________________________________________________
C     |                                                        |
C     | USE JENKENS-TRAUB ALGORITHM TO COMPUTE ZEROS OF A MONIC|
C     |  POLYNOMIAL P(1) + P(2)X +...+ P(ND)X**(ND-1) + X**ND  |
C     |                                                        |
C     |    INPUT:                                              |
C     |                                                        |
C     |         P     --THE FIRST ND COEFFICIENTS OF POLYNOMIAL|
C     |                 (COMPLEX DOUBLE PRECISION ARRAY)       |
C     |                                                        |
C     |         ND    --DEGREE OF POLYNOMIAL                   |
C     |                                                        |
C     |         W     --WORK ARRAY (REQUIRED LENGTH DEPENDS ON |
C     |                 THE DIFFICULTY IN COMPUTING ZEROS. IN  |
C     |                 MANY CASES, BETWEEN 4ND+32 AND 6ND     |
C     |                 COMPLEX DOUBLE PRECISION ENTRIES ARE   |
C     |                 SUFFICIENT)                            |
C     |                                                        |
C     |    OUTPUT:                                             |
C     |                                                        |
C     |         Z     --ZEROS (COMPLEX DOUBLE PRECISION ARRAY) |
C     |                                                        |
C     |         W(1)  --NUMBER OF COMPUTED ZEROS               |
C     |                                                        |
C     |    NOTE:                                               |
C     |                                                        |
C     |      THIS ALGORITHM CAN GENERATE UNDERFLOWS SINCE THE  |
C     |      FUNCTION VALUES TEND TO ZERO AS THE ITERATIONS    |
C     |      APPROACH A ROOT. THUS THE APPROPRIATE COMPILER    |
C     |      OPTIONS MUST BE INVOKED SO THAT AN UNDERFLOW IS   |
C     |      REPLACED BY ZERO AND THE EXECUTION CONTINUES. FOR |
C     |      IBM FORTRAN, TRY THE FOLLOWING STATEMENT:         |
C     |               CALL ERRSET(208,256,-1,-1,0)             |
C     |                                                        |
C     |    BUILTIN FUNCTIONS: CDABS,CDEXP,DACOS,DCMPLX,DLOG    |
C     |    PACKAGE FUNCTIONS: DMAG,DREAL,DIMAG                 |
C     |    PACKAGE SUBROUTINES: CC,CEF,ERR,FF,JT,QAD,RAD,RFN   |
C     |________________________________________________________|
C
      SUBROUTINE CZERO(Z,P,ND,W)
      COMPLEX*16 P(1),W(1),Z(1),PI,R,S,U,V,X
      DOUBLE PRECISION A,B,C,D,E,DMAG,DREAL,S1,T,T1,T2,V1,Y,Y0
      INTEGER F,G,H,I,J,K,L,M,N,NC,ND,NZ,N1,O,O1,O2,O3,O4,Q
      N = ND
      W(1) = N
      DO 10 I = 1,N
10         Z(I) = P(I)
20    IF ( N .GT. 1 ) GOTO 25
      Z(1) = -P(1)
      RETURN
25    IF ( N .GT. 2 ) GOTO 30
      CALL QAD(Z(2),Z(1))
      RETURN
C     --------------------------------
C     |*** REMOVE ZEROS AT ORIGIN ***|
C     --------------------------------
30    IF ( DMAG(Z(1)) .NE. 0.D0 ) GOTO 50
      DO 40 I = 2,N
40         Z(I-1) = Z(I)
      Z(N) = (0.D0,0.D0)
      N = N - 1
      GOTO 20
50    NZ = N
      E = 1.D0
C     ---------------------------------
C     |*** COMPUTE MACHINE EPSILON ***|
C     ---------------------------------
60    E = .5D0*E
      A = 1.D0 + E
      IF( A .GT. 1.D0 ) GOTO 60
      E = E + E
      M = 5
      PI = DCMPLX(0.D0,2.D0*DACOS(-1.D0))
C     ----------------------------------------------------
C     |*** SOME COMPILERS USE DARCOS INSTEAD OF DACOS ***|
C     ----------------------------------------------------
      X = (1.D0,0.D0)
      Y = 65536.
      Y0 = 1./Y
70    N1 = N + 1
      O1 = N1 + N
      O2 = O1 + N
      O3 = O2 + N1/2
      O4 = O3 + N1/2
      O = O4 - 1
      CALL CC(Z,W(O1),W(O2),N,Y,Y0)
      G = O4 + N1
      NC = 32*N
      K = 10
      L = 8
      Q = 0
      J = 0
      F = 0
80    F = F + 1
      Q = Q + J
      IF ( F .LT. 6 ) GOTO 90
      IF ( L .GT. NC ) GOTO 370
90    L = L + L
C     --------------------------------------------------------------
C     |*** ESTIMATE RADIUS OF SMALLEST CIRCLE CONTAINING A ZERO ***|
C     --------------------------------------------------------------
      CALL RAD(Z,W,W(N1),K,N,A)
      CALL CEF(W(O4),Z,W(O1),W(O2),W(O3),A,X,N,N1,T,T1,Y,Y0)
      T = N*T*E
      J = 1
      I = N1
      H = I
100   I = I/2
      J = J + J
      IF ( I .GT. 0 ) GOTO 100
      IF ( H+H .EQ. J ) GOTO 120
      H = J
      J = H + O
      DO 110 I = G,J
110        W(I) = (0.D0,0.D0)
120   J = L + O4
      CALL FF(W(O4),W(J),H,L)
      C = DMAG(W(O4))
      J = 0
      H = L - 1
C     ----------------------------------
C     |*** DETERMINE STARTING GUESS ***|
C     ----------------------------------
      DO 130 I=1,H
           D = DMAG(W(O4+I))
           IF ( D .GT. C ) GOTO 130
           J = I
           C = D
130   CONTINUE
      S = A*X*CDEXP(J*PI/L)
      S1 = CDABS(S)
      IF ( S1 .LT. 1. ) S1 = 1.
      X = X*CDEXP(PI/(4*L))
      X = X/CDABS(X)
C     ----------------------------------
C     |*** JENKINS-TRAUB ITERATIONS ***|
C     ----------------------------------
      DO 140 I=1,M
140        CALL JT(Z,W,N,V,S,S1)
      C = 0.D0
      J = 0
150   CALL JT(Z,W,N,V,S,S1)
      V1 = DMAG(V)
      IF ( V1 .EQ. 0. ) GOTO 240
      IF ( S1 .GT. 1. ) GOTO 170
      U = (1.D0,0.D0)
      I = N
160   I = I - 1
      U = W(I) + U*S
      IF ( I .GT. 1 ) GOTO 160
      IF ( DMAG(U) .EQ. 0.D0 ) GOTO 80
      GOTO 190
170   U = W(1)
      R = 1./S
      DO 180 I = 2,N
180        U = W(I) + U*R
      IF ( DMAG(U) .EQ. 0.D0 ) GOTO 80
      U = U*R
190   R = S - V/U
      J = J + 1
      D = DMAG(R-S)
      S = R
      IF ( S1 .LE. A ) GOTO 200
      CALL ERR(Z,N,S1,T2)
      T2 = N*T2*E
      IF ( V1 .LE. T2 ) GOTO 240
      GOTO 210
200   IF ( N*DLOG(T1/S1)+DLOG(T/V1) .GE. 0. ) GOTO 230
210   S1 = CDABS(S)
      IF ( S1 .LT. 1. ) S1 = 1.
      IF ( J .LT. 3 ) GOTO 220
      IF ( .5D0*D .GT. B+C ) GOTO 80
      IF ( J .LT. M ) GOTO 220
      IF ( 4.D0*D .GT. B+C ) GOTO 80
220   B = C
      C = D
      GOTO 150
230   CALL ERR(Z,N,S1,T2)
      T2 = N*E*T2
      IF ( V1 .GT. T2 ) GOTO 210
240   R = Z(N) + S
      I = N
250   I = I - 1
      V = Z(I) + R*S
      Z(I) = R
      R = V
      IF ( I .GT. 1 ) GOTO 250
      Z(N) = S
      Q = Q + J
      N = N - 1
      IF ( N .GT. 2 ) GOTO 70
      CALL QAD(Z(2),Z(1))
      NC = ND
      N = 1
260   L = N
      N = ND
      M = N + N
      J = 1
      K = N - 1
      DO 270 I = 2,K
           W(J) = J*P(I)
           W(J+N) = J*I*P(I+1)
           W(J+M) = DMAG(P(J))
270        J = I
      W(K) = K*P(N)
      W(N) = N
      W(M-1) = N*K
      W(M) = 0.D0
      W(M+K) = DMAG(P(K))
      W(M+N) = DMAG(P(N))
      O1 = M + 1
      N1 = N + 1
C     --------------------------
C     |*** REFINE THE ZEROS ***|
C     --------------------------
      DO 340 K = L,NZ
           S = Z(K)
           DO 280 J = 1,3
                CALL RFN(P,W,W(N1),W(O1),N,S,V,C,T)
                T = N*T*E
280             IF ( C .LE. T ) GOTO 330
           A = CDABS(S)
           IF ( A .GT. 1.D0 ) GOTO 300
           V = (1.D0,0.D0)
           T = 1.
           I = N
290        V = V*S + P(I)
           T = T*A + DREAL(W(I+M))
           I = I - 1
           IF ( I .GT. 0 ) GOTO 290
           GOTO 320
300        R = 1./S
           A = 1./A
           V = (0.D0,0.D0)
           T = 0.
           DO 310 I = 1,N
                V = V*R + P(I)
310             T = T*A + DREAL(W(I+M))
           V = V*R + DCMPLX(1.D0,0.D0)
           T = T*A + 1.
320        A = DMAG(V)
           IF ( A .GT. N*T*E ) GOTO 340
330        Z(K) = S
340   CONTINUE
      DO 350 I = L,N
350        W(N1-I) = Z(I)
      L = N1 - L
      DO 360 I = 1,L
360        Z(I) = W(I)
      W(1) = NC
      RETURN
370   NC = ND - N + 1
      WRITE(6,*) 'SINCE THE STOPPING CRITERION HAS NOT BEEN SATISFIED'
      WRITE(6,*) 'AFTER',Q,'ITERATIONS, WE STOP WHILE COMPUTING'
      WRITE(6,*) 'ROOT NUMBER',NC
      Z(N) = S
      GOTO 260
      END
C%
C     ------------------------------
C     |*** ROOTS OF A QUADRATIC ***|
C     ------------------------------
      SUBROUTINE QAD(B,C)
      COMPLEX*16 B,C,R,S
      DOUBLE PRECISION DMAG,T
      T = DMAG(B)
      R = DCMPLX(4.D0,0.D0)*C
      IF ( T .GT. DSQRT(DMAG(R)) ) GOTO 10
      R = CDSQRT(B*B-R)
      GOTO 20
10    S = 1.D0/T
      R = T*CDSQRT((S*B)**2-R*S*S)
20    S = R + B
      R = R - B
      T = DMAG(R)
      IF ( T .GE. DMAG(S) ) GOTO 30
      C = -(C+C)/S
      B = DCMPLX(-.5D0,0.D0)*S
      RETURN
30    IF ( T .EQ. 0.D0 ) RETURN
      C = (C+C)/R
      B = DCMPLX(.5D0,0.D0)*R
      RETURN
      END
C%
C     ---------------------------------
C     |*** JENKINS-TRAUB ITERATION ***|
C     ---------------------------------
      SUBROUTINE JT(P,Q,N,V,S,S1)
      COMPLEX*16 P(1),Q(1),R,S,T,U,V,W
      DOUBLE PRECISION DMAG,S1
      INTEGER I,M,N
      IF ( S1 .GT. 1. ) GOTO 70
      U = (0.D0,0.D0)
      V = (1.D0,0.D0)
      I = N
10    V = P(I) + V*S
      U = Q(I) + U*S
      I = I - 1
      IF ( I .GT. 0 ) GOTO 10
      IF ( DMAG(U) .NE. 0.D0 ) GOTO 50
      M = N
20    R = Q(M)
      Q(M) = (0.D0,0.D0)
      I = M
30    I = I - 1
      U = Q(I) + R*S
      Q(I) = R
      R = U
      IF ( I .GT. 1 ) GOTO 30
      U = (0.D0,0.D0)
      I = M
40    U = Q(I) + U*S
      I = I - 1
      IF ( I .GT. 0 ) GOTO 40
      M = M - 1
      IF ( DMAG(U) .EQ. 0.D0 ) GOTO 20
50    R = V/U
      T = P(N) - R*Q(N) + S
      Q(N) = (1.D0,0.D0)
      I = N
60    I = I - 1
      U = P(I) - R*Q(I) + T*S
      Q(I) = T
      T = U
      IF ( I .GT. 1 ) GOTO 60
      RETURN
70    W = 1./S
      U = Q(1)
      V = P(1)
      DO 80 I = 2,N
           U = U*W + Q(I)
80         V = V*W + P(I)
      U = U*W
      V = V*W + 1.
      IF ( DMAG(U) .NE. 0. ) GOTO 120
      M = N
90    R = Q(M)
      Q(M) = (0.D0,0.D0)
      I = M
100   I = I - 1
      U = Q(I) + R*S
      Q(I) = R
      R = U
      IF ( I .GT. 1 ) GOTO 100
      U = (0.D0,0.D0)
      M = M - 1
      DO 110 I = 1,N
110        U = U*W + Q(I)
      U = U*W
      IF ( DMAG(U) .EQ. 0. ) GOTO 90
120   R = V/U
      U = P(1) - R*Q(1)
      Q(1) = -W*U
      DO 130 I = 2,N
           U = P(I) - R*Q(I) + U*W
130        Q(I) = -W*U
      Q(N) = (1.D0,0.D0)
      RETURN
      END
C%
C     -----------------------
C     |*** A NEWTON-STEP ***|
C     -----------------------
      SUBROUTINE RFN(P,Q,R,Z,N,S,V,B,C)
      COMPLEX*16 P(1),Q(1),R(1),Z(1),M,S,T,U,V,W
      DOUBLE PRECISION A,B,C,DMAG,DREAL
      INTEGER I,N
      A = CDABS(S)
      IF ( A .GT. 1.D0 ) GOTO 20
      V = (1.D0,0.D0)
      U = (0.D0,0.D0)
      W = (0.D0,0.D0)
      C = 1.
      I = N
10    V = V*S + P(I)
      U = U*S + Q(I)
      W = W*S + R(I)
      C = C*A + DREAL(Z(I))
      I = I - 1
      IF ( I .GT. 0 ) GOTO 10
      GOTO 40
20    T = DCMPLX(1.D0,0.D0)/S
      A = 1./A
      V = P(1)
      U = Q(1)
      W = R(1)
      C = DREAL(Z(1))
      DO 30 I = 2,N
         V = V*T + P(I)
         U = U*T + Q(I)
         W = W*T + R(I)
30       C = C*A + DREAL(Z(I))
      V = V*T + DCMPLX(1.D0,0.D0)
      U = U*T
      W = W*T
      C = C*A + 1.
40    B = DMAG(V)
      IF ( B .EQ. 0.D0 ) RETURN
      A = DMAG(U)
      IF ( A .EQ. 0.D0 ) RETURN
      T = V/U
      A = DMAG(W)
      IF ( A .EQ. 0.D0 ) GOTO 50
      U = U/W
      M = U - T
      A = DMAG(M)
      IF ( A .EQ. 0.D0 ) RETURN
      M = U/M
      S = S - M*T
      RETURN
50    S = S - T
      RETURN
      END
C%
C     --------------------------------------------------------------
C     |*** ESTIMATE RADIUS OF SMALLEST CIRCLE CONTAINING A ZERO ***|
C     --------------------------------------------------------------
      SUBROUTINE RAD(P,Q,R,K,N,A)
      COMPLEX*16 P(1),Q(1),R(1),T
      DOUBLE PRECISION A,C,D,F,DMAG,S
      INTEGER I,J,K,L,M,N
      SAVE J,C,F
      DATA J/1/
      IF ( J .EQ. K ) GOTO 20
      J = 1
      L = 1
      DO 10 I = 2,N
           Q(L) = L*P(I)/N
10         L = I
      Q(N) = 1
      C = N*CDABS(P(1))
      F = 1.D0/N
      A = CDABS(P(1))**(1.D0/DFLOAT(N))
      GOTO 40
20    DO 30 I = 1,N
30         Q(I) = R(I)
40    L = J + K
50    IF ( DMAG(Q(1)) .EQ. 0.D0 ) GOTO 90
      J = J+1
      T = P(1)/Q(1)
      DO 60 I = 2,N
60         Q(I-1) = P(I) - T*Q(I)
      S = 1.D0/DFLOAT(J)
      F = F**(1.D0-S)*CDABS(T)**S
      D = F*(C/CDABS(Q(1)))**S
      IF ( D .LT. A ) A = D
70    IF ( J .LT. L ) GOTO 50
      K = J
      DO 80 I = 1,N
80         R(I) = Q(I)
      RETURN
90    M = J
100   J = J + 1
      DO 110 I = 2,N
110        Q(I-1) = Q(I)
      Q(N) = (0.D0,0.D0)
      IF ( DMAG(Q(1)) .EQ. 0.D0 ) GOTO 100
      S = 1.D0/J
      F = F**(DFLOAT(M)*S)
      D = F*(C/CDABS(Q(1)))**S
      IF ( D .LT. A ) A = D
      T = P(1)/Q(1)
      DO 120 I = 2,N
120        Q(I-1) = P(I) - T*Q(I)
      Q(N) = (1.D0,0.D0)
      F = F*CDABS(T)**S
      GOTO 70
      END
C%
C     --------------------------------
C     |*** FAST FOURIER TRANSFORM ***|
C     --------------------------------
      SUBROUTINE FF(A,B,N,L)
      COMPLEX*16 A(1),B(1),S,T,W
      DOUBLE PRECISION P
      INTEGER I,IL,I1,I2,J,K,L,L1,L2,M,M0,M1,N
      M = N
      IF ( L .GE. M ) GOTO 20
      M = M - 1
      DO 10 J = L,M,L
           DO 10 I = 1,L
10              A(I) = A(I) + A(I+J)
      M = L
20    P = DACOS(-1.D0)
C     ----------------------------------------------------
C     |*** SOME COMPILERS USE DARCOS INSTEAD OF DACOS ***|
C     ----------------------------------------------------
      W = CDEXP(M*DCMPLX(0.D0,P)/L)
      L1 = L - 1
      L2 = L/2
      IF ( M .EQ. L ) GOTO 40
      DO 30 J = M,L1,M
           DO 30 I = 1,M
30         A(I+J) = A(I)
40    M0 = M
      M = M/2
      M1 = M - 1
      IF ( M .EQ. 0 ) RETURN
      S = 1
      T = -1
      I1 = 0
      DO 60 J = 1,L,M0
           IL = J + M1
           I2 = I1 + L2
           DO 50 I = J,IL
                K = I + M
                B(I+I1) = A(I) + S*A(K)
50              B(I+I2) = A(I) + T*A(K)
           S = W*S
           T = W*T
           I1 = I1 - M
60    CONTINUE
      DO 70 I = 1,L
70         A(I) = B(I)
      W = CDSQRT(W)
      GOTO 40
      END
C%
      DOUBLE PRECISION FUNCTION DMAG(Z)
      COMPLEX*16 Z
      DOUBLE PRECISION DIMAG,DREAL
      DMAG = DABS(DREAL(Z)) + DABS(DIMAG(Z))
      RETURN
      END
      DOUBLE PRECISION FUNCTION DIMAG(A)
      DOUBLE PRECISION C(2)
      COMPLEX*16 A,B
      EQUIVALENCE (B,C)
      B = A
      DIMAG = C(2)
      RETURN
      END
      DOUBLE PRECISION FUNCTION DREAL(A)
      DOUBLE PRECISION C(2)
      COMPLEX*16 A,B
      EQUIVALENCE (B,C)
      B = A
      DREAL = C(1)
      RETURN
      END
C%
      SUBROUTINE CC(Z,F,E,N,Y,Y0)
      COMPLEX*16 F(1),Z(1),X
      DOUBLE PRECISION DMAG,E(1),S,T,Y,Y0
      INTEGER I,N
      DO 40 I = 1,N
           S = 0.
           X = Z(I)
           T = DMAG(X)
           IF ( T .GT. 1. ) GOTO 20
           IF ( T .GT. Y0 ) GOTO 30
           IF ( T .EQ. 0. ) GOTO 30
10         T = T*Y
           X = X*Y
           S = S - 1.
           IF ( T .LE. Y0 ) GOTO 10
           GOTO 30
20         T = T*Y0
           X = X*Y0
           S = S + 1.
           IF ( T .GT. 1. ) GOTO 20
30         E(I) = S
           F(I) = X
40    CONTINUE
      RETURN
      END
C%
      SUBROUTINE CEF(W,Z,F,E,D,R,X,N,N1,T,T1,Y,Y0)
      COMPLEX*16 F(1),W(1),Z(1),C,X,P
      DOUBLE PRECISION D(1),DMAG,E(1),A,B,Q,R,S,T,T1,Y,Y0
      INTEGER I,N,N1
      P = 1.
      Q = 0.
      S = E(1)
      C = X*R
      IF ( R .GT. 1. ) GOTO 50
      T = 1.
      T1 = 1.
      DO 20 I = 1,N
           A = E(I) + Q
           IF ( A .GT. S ) S = A
           W(I) = P*F(I)
           D(I) = A
           T = R*T + DMAG(Z(N1-I))
           P = P*C
           IF ( DMAG(P) .GT. Y0 ) GOTO 20
10         P = P*Y
           Q = Q - 1.
           IF ( DMAG(P) .LE. Y0 ) GOTO 10
20    CONTINUE
30    IF ( Q .GT. S ) S = Q
      DO 40 I = 1,N
           A = D(I) - S
40         W(I) = W(I)*Y**A
      A = Q - S
      W(N1) = P*Y**A
      RETURN
50    T1 = R
      B = 1./R
      T = 0.
      DO 70 I = 1,N
           A = E(I) + Q
           IF ( A .GT. S ) S = A
           W(I) = P*F(I)
           D(I) = A
           T = B*T + DMAG(Z(I))
           P = P*C
           IF ( DMAG(P) .LE. 1. ) GOTO 70
60         P = P*Y0
           Q = Q + 1.
           IF ( DMAG(P) .GT. 1. ) GOTO 60
70    CONTINUE
      T = B*T + 1.
      GOTO 30
      END
C%
      SUBROUTINE ERR(Z,N,S,T)
      COMPLEX*16 Z(1)
      DOUBLE PRECISION DMAG,R,S,T
      INTEGER I,N
      IF ( S .GT. 1. ) GOTO 20
      T = 1.
      I = N
10    T = T*S + DMAG(Z(I))
      I = I - 1
      IF ( I .GT. 0 ) GOTO 10
      RETURN
20    R = 1./S
      T = DMAG(Z(1))
      DO 30 I = 2,N
30         T = T*R + DMAG(Z(I))
      T = T*R + 1.
      RETURN
      END