subroutine rpoly ( op, degree, zeror, zeroi, fail )

c*********************************************************************72
c
cc rpoly() finds the zeros of a real polynomial.
c
c  Discussion:
c
c    In the original code, the input quantity DEGREE could be altered
c    before output.  This is an unrighteous idea and has been quashed.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    21 March 2024
c
c  Author:
c
c    Original Fortran77 version by Michael Jenkins.
c    This version by John Burkardt.
c
c  Reference:
c
c    Michael Jenkins,
c    Algorithm 493: Zeros of a Real Polynomial,
c    ACM Transactions on Mathematical Software,
c    Volume 1, Number 2, June 1975, pages 178-189.
c
c  Input:
c
c    double precision op(degree+1): 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): real and imaginary parts 
c    of the zeros.
c
c    logical fail: True only if leading coefficient is zero or if RPOLY
c    has found fewer than DEGREE zeros.  
c
      implicit none

      integer degree

      double precision a
      double precision a1
      double precision a3
      double precision a7
      double precision aa
      double precision are
      double precision b
      double precision base
      double precision bb
      double precision bnd
      double precision c
      double precision cc
      integer cnt
      double precision cosr
      double precision d
      double precision df
      double precision dx
      double precision e
      double precision eta
      double precision f
      double precision factor
      logical fail
      double precision ff
      double precision g
      double precision h
      integer i
      double precision infin
      integer j
      integer jj
      double precision k(degree+1)
      integer l
      integer l2
      double precision lo
      double precision lzi
      double precision lzr
      double precision mre
      integer n
      integer nn
      integer nz
      double precision op(degree+1)
      double precision p(degree+1)
      double precision pt(degree+1)
      double precision qk(degree+1)
      double precision qp(degree+1)
      double precision smalno
      double precision svk(degree+1)
      double precision t
      double precision temp(degree+1)
      double precision sc
      double precision si
      double precision sinr
      double precision sr
      double precision szi
      double precision szr
      double precision u
      double precision v
      double precision x
      double precision xm
      double precision xmax
      double precision xmin
      double precision xx
      double precision xxx
      double precision yy
      double precision zeroi(degree)
      logical zerok
      double precision zeror(degree)
c
c THE FOLLOWING STATEMENTS SET MACHINE CONSTANTS USED
c IN VARIOUS PARTS OF THE PROGRAM. THE MEANING OF THE
c FOUR CONSTANTS ARE...
c ETA     THE MAXIMUM RELATIVE REPRESENTATION ERROR
c         WHICH CAN BE DESCRIBED AS THE SMALLEST
c         POSITIVE FLOATING POINT NUMBER SUCH THAT
c         1.0+ETA IS GREATER THAN 1.
c INFINY  THE LARGEST FLOATING-POINT NUMBER.
c SMALNO  THE SMALLEST POSITIVE FLOATING-POINT NUMBER
c         IF THE EXPONENT RANGE DIFFERS IN SINGLE AND
c         DOUBLE PRECISION THEN SMALNO AND INFIN
c         SHOULD INDICATE THE SMALLER RANGE.
c
c BASE    THE BASE OF THE FLOATING-POINT NUMBER SYSTEM USED.
c
      base = 2.0
      eta = epsilon ( eta )
      infin = huge ( infin )
      smalno = tiny ( smalno )
c
c  ARE and MRE are the unit error in addition and multiplication.
c
      are = eta
      mre = eta
      lo = smalno / eta
c
c  Initialization of constants for shift rotation.
c
      xx = 0.70710678
      yy = - xx
      cosr = -0.069756474
      sinr = 0.99756405
      fail = .false.
      n = degree
      nn = n + 1
c
c  Algorithm fails if the leading coefficient is zero.
c
      if ( op(1) == 0.0 ) then
        fail = .true.
        return
      end if
c
c  Remove the zeros at the origin if any.
c
      do while ( op(nn) == 0.0 )
        j = degree - n + 1
        zeror(j) = 0.0
        zeroi(j) = 0.0
        nn = nn - 1
        n = n - 1
      end do
c
c  Copy the coefficients.
c
      do i = 1, nn
        p(i) = op(i)
      end do
c
c  Start the algorithm for one zero.
c
   40 continue

      if ( n < 1 ) then
        return
      else if ( n == 1 ) then
        zeror(degree) = - p(2) / p(1)
        zeroi(degree) = 0.0
        return
      else if ( n == 2 ) then
        call quad ( p(1), p(2), p(3), zeror(degree-1), 
     &    zeroi(degree-1), zeror(degree), zeroi(degree) )
        return
      end if
