cpoly Subroutine

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

    implicit none

    integer,intent(in) :: degree !! degree of polynomial
    real(wp), dimension(degree+1), intent(in) :: opr !! vectors of real parts of the coefficients in
                                                     !! order of decreasing powers.
    real(wp), dimension(degree+1), intent(in) :: opi !! vectors of imaginary parts of the coefficients in
                                                     !! order of decreasing powers.
    real(wp), dimension(degree), intent(out) :: zeror !! real parts of the zeros
    real(wp), dimension(degree), intent(out) :: zeroi !! imaginary parts of the zeros
    logical,intent(out) :: fail !! true only if leading coefficient is zero or if cpoly
                                !! has found fewer than `degree` zeros.

    real(wp) :: sr , si , tr , ti , pvr , pvi, xxx , zr , zi , bnd , xx , yy
    real(wp), dimension(degree+1) :: pr , pi , hr, hi , qpr, qpi, qhr , qhi , shr , shi
    logical :: conv
    integer :: cnt1 , cnt2, i , idnn2 , nn

    real(wp), parameter :: base = radix(1.0_wp)
    real(wp), parameter :: eta = eps
    real(wp), parameter :: infin = huge(1.0_wp)
    real(wp), parameter :: smalno = tiny(1.0_wp)
    real(wp), parameter :: are = eta
    real(wp), parameter :: cosr = cos(94.0_wp*deg2rad)  !! -.069756474
    real(wp), parameter :: sinr = sin(86.0_wp*deg2rad)  !! .99756405
    real(wp), parameter :: mre = 2.0_wp*sqrt(2.0_wp)*eta
    real(wp), parameter :: cos45 = cos(45.0_wp*deg2rad) !! .70710678

    if ( opr(1)==0.0_wp .and. opi(1)==0.0_wp ) then
        ! algorithm fails if the leading coefficient is zero.
        fail = .true.
        return
    end if

    xx = cos45
    yy = -xx
    fail = .false.
    nn = degree + 1

    ! remove the zeros at the origin if any.
    do
        if ( opr(nn)/=0.0_wp .or. opi(nn)/=0.0_wp ) then
            exit
        else
            idnn2 = degree - nn + 2
            zeror(idnn2) = 0.0_wp
            zeroi(idnn2) = 0.0_wp
            nn = nn - 1
        endif
    end do

    ! make a copy of the coefficients.
    do i = 1 , nn
        pr(i) = opr(i)
        pi(i) = opi(i)
        shr(i) = cmod(pr(i),pi(i))
    enddo
    ! scale the polynomial.
    bnd = scale(nn,shr,eta,infin,smalno,base)
    if ( bnd/=1.0_wp ) then
        do i = 1 , nn
            pr(i) = bnd*pr(i)
            pi(i) = bnd*pi(i)
        enddo
    endif

    ! start the algorithm for one zero.
    main : do
        if ( nn>2 ) then
            ! calculate bnd, a lower bound on the modulus of the zeros.
            do i = 1 , nn
                shr(i) = cmod(pr(i),pi(i))
            enddo
            bnd = cauchy(nn,shr,shi)
            ! outer loop to control 2 major passes with different sequences
            ! of shifts.
            do cnt1 = 1 , 2
                ! first stage calculation, no shift.
                call noshft(5)
                ! inner loop to select a shift.
                do cnt2 = 1 , 9
                    ! shift is chosen with modulus bnd and amplitude rotated by
                    ! 94 degrees from the previous shift
                    xxx = cosr*xx - sinr*yy
                    yy = sinr*xx + cosr*yy
                    xx = xxx
                    sr = bnd*xx
                    si = bnd*yy
                    ! second stage calculation, fixed shift.
                    call fxshft(10*cnt2,zr,zi,conv)
                    if ( conv ) then
                        ! the second stage jumps directly to the third stage iteration.
                        ! if successful the zero is stored and the polynomial deflated.
                        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)
                        enddo
                        cycle main
                    endif
                    ! if the iteration is unsuccessful another shift is chosen.
                enddo
                ! if 9 shifts fail, the outer loop is repeated with another
                ! sequence of shifts.
            enddo
            ! the zerofinder has failed on two major passes.
            ! return empty handed.
            fail = .true.
            return
        else
            exit
        endif

    end do main

    ! calculate the final zero and return.
    call cdivid(-pr(2),-pi(2),pr(1),pi(1),zeror(degree),zeroi(degree))

