subroutine cpoly ( opr, opi, degree, zeror, zeroi, fail )

c*********************************************************************72
c
cc cpoly() finds the zeros of a complex polynomial.
c
c  Discussion:
c
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
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    14 March 2024
c
c  Author:
c
c    Original Fortran77 version by Michael Jenkins, Joseph Traub.
c    This version by John Burkardt.
c
c  Reference:
c
c    Michael Jenkins, Joseph Traub,
c    Algorithm 419: Zeros of a Complex Polynomial,
c    Communications of the Association for Computing Machinery,
c    Volume 15, Number 2, February 1972, pages 97-99.
c
c  Input:
c
c    double precision opr(degree+1), opi(degree+1): the real and imaginary
c    parts of the coefficients in order of decreasing powers.
c
c    integer degree: the degree of the polynomial.
c
c  Output:
c
c    double precision zeror(degree), zeroi(degree): the real and imaginary parts 
c    of the zeros.
c
c    logical fail:  true if the leading coefficient is zero or if cpoly
c    has found fewer than degree zeros.
c
      implicit none

      integer degree

      double precision are
      double precision base
      double precision bnd
      double precision cmod
      double precision cosr
      double precision eta
      logical fail
      double precision hi(degree+1)
      double precision hr(degree+1)
      integer i
      integer idnn2
      double precision infin
      double precision mre
      integer n
      integer nn
      double precision opi(degree+1)
      double precision opr(degree+1)
      double precision pi(degree+1)
      double precision pr(degree+1)
      double precision pvi
      double precision pvr
      double precision qhi(degree+1)
      double precision qhr(degree+1)
      double precision qpi(degree+1)
      double precision qpr(degree+1)
      double precision rescale
      double precision shr(degree+1)
      double precision si
      double precision sinr
      double precision smalno
      double precision sr
      double precision ti
      double precision tr
      double precision xx
      double precision yy
      double precision zeroi(degree)
      double precision zeror(degree)
      double precision zi
      double precision zr
c
c  Set constants.
c
      call mcon ( eta, infin, smalno, base )
      are = eta
      mre = 2.0d0 * sqrt ( 2.0d0 ) * eta
      xx = 0.70710678
      yy = -xx
      cosr = -0.060756474
      sinr = 0.99756405
      fail = .false.
      nn = degree + 1
      n = degree
c
c  Algorithm fails if the leading coefficient is zero.
c
      if ( opr(1) == 0.0d0 .and. opi(1) == 0.0d0 ) then
        fail = .true.
        return
      end if
c
c  Remove the zeros at the origin if any.
c
      do while ( opr(nn) == 0.0d0 .and. opi(nn) == 0.0d0 )
        idnn2 = degree - nn + 2
        zeror(idnn2) = 0.0d0
        zeroi(idnn2) = 0.0d0
        nn = nn - 1
      end do
c
c  Copy the coefficients.
c
      do i = 1, nn
        pr(i) = opr(i)
        pi(i) = opi(i)
      end do
c
c  Rescale the polynomial.
c
      do i = 1, nn
        shr(i) = cmod ( pr(i), pi(i) )
      end do

      bnd = rescale ( nn, shr )

      do i = 1, nn
        pr(i) = bnd * pr(i)
        pi(i) = bnd * pi(i)
      end do
c
c  Start the algorithm for one zero.
c
      fail = .false.

      do while ( .not. fail )

        if ( nn <= 2 ) then
          call cdivid ( -pr(2), -pi(2), pr(1), pi(1), zeror(degree),
     &      zeroi(degree) )
          return
        end if

        call loopy ( degree, nn, hr, hi, qhr, qhi, pr, pi, qpr, 
     &    qpi, shr, pvr, pvi, sr, si, tr, ti, zr, zi, cosr, sinr, xx, 
     &    yy, are, mre, zeror, zeroi, fail )

      end do

      return
      end
      subroutine loopy ( degree, nn, hr, hi, qhr, qhi, pr, pi, qpr, 
     &  qpi, shr, pvr, pvi, sr, si, tr, ti, zr, zi, cosr, sinr, xx, 
     &  yy, are, mre, zeror, zeroi, fail )