c
c  Find largest and smallest moduli of coefficients.
c
      xmax = 0.0
      xmin = infin
      do i = 1, nn
        x = abs ( p(i) )
        xmax = max ( xmax, x )
        if ( x /= 0.0 ) then
          xmin = min ( xmin, x )
        end if
      end do
c
c  Scale if there are large or very small coefficients.
c
      sc = lo / xmin

      if ( 1.0 < sc ) then
        if ( xmax <= infin ) then
          l = int ( log ( sc ) / log ( base ) + 0.5 )
          factor = ( base * 1.0 ) ** l

          if ( factor /= 1.0 ) then
            do i = 1, nn
              p(i) = factor * p(i)
            end do
          end if
        end if
        go to 110
      end if

      if ( xmax < 10.0 ) then
        go to 110
      end if

      if ( sc == 0.0 ) then
        sc = smalno
      end if

      l = int ( log ( sc ) / log ( base ) + 0.5 )
      factor = ( base * 1.0 ) ** l

      if ( factor /= 1.0 ) then
        do i = 1, nn
          p(i) = factor * p(i)
        end do
      end if
c
c  Compute lower bound on moduli of zeros.
c
110   continue

      do i = 1, nn
        pt(i) = abs ( p(i) )
      end do
      pt(nn) = - pt(nn)
c
c  Compute upper estimate of bound.
c
      x = exp (( log ( -pt(nn) )- log ( pt(1) ) ) / dble ( n ) )
c
c  If Newton step at the origin is better, use it.
c
      if ( pt(n) /= 0.0 ) then
        xm = - pt(nn) / pt(n)
        x = min ( x, xm )
      end if
c
c  Chop the interval (0,X) until FF <= 0
c
      do

        xm = x * 0.1D+00
        ff = pt(1)
        do i = 2, nn
          ff = ff * xm + pt(i)
        end do

        if ( ff <= 0.0 ) 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.005 < abs ( dx / x ) )
        ff = pt(1)
        df = ff
        do i = 2, n
          ff = ff * x + pt(i)
          df = df * x + ff
        end do
        ff = ff * x + pt(nn)
        dx = ff / df
        x = x - dx
      end do

      bnd = x
c
c  Compute the derivative as the initial K polynomial.
c
      do i = 2, n
        k(i) = dble ( nn - i ) * p(i) / dble ( n )
      end do

      k(1) = p(1)
      aa = p(nn)
      bb = p(n)
      zerok = k(n) == 0.0
c
c  Do 5 steps with no shift.
c
      do jj = 1, 5

        cc = k(n)
c
c  Use scaled form of recurrence if value of k at 0 is nonzero.
c
        if ( .not. zerok ) then

          t = - aa / cc
          do i = 1, n - 1
            j = nn - i
            k(j) = t * k(j-1) + p(j)
          end do

          k(1) = p(1)
          zerok = abs ( k(n) ) <= abs ( bb ) * eta * 10.0
c
c  Use unscaled form of recurrence.
c
        else

          do i = 1, n - 1
            j = nn - i
            k(j) = k(j-1)
          end do

          k(1) = 0.0
          zerok = k(n) == 0.0

        end if

      end do
c
c  Save K for restarts with new shifts.
c
      do i = 1, n
        temp(i) = k(i)
      end do
c
c  Loop to select the quadratic corresponding to each new shift.
c
      do cnt = 1, 20
c
c  Quadratic corresponds to a double shift to a
c  non-real point and its complex conjugate. 
c  The point has modulus BND and amplitude rotated by 94 degrees
c  from the previous shift.
c
        xxx = cosr * xx - sinr * yy
        yy =  sinr * xx + cosr * yy
        xx = xxx
        sr = bnd * xx
        si = bnd * yy
        u = -2.0 * sr
        v = bnd
c
c  Second stage calculation, fixed quadratic.
c
        l2 = 20 * cnt
        call fxshfr ( a, a1, a3, a7, are, b, c, d, e, eta, 
     &    f, g, h, k, l2, lzr, lzi, mre, n, nn, p, qk, qp, sr, svk, 
     &    szr, szi, u, v, nz )