contains

subroutine noshft(l1)

    ! computes the derivative polynomial as the initial h
    ! polynomial and computes l1 no-shift h polynomials.

    implicit none

    integer,intent(in) :: l1

    integer :: i , j , jj , n , nm1
    real(wp) :: xni , t1 , t2

    n = nn - 1
    nm1 = n - 1
    do i = 1 , n
       xni = nn - i
       hr(i) = xni*pr(i)/real(n,wp)
       hi(i) = xni*pi(i)/real(n,wp)
    enddo
    do jj = 1 , l1
       if ( cmod(hr(n),hi(n))<=eta*10.0_wp*cmod(pr(n),pi(n)) ) then
          ! if the constant term is essentially zero, shift h coefficients.
          do i = 1 , nm1
             j = nn - i
             hr(j) = hr(j-1)
             hi(j) = hi(j-1)
          enddo
          hr(1) = 0.0_wp
          hi(1) = 0.0_wp
       else
          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)
          enddo
          hr(1) = pr(1)
          hi(1) = pi(1)
       endif
    enddo
end subroutine noshft

subroutine fxshft(l2,zr,zi,conv)

    ! computes l2 fixed-shift h polynomials and tests for
    ! convergence.
    ! initiates a variable-shift iteration and returns with the
    ! approximate zero if successful.

    implicit none

    integer,intent(in) :: l2 !! limit of fixed shift steps
    real(wp) :: zr , zi !! approximate zero if conv is .true.
    logical :: conv !! logical indicating convergence of stage 3 iteration

    integer :: i , j , n
    real(wp) :: otr , oti , svsr , svsi
    logical :: test , pasd , bool

    n = nn - 1
    ! evaluate p at s.
    call polyev(nn,sr,si,pr,pi,qpr,qpi,pvr,pvi)
    test = .true.
    pasd = .false.
    ! calculate first t = -p(s)/h(s).
    call calct(bool)
    ! main loop for one second stage step.
    do j = 1 , l2
       otr = tr
       oti = ti
        ! compute next h polynomial and new t.
       call nexth(bool)
       call calct(bool)
       zr = sr + tr
       zi = si + ti
        ! test for convergence unless stage 3 has failed once or this
        ! is the last h polynomial .
       if ( .not.(bool .or. .not.test .or. j==l2) ) then
          if ( cmod(tr-otr,ti-oti)>=0.5_wp*cmod(zr,zi) ) then
             pasd = .false.
          elseif ( .not.pasd ) then
             pasd = .true.
          else
            ! the weak convergence test has been passed twice, start the
            ! third stage iteration, after saving the current h polynomial
            ! and shift.
             do i = 1 , n
                shr(i) = hr(i)
                shi(i) = hi(i)
             enddo
             svsr = sr
             svsi = si
             call vrshft(10,zr,zi,conv)
             if ( conv ) return
            ! the iteration failed to converge. turn off testing and restore
            ! h,s,pv and t.
             test = .false.
             do i = 1 , n
                hr(i) = shr(i)
                hi(i) = shi(i)
             enddo
             sr = svsr
             si = svsi
             call polyev(nn,sr,si,pr,pi,qpr,qpi,pvr,pvi)
             call calct(bool)
          endif
       endif
    enddo
    ! attempt an iteration with final h polynomial from second stage.
    call vrshft(10,zr,zi,conv)

end subroutine fxshft