c*********************************************************************72
c
cc loopy() applies several shifts to the polynomial to search for zeros.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    14 March 2024
c
c  Author:
c
c    Original Fortran77 version by Michael Jenkins, Joseph Traub.
c    This version by John Burkardt.
c
c  Reference:
c
c    Michael Jenkins, Joseph Traub,
c    Algorithm 419: Zeros of a Complex Polynomial,
c    Communications of the Association for Computing Machinery,
c    Volume 15, Number 2, February 1972, pages 97-99.
c
c  Input:
c
c  Output:
c
      implicit none

      integer nn

      double precision are
      double precision bnd
      double precision cauchy
      double precision cmod
      integer cnt1
      integer cnt2
      logical conv
      double precision cosr
      integer degree
      logical fail
      double precision hi(nn)
      double precision hr(nn)
      integer i
      integer idnn2
      double precision mre
      double precision pi(nn)
      double precision pr(nn)
      double precision pvi
      double precision pvr
      double precision qhi(nn)
      double precision qhr(nn)
      double precision qpi(nn)
      double precision qpr(nn)
      double precision shr(nn)
      double precision si
      double precision sinr
      double precision sr
      double precision ti
      double precision tr
      double precision xx
      double precision xxx
      double precision yy
      double precision zeroi(nn)
      double precision zeror(nn)
      double precision zi
      double precision zr

      fail = .false.
c
c  Calculate a lower bound on the modulus of the zeros.
c
      do i = 1, nn
        shr(i) = cmod ( pr(i), pi(i) )
      end do

      bnd = cauchy ( nn, shr )
c
c  Outer loop to control 2 major passes with different sequences of shifts.
c
      do cnt1 = 1, 2
c
c  First stage calculation, no shift.
c
        call noshft ( nn, 5, pr, pi, tr, ti, hr, hi )
c
c  Inner loop to select a shift.
c
        do cnt2 = 1, 9
c
c  Shift is chosen with modulus bnd and amplitude rotated by
c  94 degrees from the previous shift.
c
          xxx = cosr * xx - sinr * yy
          yy  = sinr * xx + cosr * yy
          xx = xxx
          sr = bnd * xx
          si = bnd * yy
c
c  Second stage calculation, fixed shift.
c
          call fxshft ( nn, 10 * cnt2, are, mre, pr, pi, hr, 
     &      hi, qpr, qpi, qhr, qhi, sr, si, tr, ti, pvr, pvi, 
     &      zr, zi, conv )
c
c  The second stage jumps directly to the third stage iteration.
c  If successful the zero is stored and the polynomial deflated.
c
          if ( conv ) then
            idnn2 = degree - nn + 2
            zeror(idnn2) = zr
            zeroi(idnn2) = zi
            nn = nn - 1
            do i = 1, nn
              pr(i) = qpr(i)
              pi(i) = qpi(i)
            end do
            return
          end if
c
c  If the iteration is unsuccessful another shift is chosen.
c
        end do
c
c  If 9 shifts fail, the outer loop is repeated with another
c  sequence of shifts.
c
      end do
c
c  The zerofinder has failed on two major passes.
c
      fail = .true.

      return
      end
      subroutine noshft ( nn, l1, pr, pi, tr, ti, hr, hi )

