C     ALGORITHM 419 COLLECTED ALGORITHMS FROM ACM.
C     ALGORITHM APPEARED IN COMM. ACM, VOL. 15, NO. 02,
C     P. 097.
C
C Added changes from Remark on Algorithm 419 by David H. Withers
C CACM (March 1974) Vol 17 No 3 p. 157 
      PROGRAM CPOLYDR
C
C              DRIVER TO TEST CPOLY
C
C
      DIMENSION MULT(50)
      LOGICAL FLAG(50)
      REAL*8 FZR(50),FZI(50),BND(50)
      LOGICAL FAIL
      DOUBLE PRECISION  P(50),PI(50),ZR(50),ZI(50)
      WRITE(6,100)
  100 FORMAT('1EXAMPLE 1.  POLYNOMIAL WITH ZEROS 1,2,...,10.')
      P(1)=1
      P(2)=-55
      P(3)=1320
      P(4)=-18150
      P(5)=157773
      P(6)=-902055
      P(7) = 3416930
      P(8)=-8409500
      P(9)=12753576
      P(10)=-10628640
      P(11)=3628800
      DO 10 I=1,11
   10 PI(I)=0
      CALL PRTC(11,P,PI)
      CALL CPOLY(P,PI,10,ZR,ZI,FAIL)
      IF(FAIL) GO TO 95
      CALL PRTZ (10,ZR,ZI)
    2 WRITE(6,101)
  101 FORMAT('1EXAMPLE 2. ZEROS ON IMAGINARY AXIS DEGREE 3.')
      P(1)=1
      P(2)=0
      P(3)=-10001.0001D0
      P(4)=0
      PI(1)=0
      PI(2)=-10001.0001D0
      PI(3)=0
      PI(4)=1
      CALL PRTC(4,P,PI)
      CALL CPOLY(P,PI,3,ZR,ZI,FAIL)
      IF (FAIL) GO TO 96
      CALL PRTZ (3,ZR,ZI)
    3 WRITE(6,102)
  102 FORMAT('1EXAMPLE 3. ZEROS AT 1+I,1/2*(1+I)....1/(2**-9)*(1+I)')
      P(1)=1.0
      P(2)=-1.998046875
      P(3)=0.0
      P(4)=.7567065954208374D0
      P(5)=-.2002119533717632D0
      P(6)=1.271507365163416D-2
      P(7)=0
      P(8)=-1.154642632172909D-5
      P(9)=1.584803612786345D-7
      P(10)=-4.652065399568528D-10
      P(11)=0
      PI(1)=0
      PI(2)=P(2)
      PI(3)=2.658859252929688D0
      PI(4)=-7.567065954208374D-1
      PI(5)=0
      PI(6)=P(6)
      PI(7)=-7.820779428584501D-4
      PI(8)=-P(8)
      PI(9)=0
      PI(10)=P(10)
      PI(11)=9.094947017729282D-13
      CALL PRTC(11,P,PI)
      CALL CPOLY(P,PI,10,ZR,ZI,FAIL)
      IF (FAIL) GO TO 97
      CALL PRTZ(10,ZR,ZI)
    4 WRITE(6,103)
  103 FORMAT('1EXAMPLE 4. MULTIPLE ZEROS')
      P(1)=1
      P(2)=-10
      P(3)=3
      P(4)=284
      P(5)=-1293
      P(6)=2374
      P(7)=-1587
      P(8)=-920
      P(9)=2204
      P(10)=-1344
      P(11)=288
      PI(1)=0
      PI(2)=-10
      PI(3)=100
      PI(4)=-334
      PI(5)=200
      PI(6)=1394
      PI(7) =-3836
      PI(8)=4334
      PI(9)=-2352
      PI(10)=504
      PI(11)=0
      CALL PRTC(11,P,PI)
      CALL CPOLY(P,PI,10,ZR,ZI,FAIL)
      IF (FAIL) GO TO 98
      CALL PRTZ(10,ZR,ZI)
    5 WRITE(6,104)
  104 FORMAT('1EXAMPLE 5. 12 ZEROS EVENLY DISTRIBUTE ON A CIRCLE OF RADI
     -US 1. CENTERED AT 0+2I.')
      P(1)=1
      P(2)=0
      P(3)=-264
      P(4)=0
      P(5)=7920
      P(6)=0
      P(7)=-59136
      P(8)=0
      P(9)=126720
      P(10)=0
      P(11)=-67584
      P(12)=0
      P(13)=4095
      PI(1)=0
      PI(2)=-24
      PI(3)=0
      PI(4)=1760
      PI(5)=0
      PI(6)=-25344
      PI(7)=0
      PI(8)=101376
      PI(9)=0
      PI(10)=-112640
      PI(11)=0
      PI(12)=24576
      PI(13)=0
      CALL PRTC(13,P,PI)
      CALL CPOLY(P,PI,12,ZR,ZI,FAIL)
      IF(FAIL) GO TO 99
      CALL PRTZ(12,ZR,ZI)
      RETURN
   95 WRITE(6,105)
      GO TO 2
   96 WRITE(6,105)
      GO TO 3
   97 WRITE(6,105)
      GO TO 4
   98 WRITE(6,105)
      GO TO 5
   99 WRITE(6,105)
      RETURN
  105 FORMAT(//' CPOLY HAS FAILED ON THIS EXAMPLE')
      END
      SUBROUTINE PRTC(N,P,Q)
      DOUBLE PRECISION P(50),Q(50)
      WRITE(6,10) (P(I),Q(I) ,I=1,N)
   10 FORMAT(//' COEFFICIENTS' /50(2D26.16/))
      RETURN
      END
      SUBROUTINE PRTZ(N,ZR,ZI)
      DOUBLE PRECISION ZR(50),ZI(50)
      WRITE(6,10) (ZR(I),ZI(I) ,I=1,N)
   10 FORMAT(//' ZEROS'/ 50(2D26.16/))
      RETURN
      END
      SUBROUTINE CPOLY(OPR,OPI,DEGREE,ZEROR,ZEROI,FAIL)                 CPOL  10
C
C Added changes from Remark on Algorithm 419 by David H. Withers
C CACM (March 1974) Vol 17 No 3 p. 157 
C
C FINDS THE ZEROS OF A COMPLEX POLYNOMIAL.
C OPR, OPI  -  DOUBLE PRECISION VECTORS OF REAL AND
C IMAGINARY PARTS OF THE COEFFICIENTS IN
C ORDER OF DECREASING POWERS.
C DEGREE    -  INTEGER DEGREE OF POLYNOMIAL.
C ZEROR, ZEROI  -  OUTPUT DOUBLE PRECISION VECTORS OF
C REAL AND IMAGINARY PARTS OF THE ZEROS.
C FAIL      -  OUTPUT LOGICAL PARAMETER,  TRUE  ONLY IF
C LEADING COEFFICIENT IS ZERO OR IF CPOLY
C HAS FOUND FEWER THAN DEGREE ZEROS.
C THE PROGRAM HAS BEEN WRITTEN TO REDUCE THE CHANCE OF OVERFLOW
C OCCURRING. IF IT DOES OCCUR, THERE IS STILL A POSSIBILITY THAT
C THE ZEROFINDER WILL WORK PROVIDED THE OVERFLOWED QUANTITY IS
C REPLACED BY A LARGE NUMBER.
C COMMON AREA
      COMMON/GLOBAL/PR,PI,HR,HI,QPR,QPI,QHR,QHI,SHR,SHI,
     *    SR,SI,TR,TI,PVR,PVI,ARE,MRE,ETA,INFIN,NN
      DOUBLE PRECISION SR,SI,TR,TI,PVR,PVI,ARE,MRE,ETA,INFIN,
     *    PR(50),PI(50),HR(50),HI(50),QPR(50),QPI(50),QHR(50),
     *    QHI(50),SHR(50),SHI(50)
C TO CHANGE THE SIZE OF POLYNOMIALS WHICH CAN BE SOLVED, REPLACE
C THE DIMENSION OF THE ARRAYS IN THE COMMON AREA.
      DOUBLE PRECISION XX,YY,COSR,SINR,SMALNO,BASE,XXX,ZR,ZI,BND,
     *    OPR(1),OPI(1),ZEROR(1),ZEROI(1),
     *    CMOD,SCALE,CAUCHY,DSQRT
      LOGICAL FAIL,CONV
      INTEGER DEGREE,CNT1,CNT2
C INITIALIZATION OF CONSTANTS
      CALL MCON(ETA,INFIN,SMALNO,BASE)
      ARE = ETA
      MRE = 2.0D0*DSQRT(2.0D0)*ETA
      XX = .70710678
      YY = -XX
      COSR = -.069756474
      SINR = .99756405
      FAIL = .FALSE.
      NN = DEGREE+1
C ALGORITHM FAILS IF THE LEADING COEFFICIENT IS ZERO.
      IF (OPR(1) .NE. 0.0D0 .OR. OPI(1) .NE. 0.0D0) GO TO 10
          FAIL = .TRUE.
          RETURN
C REMOVE THE ZEROS AT THE ORIGIN IF ANY.
   10 IF (OPR(NN) .NE. 0.0D0 .OR. OPI(NN) .NE. 0.0D0) GO TO 20
          IDNN2 = DEGREE-NN+2
          ZEROR(IDNN2) = 0.0D0
          ZEROI(IDNN2) = 0.0D0
          NN = NN-1
          GO TO 10
C MAKE A COPY OF THE COEFFICIENTS.
   20 DO 30 I = 1,NN
          PR(I) = OPR(I)
          PI(I) = OPI(I)
          SHR(I) = CMOD(PR(I),PI(I))
   30 CONTINUE
C SCALE THE POLYNOMIAL.
      BND = SCALE (NN,SHR,ETA,INFIN,SMALNO,BASE)
      IF (BND .EQ. 1.0D0) GO TO 40
      DO 35 I = 1,NN
          PR(I) = BND*PR(I)
          PI(I) = BND*PI(I)
   35 CONTINUE
C START THE ALGORITHM FOR ONE ZERO .
   40 IF (NN.GT. 2) GO TO 50
C CALCULATE THE FINAL ZERO AND RETURN.
          CALL CDIVID(-PR(2),-PI(2),PR(1),PI(1),ZEROR(DEGREE),
     *    ZEROI(DEGREE))
          RETURN
C CALCULATE BND, A LOWER BOUND ON THE MODULUS OF THE ZEROS.
   50 DO 60 I = 1,NN
          SHR(I) = CMOD(PR(I),PI(I))
   60 CONTINUE
      BND = CAUCHY(NN,SHR,SHI)
C OUTER LOOP TO CONTROL 2 MAJOR PASSES WITH DIFFERENT SEQUENCES
C OF SHIFTS.
      DO 100 CNT1 = 1,2
C FIRST STAGE CALCULATION, NO SHIFT.
          CALL NOSHFT(5)
C INNER LOOP TO SELECT A SHIFT.
          DO 90 CNT2 = 1,9
C SHIFT IS CHOSEN WITH MODULUS BND AND AMPLITUDE ROTATED BY
C 94 DEGREES FROM THE PREVIOUS SHIFT
               XXX = COSR*XX-SINR*YY
               YY = SINR*XX+COSR*YY
               XX = XXX
               SR = BND*XX
               SI = BND*YY
C SECOND STAGE CALCULATION, FIXED SHIFT.
               CALL FXSHFT(10*CNT2,ZR,ZI,CONV)
               IF (.NOT. CONV) GO TO 80
C THE SECOND STAGE JUMPS DIRECTLY TO THE THIRD STAGE ITERATION.
C IF SUCCESSFUL THE ZERO IS STORED AND THE POLYNOMIAL DEFLATED.
                    IDNN2 = DEGREE-NN+2
                    ZEROR(IDNN2) = ZR
                    ZEROI(IDNN2) = ZI
                    NN = NN-1
                    DO 70 I = 1,NN
                         PR(I) = QPR(I)
                         PI(I) = QPI(I)
   70               CONTINUE
                    GO TO 40
   80          CONTINUE
C IF THE ITERATION IS UNSUCCESSFUL ANOTHER SHIFT IS CHOSEN.
   90     CONTINUE
C IF 9 SHIFTS FAIL, THE OUTER LOOP IS REPEATED WITH ANOTHER
C SEQUENCE OF SHIFTS.
  100 CONTINUE
C THE ZEROFINDER HAS FAILED ON TWO MAJOR PASSES.
C RETURN EMPTY HANDED.
      FAIL = .TRUE.
      RETURN
      END
      SUBROUTINE  NOSHFT(L1)                                            NOSH1130
C COMPUTES  THE DERIVATIVE  POLYNOMIAL AS THE INITIAL H
C POLYNOMIAL AND COMPUTES L1 NO-SHIFT H POLYNOMIALS.
C COMMON AREA
      COMMON/GLOBAL/PR,PI,HR,HI,QPR,QPI,QHR,QHI,SHR,SHI,
     *    SR,SI,TR,TI,PVR,PVI,ARE,MRE,ETA,INFIN,NN
      DOUBLE PRECISION SR,SI,TR,TI,PVR,PVI,ARE,MRE,ETA,INFIN,
     *    PR(50),PI(50),HR(50),HI(50),QPR(50),QPI(50),QHR(50),
     *    QHI(50),SHR(50),SHI(50)
      DOUBLE PRECISION XNI,T1,T2,CMOD
      N = NN-1
      NM1 = N-1
      DO 10 I = 1,N
          XNI = NN-I
          HR(I) = XNI*PR(I)/FLOAT(N)
          HI(I) = XNI*PI(I)/FLOAT(N)
   10 CONTINUE
      DO 50 JJ = 1,L1
          IF (CMOD(HR(N),HI(N)) .LE. ETA*10.0D0*CMOD(PR(N),PI(N)))
     *    GO TO 30
          CALL CDIVID(-PR(NN),-PI(NN),HR(N),HI(N),TR,TI)
          DO 20 I = 1,NM1
               J = NN-I
               T1 = HR(J-1)
               T2 = HI(J-1)
               HR(J) = TR*T1-TI*T2+PR(J)
               HI(J) = TR*T2+TI*T1+PI(J)
   20     CONTINUE
          HR(1) = PR(1)
          HI(1) = PI(1)
          GO TO 50
C IF THE CONSTANT TERM IS ESSENTIALLY ZERO, SHIFT H COEFFICIENTS.
   30     DO 40 I = 1,NM1
               J = NN-I
               HR(J) = HR(J-1)
               HI(J) = HI(J-1)
   40     CONTINUE
          HR(1) = 0.0D0
          HI(1) = 0.0D0
   50 CONTINUE
      RETURN
      END
      SUBROUTINE FXSHFT(L2,ZR,ZI,CONV)                                  FXSH1550
C COMPUTES L2 FIXED-SHIFT H POLYNOMIALS AND TESTS FOR
C CONVERGENCE.
C INITIATES A VARIABLE-SHIFT ITERATION AND RETURNS WITH THE
C APPROXIMATE ZERO IF SUCCESSFUL.
C L2 - LIMIT OF FIXED SHIFT STEPS
C ZR,ZI - APPROXIMATE ZERO IF CONV IS .TRUE.
C CONV  - LOGICAL INDICATING CONVERGENCE OF STAGE 3 ITERATION
C COMMON AREA
      COMMON/GLOBAL/PR,PI,HR,HI,QPR,QPI,QHR,QHI,SHR,SHI,
     *    SR,SI,TR,TI,PVR,PVI,ARE,MRE,ETA,INFIN,NN
      DOUBLE PRECISION SR,SI,TR,TI,PVR,PVI,ARE,MRE,ETA,INFIN,
     *    PR(50),PI(50),HR(50),HI(50),QPR(50),QPI(50),QHR(50),
     *    QHI(50),SHR(50),SHI(50)
      DOUBLE PRECISION ZR,ZI,OTR,OTI,SVSR,SVSI,CMOD
          LOGICAL CONV,TEST,PASD,BOOL
      N = NN-1
C EVALUATE P AT S.
      CALL POLYEV(NN,SR,SI,PR,PI,QPR,QPI,PVR,PVI)
      TEST = .TRUE.
      PASD = .FALSE.
C CALCULATE FIRST T = -P(S)/H(S).
      CALL CALCT(BOOL)
C MAIN LOOP FOR ONE SECOND STAGE STEP.
      DO 50 J = 1,L2
          OTR = TR
          OTI = TI
C COMPUTE NEXT H POLYNOMIAL AND NEW T.
          CALL NEXTH(BOOL)
          CALL CALCT(BOOL)
          ZR = SR+TR
          ZI = SI+TI
C TEST FOR CONVERGENCE UNLESS STAGE 3 HAS FAILED ONCE OR THIS
C IS THE LAST H POLYNOMIAL .
          IF ( BOOL .OR. .NOT. TEST .OR. J .EQ. L2) GO TO 50
          IF (CMOD(TR-OTR,TI-OTI) .GE. .5D0*CMOD(ZR,ZI)) GO TO 40
               IF (.NOT. PASD) GO TO 30
C THE WEAK CONVERGENCE TEST HAS BEEN PASSED TWICE, START THE
C THIRD STAGE ITERATION, AFTER SAVING THE CURRENT H POLYNOMIAL
C AND SHIFT.
                    DO 10 I = 1,N
                         SHR(I) = HR(I)
                         SHI(I) = HI(I)
   10               CONTINUE
                    SVSR = SR
                    SVSI = SI
                    CALL VRSHFT(10,ZR,ZI,CONV)
                    IF (CONV) RETURN
C THE ITERATION FAILED TO CONVERGE. TURN OFF TESTING AND RESTORE
C H,S,PV AND T.
                    TEST = .FALSE.
                    DO 20 I = 1,N
                         HR(I) = SHR(I)
                         HI(I) = SHI(I)
   20               CONTINUE
                    SR = SVSR
                    SI = SVSI
                    CALL POLYEV(NN,SR,SI,PR,PI,QPR,QPI,PVR,PVI)
                    CALL CALCT(BOOL)
                    GO TO 50
   30          PASD = .TRUE.
               GO TO 50
   40     PASD = .FALSE.
   50 CONTINUE
C ATTEMPT AN ITERATION WITH FINAL H POLYNOMIAL FROM SECOND STAGE.
      CALL VRSHFT(10,ZR,ZI,CONV)
      RETURN
      END
      SUBROUTINE VRSHFT(L3,ZR,ZI,CONV)                                  VRSH2230
C CARRIES OUT THE THIRD STAGE ITERATION.
C L3 - LIMIT OF STEPS IN STAGE 3.
C ZR,ZI   - ON ENTRY CONTAINS THE INITIAL ITERATE, IF THE
C ITERATION CONVERGES IT CONTAINS THE FINAL ITERATE
C ON EXIT.
C CONV    -  .TRUE. IF ITERATION CONVERGES
C COMMON AREA
      COMMON/GLOBAL/PR,PI,HR,HI,QPR,QPI,QHR,QHI,SHR,SHI,
     *    SR,SI,TR,TI,PVR,PVI,ARE,MRE,ETA,INFIN,NN
      DOUBLE PRECISION SR,SI,TR,TI,PVR,PVI,ARE,MRE,ETA,INFIN,
     *    PR(50),PI(50),HR(50),HI(50),QPR(50),QPI(50),QHR(50),
     *    QHI(50),SHR(50),SHI(50)
      DOUBLE PRECISION ZR,ZI,MP,MS,OMP,RELSTP,R1,R2,CMOD,DSQRT,ERREV,TP
      LOGICAL CONV,B,BOOL
      CONV = .FALSE.
      B = .FALSE.
      SR = ZR
      SI = ZI
C MAIN LOOP FOR STAGE THREE
      DO 60 I = 1,L3
C EVALUATE P AT S AND TEST FOR CONVERGENCE.
          CALL POLYEV(NN,SR,SI,PR,PI,QPR,QPI,PVR,PVI)
          MP = CMOD(PVR,PVI)
          MS = CMOD(SR,SI)
          IF (MP .GT. 20.0D0*ERREV(NN,QPR,QPI,MS,MP,ARE,MRE))
     *       GO TO 10
C POLYNOMIAL VALUE IS SMALLER IN VALUE THAN A BOUND ON THE ERROR
C IN EVALUATING P, TERMINATE THE ITERATION.
               CONV = .TRUE.
               ZR = SR
               ZI = SI
               RETURN
   10     IF (I .EQ. 1) GO TO 40
               IF (B .OR. MP .LT.OMP .OR. RELSTP .GE. .05D0)
     *            GO TO 30
C ITERATION HAS STALLED. PROBABLY A CLUSTER OF ZEROS. DO 5 FIXED
C SHIFT STEPS INTO THE CLUSTER TO FORCE ONE ZERO TO DOMINATE.
                    TP = RELSTP
                    B = .TRUE.
                    IF (RELSTP .LT. ETA) TP = ETA
                    R1 = DSQRT(TP)
                    R2 = SR*(1.0D0+R1)-SI*R1
                    SI = SR*R1+SI*(1.0D0+R1)
                    SR = R2
                    CALL POLYEV(NN,SR,SI,PR,PI,QPR,QPI,PVR,PVI)
                    DO 20 J = 1,5
                         CALL CALCT(BOOL)
                         CALL NEXTH(BOOL)
   20               CONTINUE
      OMP = INFIN
                    GO TO 50
C EXIT IF POLYNOMIAL VALUE INCREASES SIGNIFICANTLY.
   30          IF (MP*.1D0 .GT. OMP) RETURN
   40     OMP = MP
C CALCULATE NEXT ITERATE.
   50     CALL CALCT(BOOL)
          CALL NEXTH(BOOL)
          CALL CALCT(BOOL)
          IF (BOOL) GO TO 60
          RELSTP = CMOD(TR,TI)/CMOD(SR,SI)
          SR = SR+TR
          SI = SI+TI
   60 CONTINUE
      RETURN
      END
      SUBROUTINE CALCT(BOOL)                                            CALC2890
C COMPUTES  T = -P(S)/H(S).
C BOOL   - LOGICAL, SET TRUE IF H(S) IS ESSENTIALLY ZERO.
C COMMON AREA
      COMMON/GLOBAL/PR,PI,HR,HI,QPR,QPI,QHR,QHI,SHR,SHI,
     *    SR,SI,TR,TI,PVR,PVI,ARE,MRE,ETA,INFIN,NN
      DOUBLE PRECISION SR,SI,TR,TI,PVR,PVI,ARE,MRE,ETA,INFIN,
     *    PR(50),PI(50),HR(50),HI(50),QPR(50),QPI(50),QHR(50),
     *    QHI(50),SHR(50),SHI(50)
      DOUBLE PRECISION HVR,HVI,CMOD
      LOGICAL BOOL
      N = NN-1
C EVALUATE H(S).
      CALL POLYEV(N,SR,SI,HR,HI,QHR,QHI,HVR,HVI)
      BOOL = CMOD(HVR,HVI) .LE. ARE*10.0D0*CMOD(HR(N),HI(N))
      IF (BOOL) GO TO 10
          CALL CDIVID(-PVR,-PVI,HVR,HVI,TR,TI)
          RETURN
   10 TR = 0.0D0
      TI = 0.0D0
      RETURN
      END
      SUBROUTINE NEXTH(BOOL)                                            NEXT3110
C CALCULATES THE NEXT SHIFTED H POLYNOMIAL.
C BOOL   -  LOGICAL, IF .TRUE. H(S) IS ESSENTIALLY ZERO
C COMMON AREA
      COMMON/GLOBAL/PR,PI,HR,HI,QPR,QPI,QHR,QHI,SHR,SHI,
     *    SR,SI,TR,TI,PVR,PVI,ARE,MRE,ETA,INFIN,NN
      DOUBLE PRECISION SR,SI,TR,TI,PVR,PVI,ARE,MRE,ETA,INFIN,
     *    PR(50),PI(50),HR(50),HI(50),QPR(50),QPI(50),QHR(50),
     *    QHI(50),SHR(50),SHI(50)
      DOUBLE PRECISION T1,T2
      LOGICAL BOOL
      N = NN-1
      NM1 = N-1
      IF (BOOL) GO TO 20
          DO 10 J = 2,N
               T1 = QHR(J-1)
               T2 = QHI(J-1)
               HR(J) = TR*T1-TI*T2+QPR(J)
               HI(J) = TR*T2+TI*T1+QPI(J)
   10     CONTINUE
          HR(1) = QPR(1)
          HI(1) = QPI(1)
          RETURN
C IF H(S) IS ZERO REPLACE H WITH QH.
   20 DO 30 J = 2,N
          HR(J) = QHR(J-1)
          HI(J) = QHI(J-1)
   30 CONTINUE
      HR(1) = 0.0D0
      HI(1) = 0.0D0
      RETURN
      END
      SUBROUTINE POLYEV(NN,SR,SI,PR,PI,QR,QI,PVR,PVI)                   POLY3430
C EVALUATES A POLYNOMIAL  P  AT  S  BY THE HORNER RECURRENCE
C PLACING THE PARTIAL SUMS IN Q AND THE COMPUTED VALUE IN PV.
      DOUBLE PRECISION PR(NN),PI(NN),QR(NN),QI(NN),
     *    SR,SI,PVR,PVI,T
      QR(1) = PR(1)
      QI(1) = PI(1)
      PVR = QR(1)
      PVI = QI(1)
      DO 10 I = 2,NN
          T = PVR*SR-PVI*SI+PR(I)
          PVI = PVR*SI+PVI*SR+PI(I)
          PVR = T
          QR(I) = PVR
          QI(I) = PVI
   10 CONTINUE
      RETURN
      END
      DOUBLE PRECISION FUNCTION ERREV(NN,QR,QI,MS,MP,ARE,MRE)           ERRE3610
C BOUNDS THE ERROR IN EVALUATING THE POLYNOMIAL BY THE HORNER
C RECURRENCE.
C QR,QI - THE PARTIAL SUMS
C MS    -MODULUS OF THE POINT
C MP    -MODULUS OF POLYNOMIAL VALUE
C ARE, MRE -ERROR BOUNDS ON COMPLEX ADDITION AND MULTIPLICATION
      DOUBLE PRECISION QR(NN),QI(NN),MS,MP,ARE,MRE,E,CMOD
      E = CMOD(QR(1),QI(1))*MRE/(ARE+MRE)
      DO 10 I = 1,NN
          E = E*MS+CMOD(QR(I),QI(I))
   10 CONTINUE
      ERREV = E*(ARE+MRE)-MP*MRE
      RETURN
      END
      DOUBLE PRECISION FUNCTION CAUCHY(NN,PT,Q)                         CAUC3760
C CAUCHY COMPUTES A LOWER BOUND ON THE MODULI OF THE ZEROS OF A
C POLYNOMIAL - PT IS THE MODULUS OF THE COEFFICIENTS.
      DOUBLE PRECISION Q(NN),PT(NN),X,XM,F,DX,DF,
     *   DABS,DEXP,DLOG
      PT(NN) = -PT(NN)
C COMPUTE UPPER ESTIMATE OF BOUND.
      N = NN-1
      X = DEXP( (DLOG(-PT(NN)) - DLOG(PT(1)))/FLOAT(N) )
      IF (PT(N).EQ.0.0D0) GO TO 20
C IF NEWTON STEP AT THE ORIGIN IS BETTER, USE IT.
          XM = -PT(NN)/PT(N)
          IF (XM.LT.X) X=XM
C CHOP THE INTERVAL (0,X) UNITL F LE 0.
   20 XM = X*.1D0
      F = PT(1)
      DO 30 I = 2,NN
          F = F*XM+PT(I)
   30 CONTINUE
      IF (F.LE. 0.0D0) GO TO 40
          X = XM
          GO TO 20
   40 DX = X
C DO NEWTON ITERATION UNTIL X CONVERGES TO TWO DECIMAL PLACES.
   50 IF (DABS(DX/X) .LE. .005D0) GO TO 70
          Q(1) = PT(1)
          DO 60 I = 2,NN
               Q(I) = Q(I-1)*X+PT(I)
   60     CONTINUE
          F = Q(NN)
          DF = Q(1)
          DO 65 I = 2,N
               DF = DF*X+Q(I)
   65     CONTINUE
          DX = F/DF
          X = X-DX
          GO TO 50
   70 CAUCHY = X
      RETURN
      END
      DOUBLE PRECISION FUNCTION SCALE(NN,PT,ETA,INFIN,SMALNO,BASE)      SCAL4160
C RETURNS A SCALE FACTOR TO MULTIPLY THE COEFFICIENTS OF THE
C POLYNOMIAL. THE SCALING IS DONE TO AVOID OVERFLOW AND TO AVOID
C UNDETECTED UNDERFLOW INTERFERING WITH THE CONVERGENCE
C CRITERION.  THE FACTOR IS A POWER OF THE BASE.
C PT - MODULUS OF COEFFICIENTS OF P
C ETA,INFIN,SMALNO,BASE - CONSTANTS DESCRIBING THE
C FLOATING POINT ARITHMETIC.
      DOUBLE PRECISION PT(NN),ETA,INFIN,SMALNO,BASE,HI,LO,
     *    MAX,MIN,X,SC,DSQRT,DLOG
C FIND LARGEST AND SMALLEST MODULI OF COEFFICIENTS.
      HI = DSQRT(INFIN)
      LO = SMALNO/ETA
      MAX = 0.0D0
      MIN = INFIN
      DO 10 I = 1,NN
          X = PT(I)
          IF (X .GT. MAX) MAX = X
          IF (X .NE. 0.0D0 .AND. X.LT.MIN) MIN = X
   10 CONTINUE
C SCALE ONLY IF THERE ARE VERY LARGE OR VERY SMALL COMPONENTS.
      SCALE = 1.0D0
      IF (MIN .GE. LO .AND. MAX .LE. HI) RETURN
      X = LO/MIN
      IF (X .GT. 1.0D0) GO TO 20
          SC = 1.0D0/(DSQRT(MAX)*DSQRT(MIN))
          GO TO 30
   20 SC = X
      IF (INFIN/SC .GT. MAX) SC = 1.0D0
   30 L = DLOG(SC)/DLOG(BASE) + .500
      SCALE = BASE**L
      RETURN
      END
      SUBROUTINE CDIVID(AR,AI,BR,BI,CR,CI)                              CDIV4490
C COMPLEX DIVISION C = A/B, AVOIDING OVERFLOW.
      DOUBLE PRECISION AR,AI,BR,BI,CR,CI,R,D,T,INFIN,DABS
      IF (BR .NE. 0.0D0  .OR. BI .NE. 0.0D0) GO TO 10
C DIVISION BY ZERO, C = INFINITY.
          CALL MCON (T,INFIN,T,T)
          CR = INFIN
          CI = INFIN
          RETURN
   10 IF (DABS(BR) .GE. DABS(BI)) GO TO 20
          R = BR/BI
          D = BI+R*BR
          CR = (AR*R+AI)/D
          CI = (AI*R-AR)/D
          RETURN
   20 R = BI/BR
      D = BR+R*BI
      CR = (AR+AI*R)/D
      CI = (AI-AR*R)/D
      RETURN
      END
      DOUBLE PRECISION FUNCTION CMOD(R,I)                               CMOD4700
C MODULUS OF A COMPLEX NUMBER AVOIDING OVERFLOW.
      DOUBLE PRECISION R,I,AR,AI,DABS,DSQRT
      AR = DABS(R)
      AI = DABS(I)
      IF (AR .GE. AI) GO TO 10
          CMOD = AI*DSQRT(1.0D0+(AR/AI)**2)
          RETURN
   10 IF (AR .LE. AI) GO TO 20
          CMOD = AR*DSQRT(1.0D0+(AI/AR)**2)
          RETURN
   20 CMOD = AR*DSQRT(2.0D0)
      RETURN
      END
      SUBROUTINE MCON(ETA,INFINY,SMALNO,BASE)                           MCON4840
C MCON PROVIDES MACHINE CONSTANTS USED IN VARIOUS PARTS OF THE
C PROGRAM. THE USER MAY EITHER SET THEM DIRECTLY OR USE THE
C STATEMENTS BELOW TO COMPUTE THEM. THE MEANING OF THE FOUR
C CONSTANTS ARE -
C ETA       THE MAXIMUM RELATIVE REPRESENTATION ERROR
C WHICH CAN BE DESCRIBED AS THE SMALLEST POSITIVE
C FLOATING-POINT NUMBER SUCH THAT 1.0D0 + ETA IS
C GREATER THAN 1.0D0.
C INFINY    THE LARGEST FLOATING-POINT NUMBER
C SMALNO    THE SMALLEST POSITIVE FLOATING-POINT NUMBER
C BASE      THE BASE OF THE FLOATING-POINT NUMBER SYSTEM USED
C LET T BE THE NUMBER OF BASE-DIGITS IN EACH FLOATING-POINT
C NUMBER(DOUBLE PRECISION). THEN ETA IS EITHER .5*B**(1-T)
C OR B**(1-T) DEPENDING ON WHETHER ROUNDING OR TRUNCATION
C IS USED.
C LET M BE THE LARGEST EXPONENT AND N THE SMALLEST EXPONENT
C IN THE NUMBER SYSTEM. THEN INFINY IS (1-BASE**(-T))*BASE**M
C AND SMALNO IS BASE**N.
C THE VALUES FOR BASE,T,M,N BELOW CORRESPOND TO THE IBM/360.
      DOUBLE PRECISION ETA,INFINY,SMALNO,BASE
      INTEGER M,N,T
      BASE = 16.0D0
      T = 14
      M = 63
      N = -65
      ETA = BASE**(1-T)
      INFINY = BASE*(1.0D0-BASE**(-T))*BASE**(M-1)
      SMALNO = (BASE**(N+3))/BASE**3
      RETURN
      END