c
c  The second stage jumps directly to one of the third stage iterations 
c  and returns here if successful.
c
c  Deflate the polynomial, store the zero or zeros and
c  return to the main algorithm.
c
        if ( 0 < nz ) then

          j = degree - n + 1
          zeror(j) = szr
          zeroi(j) = szi
          nn = nn - nz
          n = nn - 1
          do i = 1, nn
            p(i) = qp(i)
          end do

          if ( 1 < nz ) then
            zeror(j+1) = lzr
            zeroi(j+1) = lzi
          end if

          go to 40

        end if
c
c  If the iteration is unsuccessful another quadratic
c  is chosen after restoring K.
c
        do i = 1, n
          k(i) = temp(i)
        end do

      end do
c
c  Return with failure if no convergence with 20 shifts.
c
      fail = .true.

      return
      end
      subroutine fxshfr ( a, a1, a3, a7, are, b, c, d, e, eta, 
     &  f, g, h, k, l2, lzr, lzi, mre, n, nn, p, qk, qp, sr, svk, 
     &  szr, szi, u, v, nz )

c*********************************************************************72
c
cc fxshfr() computes up to L2 fixed shift K-polynomials.
c
c  Discussion:
c
c    In the linear or quadratic cases, a convergence test is made.
c
c    The code initiates one of the variable shift iterations and 
c    returns the number of zeros found.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    21 March 2024
c
c  Author:
c
c    Original Fortran77 version by Michael Jenkins.
c    This version by John Burkardt.
c
c  Reference:
c
c    Michael Jenkins,
c    Algorithm 493: Zeros of a Real Polynomial,
c    ACM Transactions on Mathematical Software,
c    Volume 1, Number 2, June 1975, pages 178-189.
c
c  Input:
c
c    integer L2: limit of fixed shift steps
c
c  Output:
c
c    integer NZ: the number of zeros found.
c
      implicit none

      integer nn

      double precision a
      double precision a1
      double precision a3
      double precision a7
      double precision are
      double precision b
      double precision betas
      double precision betav
      double precision c
      double precision d
      logical doit
      double precision e
      double precision eta
      double precision f
      double precision g
      double precision h
      integer i
      integer iflag
      integer j
      double precision k(nn)
      integer l2
      double precision lzi
      double precision lzr
      double precision mre
      integer n
      integer nz
      double precision oss
      double precision ots
      double precision otv
      double precision ovv
      double precision p(nn)
      double precision qk(nn)
      double precision qp(nn)
      double precision s
      logical spass
      double precision sr
      double precision ss
      logical stry
      double precision svk(nn)
      double precision svu
      double precision svv
      double precision szi
      double precision szr
      double precision ts
      double precision tss
      double precision tv
      double precision tvv
      integer type
      double precision u
      double precision ui
      double precision v
      double precision vi
      logical vpass
      logical vtry
      double precision vv

      nz = 0
      betav = 0.25
      betas = 0.25
      oss = sr
      ovv = v
c
c  Evaluate polynomial by synthetic division.
c
      call quadsd ( nn, u, v, p, qp, a, b )

      call calcsc ( nn, n, a, b, eta, k, type, u, v, a1, a3, a7, 
     &  c, d, e, f, g, h, qk )

      do j = 1, l2
c
c  Calculate next K polynomial and estimate V.
c
        call nextk ( nn, type, a, a1, a3, a7, b, eta, qk, qp, k )

        call calcsc ( nn, n, a, b, eta, k, type, u, v, a1, a3, a7, 
     &    c, d, e, f, g, h, qk )

        call newest ( nn, n, a, a1, a3, a7, b, c, d, f, g, h, 
     &    k, p, type, u, v, ui, vi )
        vv = vi
c
c  Estimate S.
c
        if ( k(n) == 0.0 ) then
          ss = 0.0
        else
          ss = - p(nn) / k(n)
        end if

        tv = 1.0
        ts = 1.0

        if ( j == 1 .or. type == 3 ) then
          ovv = vv
          oss = ss
          otv = tv
          ots = ts
          cycle
        end if
c
c  Compute relative measures of convergence of S and V sequences.
c
        if ( vv /= 0.0 ) then
          tv = abs ( ( vv - ovv ) / vv )
        end if

        if ( ss /= 0.0 ) then
          ts = abs ( ( ss - oss ) / ss )
        end if
c
c  If decreasing, multiply two most recent convergence measures.
c
        if ( tv < otv ) then
          tvv = tv * otv
        else
          tvv = 1.0
        end if

        if ( ts < ots ) then
          tss = ts * ots
        else
          tss = 1.0
        end if