c*********************************************************************72
c
cc noshft() computes the derivative polynomial as the initial h polynomial.
c
c  Discussion:
c
c    It also computes l1 no-shift h polynomials.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    13 March 2024
c
c  Author:
c
c    Original Fortran77 version by Michael Jenkins, Joseph Traub.
c    This version by John Burkardt.
c
c  Reference:
c
c    Michael Jenkins, Joseph Traub,
c    Algorithm 419: Zeros of a Complex Polynomial,
c    Communications of the Association for Computing Machinery,
c    Volume 15, Number 2, February 1972, pages 97-99.
c
c  Input:
c
c    integer l1: the number of no-shift h polynomials to compute.
c
c  Output:
c
      implicit none

      integer nn

      double precision base
      double precision cmod
      double precision eta
      double precision hi(nn)
      double precision hr(nn)
      integer i
      double precision infin
      integer j
      integer jj
      integer l1
      integer n
      integer nm1
      double precision pi(nn)
      double precision pr(nn)
      double precision smalno
      double precision t1
      double precision t2
      double precision ti
      double precision tr
      double precision xni

      call mcon ( eta, infin, smalno, base )

      n = nn - 1
      nm1 = n - 1
      do i = 1, n
        xni = nn - i
        hr(i) = xni * pr(i) / float ( n )
        hi(i) = xni * pi(i) / float ( n )
      end do

      do jj = 1, l1

        if ( eta * 10.0d0 * cmod ( pr(n), pi(n) ) 
     &    < cmod ( hr(n), hi(n) ) ) then
 
          call cdivid ( -pr(nn), -pi(nn), hr(n), hi(n), tr, ti )

          do 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)
          end do
          hr(1) = pr(1)
          hi(1) = pi(1)
c
c  If the constant term is essentially zero, shift h coefficients.
c
        else

          do i = 1, nm1
            j = nn - i
            hr(j) = hr(j-1)
            hi(j) = hi(j-1)
          end do
          hr(1) = 0.0d0
          hi(1) = 0.0d0

        end if

      end do

      return
      end
      subroutine fxshft ( nn, l2, are, mre, pr, pi, hr, hi, qpr, qpi,
     &  qhr, qhi, sr, si, tr, ti, pvr, pvi, zr, zi, conv )

c*********************************************************************72
c
cc fxshft() computes l2 fixed-shift h polynomials and tests for convergence.
c
c  Discussion:
c
c    It initiates a variable-shift iteration and returns with the
c    approximate zero if successful.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    14 March 2024
c
c  Author:
c
c    Original Fortran77 version by Michael Jenkins, Joseph Traub.
c    This version by John Burkardt.
c
c  Reference:
c
c    Michael Jenkins, Joseph Traub,
c    Algorithm 419: Zeros of a Complex Polynomial,
c    Communications of the Association for Computing Machinery,
c    Volume 15, Number 2, February 1972, pages 97-99.
c
c  Input:
c
c    integer nn: ?
c
c    integer l2: the maximum number of fixed shift steps to take.
c
c    double precision are: error bound on complex addition.
c
c    double precision mre: error bound on complex multiplication.
c
c    double precision pr(nn), pi(nn): ?
c
c    double precision hr(nn), hi(nn): ?
c
c    double precision qpr(nn), qpi(nn): ?
c
c    double precision qhr(nn), qhi(nn): ?
c
c    double precision sr, si: ?
c
c    double precision tr, ti: ?
c
c    double precision pvr, pvi: ?
c
c  Output:
c
c    double precision zr, zi: approximate zero if conv is true.
c
c    logical conv: true if the stage 3 iteration converged.
c
      implicit none

      integer n
      integer nn

      double precision are
      logical bool
      double precision cmod
      logical conv
      double precision hi(nn)
      double precision hr(nn)
      integer i
      integer j
      integer l2
      double precision mre
      double precision oti
      double precision otr
      logical pasd
      double precision pi(nn)
      double precision pr(nn)
      double precision pvi
      double precision pvr
      double precision qhi(nn)
      double precision qhr(nn)
      double precision qpi(nn)
      double precision qpr(nn)
      double precision shi(nn)
      double precision shr(nn)
      double precision si
      double precision sr
      double precision svsi
      double precision svsr
      logical test
      double precision ti
      double precision tr
      double precision zi
      double precision zr

      n = nn - 1
c
c  Evaluate polynomial at s.
c
      call polyev ( nn, sr, si, pr, pi, qpr, qpi, pvr, pvi )
      test = .true.
      pasd = .false.
c
c  Calculate first t = -p(s)/h(s).
c
      call calct ( nn, are, bool, sr, si, hr, hi, qhr, qhi, pvr, 
     &  pvi, tr, ti )
c
c  Main loop for one second stage step.
c
      do j = 1, l2

        otr = tr
        oti = ti