subroutine vrshft(l3,zr,zi,conv)

    ! carries out the third stage iteration.

    implicit none

    integer,intent(in) :: l3 !! limit of steps in stage 3.
    real(wp) :: zr , zi !! on entry contains the initial iterate, if the
                        !! iteration converges it contains the final iterate
                        !! on exit.
    logical :: conv !! .true. if iteration converges

    real(wp) :: mp , ms , omp , relstp , r1 , r2 , tp
    logical :: b , bool
    integer :: i , j

    conv = .false.
    b = .false.
    sr = zr
    si = zi

    ! main loop for stage three
    do i = 1 , l3
        ! 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>20.0_wp*errev(nn,qpr,qpi,ms,mp,are,mre) ) then
          if ( i==1 ) then
             omp = mp
          elseif ( b .or. mp<omp .or. relstp>=0.05_wp ) then
            ! exit if polynomial value increases significantly.
             if ( mp*0.1_wp>omp ) return
             omp = mp
          else
            ! iteration has stalled. probably a cluster of zeros. do 5 fixed
            ! shift steps into the cluster to force one zero to dominate.
             tp = relstp
             b = .true.
             if ( relstp<eta ) tp = eta
             r1 = sqrt(tp)
             r2 = sr*(1.0_wp+r1) - si*r1
             si = sr*r1 + si*(1.0_wp+r1)
             sr = r2
             call polyev(nn,sr,si,pr,pi,qpr,qpi,pvr,pvi)
             do j = 1 , 5
                call calct(bool)
                call nexth(bool)
             enddo
             omp = infin
          endif
            ! calculate next iterate.
          call calct(bool)
          call nexth(bool)
          call calct(bool)
          if ( .not.(bool) ) then
             relstp = cmod(tr,ti)/cmod(sr,si)
             sr = sr + tr
             si = si + ti
          endif
       else
          ! polynomial value is smaller in value than a bound on the error
          ! in evaluating p, terminate the iteration.
          conv = .true.
          zr = sr
          zi = si
          return
       endif
    enddo

end subroutine vrshft

subroutine calct(bool)

    ! computes `t = -p(s)/h(s)`.

    implicit none

    logical,intent(out) :: bool !! logical, set true if `h(s)` is essentially zero.

    real(wp) :: hvr , hvi
    integer :: n

    n = nn - 1
    ! evaluate h(s).
    call polyev(n,sr,si,hr,hi,qhr,qhi,hvr,hvi)
    bool = cmod(hvr,hvi)<=are*10.0_wp*cmod(hr(n),hi(n))
    if ( bool ) then
       tr = 0.0_wp
       ti = 0.0_wp
    else
        call cdivid(-pvr,-pvi,hvr,hvi,tr,ti)
    end if

end subroutine calct

subroutine nexth(bool)

    ! calculates the next shifted h polynomial.

    implicit none

    logical,intent(in) :: bool !! logical, if .true. `h(s)` is essentially zero

    real(wp) :: t1 , t2
    integer :: j , n , nm1

    n = nn - 1
    nm1 = n - 1
    if ( bool ) then
       ! if h(s) is zero replace h with qh.
       do j = 2 , n
          hr(j) = qhr(j-1)
          hi(j) = qhi(j-1)
       enddo
       hr(1) = 0.0_wp
       hi(1) = 0.0_wp
    else
        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)
        enddo
        hr(1) = qpr(1)
        hi(1) = qpi(1)
    end if

end subroutine nexth

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

    ! evaluates a polynomial  p  at  s  by the horner recurrence
    ! placing the partial sums in q and the computed value in pv.

    implicit none

    integer,intent(in) :: nn
    real(wp) :: pr(nn) , pi(nn) , qr(nn) , qi(nn) , sr , si , pvr , pvi

    real(wp) :: t
    integer :: i

    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
    enddo

end subroutine polyev

real(wp) function errev(nn,qr,qi,ms,mp,are,mre)

    ! bounds the error in evaluating the polynomial
    ! by the horner recurrence.

    implicit none

    integer,intent(in) :: nn
    real(wp) :: qr(nn), qi(nn) !! the partial sums
    real(wp) :: ms !! modulus of the point
    real(wp) :: mp !! modulus of polynomial value
    real(wp) :: are, mre !! error bounds on complex addition and multiplication

    integer :: i
    real(wp) :: e

    e = cmod(qr(1),qi(1))*mre/(are+mre)
    do i = 1 , nn
       e = e*ms + cmod(qr(i),qi(i))
    enddo
    errev = e*(are+mre) - mp*mre