c
c  Compare with convergence criteria.
c
        vpass = tvv < betav
        spass = tss < betas

        if ( .not. ( spass .or. vpass ) ) then
          ovv = vv
          oss = ss
          otv = tv
          ots = ts
          cycle
        end if
c
c  At least one sequence has passed the convergence test. 
c  Store variables before iterating
c
        svu = u
        svv = v
        do i = 1, n
          svk(i) = k(i)
        end do
        s = ss
c
c  Choose iteration according to the fastest converging sequence.
c
        vtry = .false.
        stry = .false.

        if ( spass .and. ( ( .not. vpass ) .or. tss < tvv ) ) then
          doit = .false.
        else
          doit = .true.
        end if

       do while ( .true. )

          if ( .not. doit ) then

            doit = .true.

          else

            call quadit ( a, a1, a3, a7, are, b, c, d, e, eta, f, g, 
     &        h, k, lzr, lzi, mre, n, nn, p, qk, qp, szr, szi, 
     &        u, v, ui, vi, nz )

            if ( 0 < nz ) then
              return
            end if
c
c  Quadratic iteration has failed. Flag that it has
c  been tried and decrease the convergence criterion.
c
            vtry = .true.
            betav = betav * 0.25
c
c  Try linear iteration if it has not been tried and
c  the S sequence is converging.
c
            if ( stry .or. ( .not. spass ) ) then

              u = svu
              v = svv
              do i = 1, n
                k(i) = svk(i)
              end do
c
c  Try quadratic iteration if it has not been tried
c  and the v sequence is converging.
c
              if ( .not. vpass .or. vtry ) then
                exit
              end if

              cycle

            end if

            do i = 1, n
              k(i) = svk(i)
            end do

          end if

          call realit ( are, eta, k, mre, n, nn, p, qk, qp, szr, szi, 
     &      s, nz, iflag )

          if ( 0 < nz ) then
            return
          end if
c
c  Linear iteration has failed. Flag that it has been
c  tried and decrease the convergence criterion.
c
          stry = .true.
          betas = betas * 0.25
c
c  If linear iteration signals an almost double real
c  zero, attempt quadratic interation.
c
          if ( iflag /= 0 ) then
            ui = - ( s + s )
            vi = s * s
            cycle
          end if
c
c  Restore variables.
c
          u = svu
          v = svv
          do i = 1, n
            k(i) = svk(i)
          end do
c
c  Try quadratic iteration if it has not been tried
c  and the v sequence is converging.
c
          if ( .not. vpass .or. vtry ) then
            exit
          end if

        end do
c
c  Recompute QP and scalar values to continue the second stage.
c
        call quadsd ( nn, u, v, p, qp, a, b )

        call calcsc ( nn, n, a, b, eta, k, type, u, v, a1, a3, a7, 
     &    c, d, e, f, g, h, qk )

        ovv = vv
        oss = ss
        otv = tv
        ots = ts

      end do

      return
      end
      subroutine quadit ( a, a1, a3, a7, are, b, c, d, e, eta, f, g, 
     &  h, k, lzr, lzi, mre, n, nn, p, qk, qp, szr, szi, u, v, uu, 
     &  vv, nz )

c*********************************************************************72
c
cc quadit() applies an iteration for a quadratic factor.
c
c  Discussion:
c
c    Variable-shift K-polynomial iteration for a quadratic factor.
c
c    Converges only if the zeros are equimodular or nearly so.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    20 March 2024
c
c  Author:
c
c    Original Fortran77 version by Michael Jenkins.
c    This version by John Burkardt.
c
c  Reference:
c
c    Michael Jenkins,
c    Algorithm 493: Zeros of a Real Polynomial,
c    ACM Transactions on Mathematical Software,
c    Volume 1, Number 2, June 1975, pages 178-189.
c
c  Input:
c
c    double precision UU, VV: coefficients of starting quadratic.
c
c  Output:
c
c    integer nz: the number of zeros found.
c
      implicit none

      integer nn

      double precision a
      double precision a1
      double precision a3
      double precision a7
      double precision are
      double precision b
      double precision c
      double precision d
      double precision e
      double precision ee
      double precision eta
      double precision f
      double precision g
      double precision h
      integer i
      integer j
      double precision k(nn)
      double precision lzi
      double precision lzr
      double precision mp
      double precision mre
      integer n
      integer nz
      double precision omp
      double precision p(nn)
      double precision qk(nn)
      double precision qp(nn)
      double precision relstp
      double precision szi
      double precision szr
      double precision t
      logical tried
      integer type
      double precision u
      double precision ui
      double precision uu
      double precision v
      double precision vi
      double precision vv
      double precision zm

      nz = 0
      tried = .false.
      u = uu
      v = vv
      j = 0

      do

        a = 1.0
        call quad ( a, u, v, szr, szi, lzr, lzi )