c
c  Compute next h polynomial and new t.
c
        call nexth ( bool, nn, qhr, qhi, qpr, qpi, tr, ti, hr, hi )
        call calct ( nn, are, bool, sr, si, hr, hi, qhr, qhi, pvr, 
     &    pvi, tr, ti )
        zr = sr + tr
        zi = si + ti
c
c  Test for convergence unless stage 3 has failed once or this
c  is the last h polynomial .
c
        if ( bool .or. .not. test .or. j == l2 ) then
          cycle
        end if

        if ( 0.5d0 * cmod ( zr, zi ) <= 
     &               cmod ( tr - otr, ti - oti ) ) then
          pasd = .false.
          cycle
        end if

        if ( .not. pasd ) then
          pasd = .true.
          cycle
        end if
c
c  The weak convergence test has been passed twice.
c  Save the current h polynomial and shift.
c
        do i = 1, n
          shr(i) = hr(i)
          shi(i) = hi(i)
        end do
        svsr = sr
        svsi = si
c
c  Start the third stage iteration. 
c
        call vrshft ( nn, 10, are, mre, sr, si, pr, pi, qpr, qpi,
     &    pvr, pvi, hr, hi, qhr, qhi, tr, ti, zr, zi, conv )

        if ( conv ) then
          return
        end if
c
c  The iteration failed to converge. 
c  Turn off testing and restore h, s, pv and t.
c
        test = .false.

        do i = 1, n
          hr(i) = shr(i)
          hi(i) = shi(i)
        end do
        sr = svsr
        si = svsi
        call polyev ( nn, sr, si, pr, pi, qpr, qpi, pvr, pvi )
        call calct ( nn, are, bool, sr, si, hr, hi, qhr, qhi, pvr, 
     &    pvi, tr, ti )

      end do
c
c  Attempt an iteration with final h polynomial from second stage.
c
      call vrshft ( nn, 10, are, mre, sr, si, pr, pi, qpr, qpi,
     &  pvr, pvi, hr, hi, qhr, qhi, tr, ti, zr, zi, conv )

      return
      end
      subroutine vrshft ( nn, l3, are, mre, sr, si, pr, pi, qpr, qpi,
     &  pvr, pvi, hr, hi, qhr, qhi, tr, ti, zr, zi, conv )

c*********************************************************************72
c
cc vrshft() carries out the third stage iteration.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    13 March 2024
c
c  Author:
c
c    Original Fortran77 version by Michael Jenkins, Joseph Traub.
c    This version by John Burkardt.
c
c  Reference:
c
c    Michael Jenkins, Joseph Traub,
c    Algorithm 419: Zeros of a Complex Polynomial,
c    Communications of the Association for Computing Machinery,
c    Volume 15, Number 2, February 1972, pages 97-99.
c
c  Input:
c
c    integer l3: the maximum number of steps to take in stage 3.
c
c    double precision are: error bound on complex addition.
c
c    double precision mre: error bound on complex multiplication.
c
c    double precision zr, zi: the real and imaginary parts of the 
c    initial iterate.
c
c  Output:
c
c    double precision zr, zi: the real and imaginary parts of the
c    final iterate.
c
c    logical conv: true if the iteration converged.
c
      implicit none

      integer nn

      double precision are
      logical b
      double precision base
      logical bool
      double precision cmod
      logical conv
      double precision errev
      double precision eta
      double precision hi(nn)
      double precision hr(nn)
      integer i
      double precision infin
      integer j
      integer l3
      double precision mp
      double precision mre
      double precision ms
      double precision omp
      double precision pi(nn)
      double precision pr(nn)
      double precision pvi
      double precision pvr
      double precision qhi(nn)
      double precision qhr(nn)
      double precision qpi(nn)
      double precision qpr(nn)
      double precision r1
      double precision r2
      double precision relstp
      double precision si
      double precision smalno
      double precision sr
      double precision tp
      double precision ti
      double precision tr
      double precision zi
      double precision zr

      call mcon ( eta, infin, smalno, base )

      conv = .false.
      b = .false.
      sr = zr
      si = zi
c
c  Main loop for stage three
c 
      do i = 1, l3
c
c  Evaluate polynomial at s and test for convergence.
c
        call polyev ( nn, sr, si, pr, pi, qpr, qpi, pvr, pvi )
        mp = cmod ( pvr, pvi )
        ms = cmod ( sr, si )