end function errev

real(wp) function cauchy(nn,pt,q)

    ! cauchy computes a lower bound on the moduli of
    ! the zeros of a polynomial

    implicit none

    integer,intent(in) :: nn
    real(wp) :: q(nn)
    real(wp) :: pt(nn) !! the modulus of the coefficients.

    integer :: i , n
    real(wp) :: x , xm , f , dx , df

    pt(nn) = -pt(nn)
    ! compute upper estimate of bound.
    n = nn - 1
    x = exp((log(-pt(nn))-log(pt(1)))/real(n,wp))
    if ( pt(n)/=0.0_wp ) then
       ! if newton step at the origin is better, use it.
       xm = -pt(nn)/pt(n)
       if ( xm<x ) x = xm
    endif

    do
        ! chop the interval (0,x) unitl f <= 0.
        xm = x*0.1_wp
        f = pt(1)
        do i = 2 , nn
           f = f*xm + pt(i)
        enddo
        if ( f<=0.0_wp ) then
           dx = x
           do
               ! newton iteration until x converges to two decimal places.
               if ( abs(dx/x)<=0.005_wp ) then
                  cauchy = x
                  exit
               end if
               q(1) = pt(1)
               do i = 2 , nn
                  q(i) = q(i-1)*x + pt(i)
               enddo
               f = q(nn)
               df = q(1)
               do i = 2 , n
                  df = df*x + q(i)
               enddo
               dx = f/df
               x = x - dx
           end do
           exit
        else
           x = xm
        endif
    end do

end function cauchy

real(wp) function scale(nn,pt,eta,infin,smalno,base)

    ! returns a scale factor to multiply the coefficients of the
    ! polynomial. the scaling is done to avoid overflow and to avoid
    ! undetected underflow interfering with the convergence
    ! criterion.  the factor is a power of the base.

    implicit none

    integer :: nn
    real(wp) :: pt(nn)  !! modulus of coefficients of p
    real(wp) :: eta , infin , smalno , base !! constants describing the
                                            !! floating point arithmetic.

    real(wp) :: hi , lo , max , min , x , sc
    integer :: i , l

    ! find largest and smallest moduli of coefficients.
    hi = sqrt(infin)
    lo = smalno/eta
    max = 0.0_wp
    min = infin
    do i = 1 , nn
       x = pt(i)
       if ( x>max ) max = x
       if ( x/=0.0_wp .and. x<min ) min = x
    enddo
    ! scale only if there are very large or very small components.
    scale = 1.0_wp
    if ( min>=lo .and. max<=hi ) return
    x = lo/min
    if ( x>1.0_wp ) then
       sc = x
       if ( infin/sc>max ) sc = 1.0_wp
    else
       sc = 1.0_wp/(sqrt(max)*sqrt(min))
    endif
    l = log(sc)/log(base) + 0.5_wp
    scale = base**l

end function scale

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

    ! complex division c = a/b, avoiding overflow.

    implicit none

    real(wp) :: ar , ai , br , bi , cr , ci , r , d

    if ( br==0.0_wp .and. bi==0.0_wp ) then
       ! division by zero, c = infinity.
       cr = infin
       ci = infin
    elseif ( abs(br)>=abs(bi) ) then
       r = bi/br
       d = br + r*bi
       cr = (ar+ai*r)/d
       ci = (ai-ar*r)/d
    else
        r = br/bi
        d = bi + r*br
        cr = (ar*r+ai)/d
        ci = (ai*r-ar)/d
    end if

end subroutine cdivid

real(wp) function cmod(r,i)

    implicit none

    real(wp) :: r , i , ar , ai

    ar = abs(r)
    ai = abs(i)
    if ( ar<ai ) then
       cmod = ai*sqrt(1.0_wp+(ar/ai)**2)
    elseif ( ar<=ai ) then
       cmod = ar*sqrt(2.0_wp)
    else
       cmod = ar*sqrt(1.0_wp+(ai/ar)**2)
    end if

end function cmod

end subroutine cpoly