c
c  Return if roots of the quadratic are real and not
c  close to multiple or nearly equal and of opposite sign.
c
        if ( 0.01 * abs ( lzr ) 
     &    < abs ( abs ( szr ) - abs ( lzr ) ) ) then
          exit
        end if
c
c  Evaluate polynomial by quadratic synthetic division.
c
        call quadsd ( nn, u, v, p, qp, a, b )
        mp = abs ( a - szr * b ) + abs ( szi * b )
c
c  Compute a rigorous bound on the rounding error in evaluating P.
c
        zm = sqrt ( abs ( v ) )
        ee = 2.0 * abs ( qp(1) )
        t = - szr * b
        do i = 2, n
          ee = ee * zm + abs ( qp(i) )
        end do
        ee = ee * zm + abs ( a + t )
        ee = ( 5.0 * mre + 4.0 * are ) * ee 
     &    - ( 5.0 * mre + 2.0 * are ) *
     &    ( abs ( a + t ) + abs ( b ) * zm ) + 2.0 * are * abs ( t )
c
c  Iteration has converged sufficiently if the
c  polynomial value is less than 20 times this bound.
c
        if ( mp <= 20.0 * ee ) then
          nz = 2
          exit
        end if

        j = j + 1
c
c  Stop iteration after 20 steps
c
        if ( 20 < j ) then
          exit
        end if
c
c  A cluster appears to be stalling the convergence.
c  Five fixed shift steps are taken with a U, V close to the cluster.
c
        if ( 2 <= j .and. 
     &    relstp <= 0.01 .and. 
     &    omp <= mp .and. 
     &    .not. tried ) then

          relstp = max ( relstp, eta )
          relstp = sqrt ( relstp )
          u = u - u * relstp
          v = v + v * relstp

          call quadsd ( nn, u, v, p, qp, a, b )

          do i = 1, 5

            call calcsc ( nn, n, a, b, eta, k, type, u, v, a1, a3, a7, 
     &        c, d, e, f, g, h, qk )

            call nextk ( nn, type, a, a1, a3, a7, b, eta, qk, qp, k )

          end do

          tried = .true.
          j = 0

        end if

        omp = mp
c
c  Calculate next K polynomial and new U and V.
c
        call calcsc ( nn, n, a, b, eta, k, type, u, v, a1, a3, a7, 
     &  c, d, e, f, g, h, qk )

        call nextk ( nn, type, a, a1, a3, a7, b, eta, qk, qp, k )

        call calcsc ( nn, n, a, b, eta, k, type, u, v, a1, a3, a7, 
     &    c, d, e, f, g, h, qk )

        call newest ( nn, n, a, a1, a3, a7, b, c, d, f, g, h, 
     &    k, p, type, u, v, ui, vi )
c
c  If VI is zero, the iteration is not converging.
c
        if ( vi == 0.0 ) then
          exit
        end if

        relstp = abs ( ( vi - v ) / vi )
        u = ui
        v = vi
      
      end do

      return
      end
      subroutine realit ( are, eta, k, mre, n, nn, p, qk, qp, szr, szi, 
     &  sss, nz, iflag )

c*********************************************************************72
c
cc realit() is a variable-shift H polynomial iteration for a real zero.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    18 March 2024
c
c  Author:
c
c    Original Fortran77 version by Michael Jenkins.
c    This version by John Burkardt.
c
c  Reference:
c
c    Michael Jenkins,
c    Algorithm 493: Zeros of a Real Polynomial,
c    ACM Transactions on Mathematical Software,
c    Volume 1, Number 2, June 1975, pages 178-189.
c
c  Input:
c
c    double precision sss: the starting iterate.
c
c  Output:
c
c    double precision sss: the final iterate.
c
c    integer nz: the number of zeros found.
c
c    integer iflag: indicates a pair of zeros near real axis.
c
      implicit none

      integer n
      integer nn

      double precision are
      double precision ee
      double precision eta
      integer i
      integer iflag
      integer j
      double precision k(nn)
      double precision kv
      double precision mp
      double precision mre
      double precision ms
      integer nz
      double precision omp
      double precision p(nn)
      double precision pv
      double precision qk(nn)
      double precision qp(nn)
      double precision s
      double precision sss
      double precision szi
      double precision szr
      double precision t

      nz = 0
      s = sss
      iflag = 0
      j = 0

      do

        pv = p(1)