c
c  Polynomial value is smaller in value than a bound on the error
c  in evaluating p, terminate the iteration.
c
        if ( mp < 20.0d0 * 
     &    errev ( nn, qpr, qpi, ms, mp, are, mre ) ) then
          conv = .true.
          zr = sr
          zi = si
          return
        end if

        if ( i == 1 ) then
          omp = mp
          cycle
        end if

        if ( b .or. mp < omp .or. 0.05d0 <= relstp ) then

          if ( omp < mp * 0.1d0 ) then
            return
          end if
          omp = mp
c
c  Iteration has stalled. probably a cluster of zeros. do 5 fixed
c  shift steps into the cluster to force one zero to dominate.
c
        else

          tp = relstp
          b = .true.
          if ( relstp < eta ) then
            tp = eta
          end if
          r1 = sqrt ( 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 j = 1, 5
            call calct ( nn, are, bool, sr, si, hr, hi, qhr, qhi, pvr, 
     &        pvi, tr, ti )
            call nexth ( bool, nn, qhr, qhi, qpr, qpi, tr, ti, hr, hi )
          end do

          omp = infin

        end if

        call calct ( nn, are, bool, sr, si, hr, hi, qhr, qhi, pvr, 
     &    pvi, tr, ti )
        call nexth ( bool, nn, qhr, qhi, qpr, qpi, tr, ti, hr, hi )
        call calct ( nn, are, bool, sr, si, hr, hi, qhr, qhi, pvr, 
     &    pvi, tr, ti )

        if ( .not. bool ) then
          relstp = cmod ( tr, ti ) / cmod ( sr, si )
          sr = sr + tr
          si = si + ti
        end if

      end do

      return
      end
      subroutine calct ( nn, are, bool, sr, si, hr, hi, qhr, qhi, pvr, 
     &  pvi, tr, ti )

c*********************************************************************72
c
cc calct() computes t = -p(s)/h(s).
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    14 March 2024
c
c  Author:
c
c    Original Fortran77 version by Michael Jenkins, Joseph Traub.
c    This version by John Burkardt.
c
c  Reference:
c
c    Michael Jenkins, Joseph Traub,
c    Algorithm 419: Zeros of a Complex Polynomial,
c    Communications of the Association for Computing Machinery,
c    Volume 15, Number 2, February 1972, pages 97-99.
c
c  Input:
c
c    double precision are: error bound on complex addition.
c
c  Output:
c
c    logical bool: set true if h(s) is essentially zero.
c
      implicit none

      integer nn

      double precision are
      logical bool
      double precision cmod
      double precision hi(nn)
      double precision hr(nn)
      double precision hvi
      double precision hvr
      integer n
      double precision pvi
      double precision pvr
      double precision qhi(nn)
      double precision qhr(nn)
      double precision si
      double precision sr
      double precision ti
      double precision tr

      n = nn - 1
c
c  Evaluate h(s).
c
      call polyev ( n, sr, si, hr, hi, qhr, qhi, hvr, hvi )

      bool = cmod ( hvr, hvi ) <= are * 10.0d0 * cmod ( hr(n), hi(n) )

      if ( .not. bool ) then
        call cdivid ( -pvr, -pvi, hvr, hvi, tr, ti )
      else
        tr = 0.0d0
        ti = 0.0d0
      end if

      return
      end
      subroutine nexth ( bool, nn, qhr, qhi, qpr, qpi, tr, ti, hr, hi )

c*********************************************************************72
c
cc nexth() calculates the next shifted h polynomial.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    14 March 2024
c
c  Author:
c
c    Original Fortran77 version by Michael Jenkins, Joseph Traub.
c    This version by John Burkardt.
c
c  Reference:
c
c    Michael Jenkins, Joseph Traub,
c    Algorithm 419: Zeros of a Complex Polynomial,
c    Communications of the Association for Computing Machinery,
c    Volume 15, Number 2, February 1972, pages 97-99.
c
c  Input:
c
c    logical bool: if true, h(s) is essentially zero.
c
c  Output:
c
      implicit none

      integer nn

      logical bool
      double precision hr(nn)
      double precision hi(nn)
      integer j
      integer n
      integer nm1
      double precision qhr(nn)
      double precision qhi(nn)
      double precision qpr(nn)
      double precision qpi(nn)
      double precision t1
      double precision t2
      double precision ti
      double precision tr

      n = nn - 1
      nm1 = n - 1

      if ( .not. bool ) then

        do 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)
        end do
        hr(1) = qpr(1)
        hi(1) = qpi(1)
c
c  If h(s) is zero, replace h with qh.
c
      else

        do j = 2, n
          hr(j) = qhr(j-1)
          hi(j) = qhi(j-1)
        end do
        hr(1) = 0.0d0
        hi(1) = 0.0d0

      end if

      return
      end
      subroutine polyev ( nn, sr, si, pr, pi, qr, qi, pvr, pvi )

c*********************************************************************72
c
cc polyev() evaluates a polynomial by the Horner recurrence.
c
c  Discussion:
c
c    It places the partial sums in q and the computed value in pv.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    13 March 2024
c
c  Author:
c
c    Original Fortran77 version by Michael Jenkins, Joseph Traub.
c    This version by John Burkardt.
c
c  Reference:
c
c    Michael Jenkins, Joseph Traub,
c    Algorithm 419: Zeros of a Complex Polynomial,
c    Communications of the Association for Computing Machinery,
c    Volume 15, Number 2, February 1972, pages 97-99.
c
c  Input:
c
c    integer nn: the polynomial order.
c
c    double precison sr, si: the real and imaginary parts of the
c    polynomial argument.
c
c    double precison pr(nn), pi(nn): the real and imaginary parts of the
c    polynomial coefficients.
c
c  Output:
c
c    double precision qr(nn), qi(nn): the partial sums.
c
c    double precision pvr, pvi: the real and imaginary parts of the
c    polynomial value.
c
      implicit none

      integer nn

      integer i
      double precision pi(nn)
      double precision pr(nn)
      double precision pvi
      double precision pvr
      double precision qi(nn)
      double precision qr(nn)
      double precision si
      double precision sr
      double precision t

      qr(1) = pr(1)
      qi(1) = pi(1)
      pvr = qr(1)
      pvi = qi(1)
      do 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
      end do

      return
      end
      function errev ( nn, qr, qi, ms, mp, are, mre )

c*********************************************************************72
c
cc errev() bounds the error evaluating the polynomial by the Horner recurrence.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    13 March 2024
c
c  Author:
c
c    Original Fortran77 version by Michael Jenkins, Joseph Traub.
c    This version by John Burkardt.
c
c  Reference:
c
c    Michael Jenkins, Joseph Traub,
c    Algorithm 419: Zeros of a Complex Polynomial,
c    Communications of the Association for Computing Machinery,
c    Volume 15, Number 2, February 1972, pages 97-99.
c
c  Input:
c
c    integer nn: the order of the polynomial.
c
c    double precision qr(nn), qi(nn): the partial sums.
c
c    double precision ms: the modulus of the points.
c
c    double precision mp: the modulus of the polynomial value.
c
c    double precision are, mre: error bounds on complex addition 
c    and multiplication.
c
c  Output:
c
c    double precision errev: the polynomial evaluation error bound.
c
      implicit none

      integer nn

      double precision are
      double precision cmod
      double precision e
      double precision errev
      integer i
      double precision mp
      double precision mre
      double precision ms
      double precision qi(nn)
      double precision qr(nn)

      e = cmod ( qr(1), qi(1) ) * mre / ( are + mre )
      do i = 1, nn
        e = e * ms + cmod ( qr(i), qi(i) )
      end do

      errev = e * ( are + mre ) - mp * mre

      return
      end
      function cauchy ( nn, pt )