c
c  Evaluate P(S).
c
        qp(1) = pv
        do i = 2, nn
          pv = pv * s + p(i)
          qp(i) = pv
        end do
        mp = abs ( pv )
c
c  Compute a rigorous bound on the error in evaluating P.
c
        ms = abs ( s )
        ee = ( mre / ( are + mre ) ) * abs ( qp(1) )
        do i = 2, nn
          ee = ee * ms + abs ( qp(i) )
        end do
c
c  Iteration has converged sufficiently if the
c  polynomial value is less than 20 times this bound.
c
        if ( mp <= 20.0 * ( ( are + mre ) * ee - mre * mp ) ) then
          nz = 1
          szr = s
          szi = 0.0
          exit
        end if

        j = j + 1
c
c  Stop iteration after 10 steps.
c
        if ( 10 < j ) then
          exit
        end if
c
c  A cluster of zeros near the real axis has been encountered.
c  Return with IFLAG set to initiate a quadratic iteration.
c
        if ( 2 <= j .and.
     &    abs ( t ) < 0.001 * abs ( s - t ) .and.
     &    omp < mp ) then
          iflag = 1
          sss = s
          exit
        end if
c
c  Return if the polynomial value has increased significantly.
c
        omp = mp
c
c  Compute T, the next polynomial, and the new iterate.
c
        kv = k(1)
        qk(1) = kv
        do i = 2, n
          kv = kv * s + k(i)
          qk(i) = kv
        end do
c
c  Use the scaled form of the recurrence if the value
c  of K at S is nonzero.
c
        if ( abs ( k(n) ) * 10.0 * eta < abs ( kv ) ) then

          t = - pv / kv
          k(1) = qp(1)
          do i = 2, n
            k(i) = t * qk(i-1) + qp(i)
          end do

        else

          k(1) = 0.0
          do i = 2, n
            k(i) = qk(i-1)
          end do

        end if

        kv = k(1)
        do i = 2, n
          kv = kv * s + k(i)
        end do

        if ( abs ( kv ) <= abs ( k(n) ) * 10.0 * eta ) then
          t = 0.0
        else
          t = - pv / kv
        end if

        s = s + t

      end do

      return
      end
      subroutine calcsc ( nn, n, a, b, eta, k, type, u, v, a1, a3, a7, 
     &  c, d, e, f, g, h, qk )

c*********************************************************************72
c
cc calcsc() calculates scalar quantities for the computation.
c
c  Discussion:
c
c    The quantities are used to compute the next K polynomial and 
c    new estimates of the quadratic coefficients.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    16 March 2024
c
c  Author:
c
c    Original Fortran77 version by Michael Jenkins.
c    This version by John Burkardt.
c
c  Reference:
c
c    Michael Jenkins,
c    Algorithm 493: Zeros of a Real Polynomial,
c    ACM Transactions on Mathematical Software,
c    Volume 1, Number 2, June 1975, pages 178-189.
c
c  Input:
c
c    integer TYPE: indicates how calculations are normalized 
c    to avoid overflow.
c
c    double precision k(nn):
c
c  Output:
c
      implicit none

      integer nn

      double precision a
      double precision a1
      double precision a3
      double precision a7
      double precision b
      double precision c
      double precision d
      double precision e
      double precision eta
      double precision f
      double precision g
      double precision h
      double precision k(nn)
      integer n
      double precision qk(nn)
      integer type
      double precision u
      double precision v
c
c  Synthetic division of K by the quadratic 1,U,V.
c
      call quadsd ( n, u, v, k, qk, c, d )

      if ( abs ( c ) <= abs ( k(n) ) * 100.0 * eta .and.
     &     abs ( d ) <= abs ( k(n-1) ) * 100.0 * eta ) then
        type = 3
c
c  TYPE=3 indicates the quadratic is almost a factor of K.
c
      else if ( abs ( c ) <= abs ( d ) ) then
        type = 2