c*********************************************************************72
c
cc cauchy() computes a lower bound on the moduli of the zeros of a polynomial.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    13 March 2024
c
c  Author:
c
c    Original Fortran77 version by Michael Jenkins, Joseph Traub.
c    This version by John Burkardt.
c
c  Reference:
c
c    Michael Jenkins, Joseph Traub,
c    Algorithm 419: Zeros of a Complex Polynomial,
c    Communications of the Association for Computing Machinery,
c    Volume 15, Number 2, February 1972, pages 97-99.
c
c  Input:
c
c    integer nn: the order of the polynomial.
c
c    double precision pt(nn): the modulus of the coefficients.
c
c  Output:
c
c    double precision cauchy: a lower bound on the moduli of the zeros.
c
      implicit none

      integer nn

      double precision cauchy
      double precision df
      double precision dx
      double precision f
      integer i
      integer n
      double precision pt(nn)
      double precision q(nn)
      double precision x
      double precision xm

      pt(nn) = - pt(nn)
c
c  Compute upper estimate of bound.
c
      n = nn - 1
      x = exp ( ( log ( - pt(nn) ) - log ( pt(1) ) ) / float ( n ) )
c
c  If Newton step at the origin is better, use it.
c
      if ( pt(n) /= 0.0d0 ) then
        xm = - pt(nn) / pt(n)
        x = min ( x, xm )
      end if
c
c  Chop the interval (0,x) unitl f <= 0.
c
      do

        xm = x * 0.1d0

        f = pt(1)
        do i = 2, nn
          f = f * xm + pt(i)
        end do

        if ( f <= 0.0d0 ) then
          exit
        end if

        x = xm
 
      end do

      dx = x
c
c  Do Newton iteration until x converges to two decimal places.
c
      do while ( 0.005D0 < abs ( dx / x ) )

        q(1) = pt(1)
        do i = 2, nn
          q(i) = q(i-1) * x + pt(i)
        end do
        f = q(nn)

        df = q(1)
        do i = 2, n
          df = df * x + q(i)
        end do

        dx = f / df
        x = x - dx

      end do

      cauchy = x

      return
      end
      function rescale ( nn, pt )

c*********************************************************************72
c
cc rescale() returns a scale factor for the coefficients of the polynomial. 
c
c  Discussion:
c
c    The scaling is done to avoid overflow and to avoid undetected underflow 
c    interfering with the convergence criterion.  
c
c    The factor is a power of the base.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    13 March 2024
c
c  Author:
c
c    Original Fortran77 version by Michael Jenkins, Joseph Traub.
c    This version by John Burkardt.
c
c  Reference:
c
c    Michael Jenkins, Joseph Traub,
c    Algorithm 419: Zeros of a Complex Polynomial,
c    Communications of the Association for Computing Machinery,
c    Volume 15, Number 2, February 1972, pages 97-99.
c
c  Input:
c
c    integer nn: the order of the polynomial.
c
c    double precision pt(nn): modulus of coefficients of the polynomial.
c
c  Output:
c
c    double precision rescale: the scale factor.
c
      implicit none

      integer nn

      double precision base
      double precision eta
      double precision hi
      integer i
      double precision infin
      integer l
      double precision lo
      double precision pt(nn)
      double precision rescale
      double precision sc
      double precision smalno
      double precision x
      double precision xmax
      double precision xmin
c
c  Get machine constants.
c
      call mcon ( eta, infin, smalno, base )
c
c  Find largest and smallest moduli of coefficients.
c
      hi = sqrt ( infin )
      lo = smalno / eta
      xmax = 0.0d0
      xmin = infin

      do i = 1, nn
        x = pt(i)
        xmax = max ( xmax, x )
        if ( x /= 0.0d0 ) then
          xmin = min ( xmin, x )
        end if
      end do
c
c  Scale only if there are very large or very small components.
c
      rescale = 1.0d0
      if ( lo <= xmin .and. xmax <= hi ) then
        return
      end if

      x = lo / xmin
      if ( x <= 1.0d0 ) then
        sc = 1.0d0 / ( sqrt ( xmax ) * sqrt ( xmin ) )
      else
        sc = x
        if ( xmax < infin / sc ) then
          sc = 1.0d0
        end if
      end if

      l = int ( log ( sc ) / log ( base ) + 0.5d0 )
      rescale = base ** l

      return
      end
      subroutine cdivid ( ar, ai, br, bi, cr, ci )