c
c  TYPE=2 indicates that all formulas are divided by D.
c
        e = a / d
        f = c / d
        g = u * b
        h = v * b
        a3 = ( a + g ) * e + h * ( b / d )
        a1 = b * f - a
        a7 = ( f + u ) * a + h
c
c  TYPE=1 indicates that all formulas are divided by C
c
      else

        type = 1
        e = a / c
        f = d / c
        g = u * e
        h = v * b
        a3 = a * e + ( h / c + g ) * b
        a1 = b - a * ( d / c )
        a7 = a + g * d + h * f

      end if

      return
      end
      subroutine nextk ( nn, type, a, a1, a3, a7, b, eta, qk, qp, k )

c*********************************************************************72
c
cc nextk() computes the next K polynomials using scalars from calcsc().
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    16 March 2024
c
c  Author:
c
c    Original Fortran77 version by Michael Jenkins.
c    This version by John Burkardt.
c
c  Reference:
c
c    Michael Jenkins,
c    Algorithm 493: Zeros of a Real Polynomial,
c    ACM Transactions on Mathematical Software,
c    Volume 1, Number 2, June 1975, pages 178-189.
c
c  Input:
c
c    integer nn:
c
c    integer type:
c
c    double precision a:
c
c    double precision a1, a3, a7:
c
c    double precision b:
c
c    double precision eta:
c
c    double precision qk(nn):
c
c    double precision qp(nn):
c
c  Output:
c
c    double precision k(nn):
c
      implicit none

      integer nn

      double precision a
      double precision a1
      double precision a3
      double precision a7
      double precision b
      double precision eta
      integer i
      double precision k(nn)
      double precision qk(nn)
      double precision qp(nn)
      double precision temp
      integer type

      if ( type == 1 ) then
        temp = b
      else if ( type == 2 ) then
        temp = a
      end if

      if ( type == 3 ) then

        k(1) = 0.0
        k(2) = 0.0
        do i = 3, nn
          k(i) = qk(i-2)
        end do

      else if ( abs ( temp ) * eta * 10.0 <= abs ( a1 ) ) then

        a7 = a7 / a1
        a3 = a3 / a1
        k(1) = qp(1)
        k(2) = qp(2) - a7 * qp(1)
        do i = 3, nn - 1
          k(i) = a3 * qk(i-2) - a7 * qp(i-1) + qp(i)
        end do

      else

        k(1) = 0.0
        k(2) = - a7 * qp(1)
        do i = 3, nn - 1
          k(i) = a3 * qk(i-2) - a7 * qp(i-1)
        end do

      end if

      return
      end
      subroutine newest ( nn, n, a, a1, a3, a7, b, c, d, f, g, h, 
     &  k, p, type, u, v, ui, vi )

c*********************************************************************72
c
cc newest() computes new estimates of the quadratic coefficients.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    16 March 2024
c
c  Author:
c
c    Original Fortran77 version by Michael Jenkins.
c    This version by John Burkardt.
c
c  Reference:
c
c    Michael Jenkins,
c    Algorithm 493: Zeros of a Real Polynomial,
c    ACM Transactions on Mathematical Software,
c    Volume 1, Number 2, June 1975, pages 178-189.
c
c  Input:
c
c    integer type:
c
c  Output:
c
c    double precision ui, vi:
c
      implicit none

      double precision a
      double precision a1
      double precision a3
      double precision a4
      double precision a5
      double precision a7
      double precision b
      double precision b1
      double precision b2
      double precision c
      double precision c1
      double precision c2
      double precision c3
      double precision c4
      double precision d
      double precision f
      double precision g
      double precision h
      double precision k(nn)
      integer n
      integer nn
      double precision p(nn)
      double precision temp
      integer type
      double precision u
      double precision ui
      double precision v
      double precision vi
c
c  Use formulas appropriate to setting of TYPE.
c
      if ( type == 3 ) then
        ui = 0.d0
        vi = 0.d0
        return
      end if

      if ( type == 1 ) then
        a4 = a + u * b + h * f
        a5 = c + ( u + v * f ) * d
      else
        a4 = ( a + g ) * f + h
        a5 = ( f + u ) * c + v * d
      end if
c
c  Evaluate new quadratic coefficients.
c
      b1 = - k(n) / p(nn)
      b2 = - ( k(n-1) + b1 * p(n) ) / p(nn)
      c1 = v * b2 * a1
      c2 = b1 * a7
      c3 = b1 * b1 * a3
      c4 = c1 - c2 - c3
      temp = a5 + b1 * a4 - c4
      if ( temp == 0.0 ) then
        ui = 0.0
        vi = 0.0
      else
        ui = u - ( u * ( c3 + c2 ) 
     &    + v * ( b1 * a1 + b2 * a7 ) ) / temp
        vi = v * ( 1.0 + c4 / temp )
      end if

      return
      end
      subroutine quadsd ( nn, u, v, p, q, a, b )

c*********************************************************************72
c
cc quadsd() divides p by the quadratic 1,u,v using synthetic division.
c
c  Discussion:
c
c    The quotient is in Q and the remainder in A, B.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    16 March 2024
c
c  Author:
c
c    Original Fortran77 version by Michael Jenkins.
c    This version by John Burkardt.
c
c  Reference:
c
c    Michael Jenkins,
c    Algorithm 493: Zeros of a Real Polynomial,
c    ACM Transactions on Mathematical Software,
c    Volume 1, Number 2, June 1975, pages 178-189.
c
c  Input:
c
c    integer nn:
c
c    double precision u, v:
c
c    double precision p(nn): the polynomial coefficients.
c
c  Output:
c
c    double precision q(nn): the quotient.
c
c    double precision a, b: the remainder.
c
      implicit none

      integer nn

      double precision a
      double precision b
      double precision c
      integer i
      double precision p(nn)
      double precision q(nn)
      double precision u
      double precision v

      b = p(1)
      q(1) = b
      a = p(2) - u * b
      q(2) = a
      do i = 3, nn
        c = p(i) - u * a - v * b
        q(i) = c
        b = a
        a = c
      end do

      return
      end
      subroutine quad ( a, b1, c, sr, si, lr, li )

c*********************************************************************72
c
cc quad() calculates the zeros of A*Z^2+B1*Z+C.
c
c  Discussion:
c
c    The quadratic formula, modified to avoid
c    overflow, is used to find the larger zero if the
c    zeros are real and both zeros are complex.
c
c    The smaller real zero is found directly from the
c    product of the zeros c/a.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    16 March 2024
c
c  Author:
c
c    Original Fortran77 version by Michael Jenkins.
c    This version by John Burkardt.
c
c  Reference:
c
c    Michael Jenkins,
c    Algorithm 493: Zeros of a Real Polynomial,
c    ACM Transactions on Mathematical Software,
c    Volume 1, Number 2, June 1975, pages 178-189.
c
c  Input:
c
c    double precision a, b1, c: the coefficients of a quadratic polynomial.
c
c  Output:
c
c    double precision sr, si: the real and imaginary parts of the smaller zero.
c
c    double precision lr, li: the real and imaginary parts of the larger zero.
c
      implicit none

      double precision a
      double precision b
      double precision b1
      double precision c
      double precision d
      double precision e
      double precision li
      double precision lr
      double precision si
      double precision sr

      if ( a == 0.0 ) then
        if ( b1 == 0.0 ) then
          sr = 0.0
        else
          sr = - c / b1
        end if
        lr = 0.0
        si = 0.0
        li = 0.0
        return
      end if

      if ( c == 0.0 ) then
        sr = 0.0
        lr = - b1 / a
        si = 0.0
        li = 0.0
        return
      end if
c
c  Compute discriminant avoiding overflow.
c
      b = b1 / 2.0

      if ( abs ( c ) <= abs ( b ) ) then
        e = 1.0 - ( a / b ) * ( c / b )
        d = sqrt ( abs ( e ) ) * abs ( b )
      else
        if ( c < 0.0 ) then
          e = b * ( b / abs ( c ) ) + a
        else
          e = b * ( b / abs ( c ) ) - a
        end if
        d = sqrt ( abs ( e ) ) * sqrt ( abs ( c ) )
      end if
c
c  Real zeros.
c
      if ( 0.0 <= e ) then

        if ( 0.0 <= b ) then
          d = - d
        end if
        lr = ( - b + d ) / a
        sr = 0.0
        if ( lr /= 0.0 ) then
          sr = ( c / lr ) / a
        end if
        si = 0.0
        li = 0.0
c
c  Complex conjugate zeros.
c
      else

        sr = - b / a
        lr = sr
        si = abs ( d / a )
        li = - si

      end if

      return
      end