c*********************************************************************72
c
cc cdivid() performs complex division c = a/b, avoiding overflow.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    13 March 2024
c
c  Author:
c
c    Original Fortran77 version by Michael Jenkins, Joseph Traub.
c    This version by John Burkardt.
c
c  Reference:
c
c    Michael Jenkins, Joseph Traub,
c    Algorithm 419: Zeros of a Complex Polynomial,
c    Communications of the Association for Computing Machinery,
c    Volume 15, Number 2, February 1972, pages 97-99.
c
c  Input:
c
c    double precision ar, ai: the real and imaginary part of the numerator.
c
c    double precision br, bi: the real and imaginary part of the denominator.
c
c  Output:
c
c    double precision cr, ci: the real and imaginary part of the result.
c
      implicit none

      double precision ar
      double precision ai
      double precision base
      double precision br
      double precision bi
      double precision cr
      double precision ci
      double precision d
      double precision eta
      double precision infin
      double precision r
      double precision smalno
c
c  Division by zero, c = infinity.
c
      if ( br == 0.0d0 .and. bi == 0.0d0 ) then

        call mcon ( eta, infin, smalno, base )

        cr = infin
        ci = infin

      else if ( abs ( br ) < abs ( bi ) ) then

        r = br / bi
        d = bi + r * br
        cr = ( ar * r + ai ) / d
        ci = ( ai * r - ar ) / d

      else

        r = bi / br
        d = br + r * bi
        cr = ( ar + ai * r ) / d
        ci = ( ai - ar * r ) / d

      end if

      return
      end
      function cmod ( r, i )

c*********************************************************************72
c
cc cmod() computes the modulus of a complex number, avoiding overflow.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    13 March 2024
c
c  Author:
c
c    Original Fortran77 version by Michael Jenkins, Joseph Traub.
c    This version by John Burkardt.
c
c  Reference:
c
c    Michael Jenkins, Joseph Traub,
c    Algorithm 419: Zeros of a Complex Polynomial,
c    Communications of the Association for Computing Machinery,
c    Volume 15, Number 2, February 1972, pages 97-99.
c
c  Input:
c
c    double precision r, i: the real and imaginary parts of a complex number.
c
c  Output:
c
c    double precision cmod: the modulus of the complex number.
c
      implicit none

      double precision ai
      double precision ar
      double precision cmod
      double precision i
      double precision r

      ar = abs ( r )
      ai = abs ( i )

      if ( ar < ai ) then
        cmod = ai * sqrt ( 1.0d0 + ( ar / ai )**2 )
      else if ( ai < ar ) then
        cmod = ar * sqrt ( 1.0d0 + ( ai / ar )**2 )
      else
        cmod = ar * sqrt ( 2.0d0 )
      end if

      return
      end
      subroutine mcon ( eta, infin, smalno, base )

c*********************************************************************72
c
cc mcon() provides machine constants.
c
c  Discussion:
c
c    The user may either set them directly or use the
c    statements below to compute them. 
c
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
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
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    13 March 2024
c
c  Author:
c
c    Original Fortran77 version by Michael Jenkins, Joseph Traub.
c    This version by John Burkardt.
c
c  Reference:
c
c    Michael Jenkins, Joseph Traub,
c    Algorithm 419: Zeros of a Complex Polynomial,
c    Communications of the Association for Computing Machinery,
c    Volume 15, Number 2, February 1972, pages 97-99.
c
c  Output:
c
c    double precision eta: the smallest positive floating-point number such 
c    that 1.0d0 + eta is greater than 1.0d0.
c
c    double precision infin: the largest floating-point number.
c
c    double precision smalno: the smallest positive floating-point number.
c
c    double precision base: the base of the floating-point number system used.
c
      implicit none

      double precision base
      double precision eta
      double precision infin
      integer m
      integer n
      double precision smalno
      integer t

      base = 16.0d0
      t = 14
      m = 63
      n = -65
      eta = base ** ( 1 - t )
      infin = base * ( 1.0d0 - base ** ( - t ) ) * base ** ( m - 1 )
      smalno = ( base ** ( n + 3 ) ) / base ** 3

      return
      end