module POLY_ZEROES
use QUADRUPLE_PRECISION
implicit none
private

private          :: ABERTH, NEWTON, START, CNVEX, CMERGE, LEFT, RIGHT, CTEST
public           :: POLZEROS

contains

!************************************************************************
!    NUMERICAL COMPUTATION OF THE ROOTS OF A POLYNOMIAL HAVING          *
!        COMPLEX COEFFICIENTS, BASED ON ABERTH'S METHOD.                *
!                      Version 1.4, June   1996                         *
!    (D. Bini, Dipartimento di Matematica, Universita' di Pisa)         *
!                         (bini@dm.unipi.it)                            *
!
!  ***************************************************************************
!  * All the software  contained in this library  is protected by copyright. *
!  * Permission  to use, copy, modify, and  distribute this software for any *
!  * purpose without fee is hereby granted, provided that this entire notice *
!  * is included  in all copies  of any software which is or includes a copy *
!  * or modification  of this software  and in all copies  of the supporting *
!  * documentation for such software.                                        *
!  ***************************************************************************
!  * THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMPLIED *
!  * WARRANTY. IN NO EVENT, NEITHER  THE AUTHORS, NOR THE PUBLISHER, NOR ANY *
!  * MEMBER  OF THE EDITORIAL BOARD OF  THE JOURNAL  "NUMERICAL ALGORITHMS", *
!  * NOR ITS EDITOR-IN-CHIEF, BE  LIABLE FOR ANY ERROR  IN THE SOFTWARE, ANY *
!  * MISUSE  OF IT  OR ANY DAMAGE ARISING OUT OF ITS USE. THE ENTIRE RISK OF *
!  * USING THE SOFTWARE LIES WITH THE PARTY DOING SO.                        *
!  ***************************************************************************
!  * ANY USE  OF THE SOFTWARE  CONSTITUTES  ACCEPTANCE  OF THE TERMS  OF THE *
!  * ABOVE STATEMENT.                                                        *
!  ***************************************************************************
!
!   AUTHOR:
!
!       DARIO ANDREA BINI
!       UNIVERSITY OF PISA, ITALY
!       E-MAIL: bini@dm.unipi.it
!
!   REFERENCE:
!
!    -  NUMERICAL COMPUTATION OF POLYNOMIAL ZEROS BY MEANS OF ABERTH'S METHOD
!       NUMERICAL ALGORITHMS, 13 (1996), PP. 179-200
!
!   Work performed under the support of the ESPRIT BRA project 6846 POSSO *
!  *************************************************************************
!
!   This version, which is compatible with Lahey's free ELF90 compiler,
!   is by Alan Miller, CSIRO Mathematical & Information Sciences,
!   Private Bag 10, Clayton South MDC, Victoria, Australia 3169.
!   Alan.Miller @ mel.dms.csiro.au    http://www.ozemail.com.au/~milleraj
!   Latest revision of ELF90 version - 19 January 1998
!
!  *************************************************************************
!
!              QUADRUPLE - PRECISION VERSION
!
!************************************************************************
!**********         SUBROUTINES AND FUNCTIONS                 ***********
!************************************************************************
!  The following routines/functions are listed:                         *
!  POLZEROS  :  computes polynomial roots by means of Aberth's method   *
!    ABERTH  :  computes the Aberth correction                          *
!    NEWTON  :  computes p(x)/p'(x) by means of Ruffini-Horner's rule   *
!    START   :  Selects N starting points by means of Rouche's theorem  *
!    CNVEX   :  Computes the convex hull, used by START                 *
!    CMERGE  :  Used by CNVEX                                           *
!    LEFT    :  Used by CMERGE                                          *
!    RIGHT   :  Used by CMERGE                                          *
!    CTEST   :  Convexity test, Used by CMERGE                          *
!************************************************************************
!                                                                       *
!                                                                       *
!************************************************************************
!********************** SUBROUTINE POLZEROS *****************************
!************************************************************************
!                        GENERAL COMMENTS                               *
!************************************************************************
!  This routine approximates the roots of   the  polynomial             *
!  p(x)=a(n+1)x^n+a(n)x^(n-1)+...+a(1), a(j)=cr(j)+I ci(j), I**2=-1,    *
!  where a(1) and a(n+1) are nonzero.                                   *
!  The coefficients are complex*16 numbers. The routine is fast, robust *
!  against overflow, and allows to deal with polynomials of any degree. *
!  Overflow situations are very unlikely and may occurr if there exist  *
!  simultaneously coefficients of moduli close to BIG and close to      *
!  SMALL, i.e., the greatest and the smallest positive real*8 numbers,  *
!  respectively. In this limit situation the program outputs a warning  *
!  message. The computation can be speeded up by performing some side   *
!  computations in single precision, thus slightly reducing the         *
!  robustness of the program (see the comments in the routine ABERTH).  *
!  Besides a set of approximations to the roots, the program delivers a *
!  set of a-posteriori error bounds which are guaranteed in the most    *
!  part of cases. In the situation where underflow does not allow to    *
!  compute a guaranteed bound, the program outputs a warning message    *
!  and sets the bound to 0. In the situation where the root cannot be   *
!  represented as a complex*16 number the error bound is set to -1.     *
!************************************************************************
!  The computation is performed by means of Aberth's method             *
!  according to the formula                                             *
!           x(i)=x(i)-newt/(1-newt*abcorr), i=1,...,n             (1)   *
!  where newt=p(x(i))/p'(x(i)) is the Newton correction and abcorr=     *
!  =1/(x(i)-x(1))+...+1/(x(i)-x(i-1))+1/(x(i)-x(i+1))+...+1/(x(i)-x(n)) *
!  is the Aberth correction to the Newton method.                       *
!************************************************************************
!  The value of the Newton correction is computed by means of the       *
!  synthetic division algorithm (Ruffini-Horner's rule) if |x|<=1,      *
!  otherwise the following more robust (with respect to overflow)       *
!  formula is applied:                                                  *
!                    newt=1/(n*y-y**2 R'(y)/R(y))                 (2)   *
!  where                                                                *
!                    y=1/x                                              *
!                    R(y)=a(1)*y**n+...+a(n)*y+a(n+1)            (2')   *
!  This computation is performed by the routine NEWTON.                 *
!************************************************************************
!  The starting approximations are complex numbers that are             *
!  equispaced on circles of suitable radii. The radius of each          *
!  circle, as well as the number of roots on each circle and the        *
!  number of circles, is determined by applying Rouche's theorem        *
!  to the functions a(k+1)*x**k and p(x)-a(k+1)*x**k, k=0,...,n.        *
!  This computation is performed by the routine START.                  *
!************************************************************************
!                              STOP CONDITION                           *
!************************************************************************
! If the condition                                                      *
!                     |p(x(j))|<EPS s(|x(j)|)                      (3)  *
! is satisfied,    where      s(x)=s(1)+x*s(2)+...+x**n * s(n+1),       *
! s(i)=|a(i)|*(1+3.8*(i-1)),  EPS is the machine precision (EPS=2**-53  *
! for the IEEE arithmetic), then the approximation x(j) is not updated  *
! and the subsequent iterations (1)  for i=j are skipped.               *
! The program stops if the condition (3) is satisfied for j=1,...,n,    *
! or if the maximum number NITMAX of  iterations   has   been reached.  *
! The condition (3) is motivated by a backward rounding error analysis  *
! of the Ruffini-Horner rule, moreover the condition (3) guarantees     *
! that the computed approximation x(j) is an exact root of a slightly   *
! perturbed polynomial.                                                 *
!************************************************************************
!             INCLUSION DISKS, A-POSTERIORI ERROR BOUNDS                *
!************************************************************************
! For each approximation x of a root, an a-posteriori absolute error    *
! bound r is computed according to the formula                          *
!                   r=n(|p(x)|+EPS s(|x|))/|p'(x)|                 (4)  *
! This provides an inclusion disk of center x and radius r containing a *
! root.                                                                 *
!************************************************************************
!************************************************************************
!*************       MEANING OF THE INPUT VARIABLES         *************
!************************************************************************
!************************************************************************
!                                                                       *
!  -- N     : degree of the polynomial.                                 *
!  -- POLY  : complex vector of N+1 components, POLY(i) is the          *
!           coefficient of x**(i-1), i=1,...,N+1 of the polynomial p(x) *
!  -- EPS   : machine precision of the floating point arithmetic used   *
!            by the computer, EPS=2**(-53)  for the IEEE standard.      *
!  -- BIG   : the max real*8, BIG=2**1023 for the IEEE standard.        *
!  -- SMALL : the min positive real*8, SMALL=2**(-1074) for the IEEE.   *
!  -- NITMAX: the max number of allowed iterations.                     *
!************************************************************************
!************************************************************************
!*************      MEANING OF THE OUTPUT VARIABLES         *************
!************************************************************************
!************************************************************************
!  ROOT   : complex vector of N components, containing the              *
!           approximations to the roots of p(x).                        *
!  RADIUS : real vector of N components, containing the error bounds to *
!           the approximations of the roots, i.e. the disk of center    *
!           ROOT(i) and radius RADIUS(i) contains a root of p(x), for   *
!           i=1,...,N. RADIUS(i) is set to -1 if the corresponding root *
!           cannot be represented as floating point due to overflow or  *
!           underflow.                                                  *
!  ERR    : vector of (N+1) components detecting an error condition;    *
!           ERR(j)=.TRUE. if after NITMAX iterations the stop condition *
!                         (3) is not satisfied for x(j)=ROOT(j);        *
!           ERR(j)=.FALSE.  otherwise, i.e., the root is reliable,      *
!                         i.e., it can be viewed as an exact root of a  *
!                         slightly perturbed polynomial.                *
!           The vector ERR is used also in the routine convex hull for  *
!           storing the abscissae of the vertices of the convex hull.   *
!           ERR(N+1) is only used for the convex hull.                  *
!  ITER   : number of iterations peformed.                              *
!************************************************************************
!************************************************************************
!************    MEANING OF THE AUXILIARY VARIABLES         *************
!************************************************************************
!************************************************************************
!  APOLY  : real vector of N+1 components used to store the moduli of   *
!           the coefficients of p(x) and the coefficients of s(x) used  *
!           to test the stop condition (3).                             *
!  APOLYR : real vector of N+1 components used to test the stop         *
!           condition                                                   *
!************************************************************************
!*****         WARNING:   2 is the output unit                    *******
!************************************************************************

subroutine POLZEROS (N, POLY, EPS, BIG, SMALL, NITMAX, ROOT, RADIUS, ERR, NITER)
integer, intent(in)                    :: N, NITMAX
integer, intent(out)                   :: NITER
type (QC), intent(in), dimension(:)    :: POLY
type (QC), intent(out), dimension(:)   :: ROOT
type (QUAD), intent(out), dimension(:) :: RADIUS
real (kind=DP), intent(in)             :: EPS, SMALL, BIG
logical, intent(out), dimension(:)     :: ERR

! Local variables
integer                     :: I, ITER, NZEROS
type (QC)                   :: CORR, ABCORR
type (QUAD)                 :: AMAX
type (QUAD), dimension(N+1) :: APOLY, APOLYR
type (QUAD), parameter      :: ZERO = QUAD( 0.0_DP, 0.0_DP ),  &
                               ONE = QUAD( 1.0_DP, 0.0_DP )

! Check consistency of data
if (abs(POLY(N+1)) == ZERO) then
  write(unit=*,fmt=*) "Inconsistent data: the leading coefficient is zero"
  stop
end if
if (abs(POLY(1)) == ZERO) then
  write(unit=*,fmt=*) "The constant term is zero: deflate the polynomial"
  stop
end if

! Compute the moduli of the coefficients
AMAX = ZERO
do I = 1, N+1
  APOLY(I) = abs(POLY(I))
  if (APOLY(I) > AMAX) then
    AMAX = APOLY(I)
  end if
  APOLYR(I) = APOLY(I)
end do
if (AMAX%HI >= BIG/(N+1)) then
  write(unit=*,fmt=*)"WARNING: COEFFICIENTS TOO BIG, OVERFLOW IS LIKELY"
  write(unit=2,fmt=*)"WARNING: COEFFICIENTS TOO BIG, OVERFLOW IS LIKELY"
end if
! Initialize
do I = 1, N
  RADIUS(I) = ZERO
  ERR(I) = .true.
end do

! Select the starting points
call START(N, APOLYR, ROOT, RADIUS, NZEROS, SMALL, BIG, ERR)

! Compute the coefficients of the backward-error polynomial
do I = 1, N+1
  APOLYR(N-I+2) = EPS*APOLY(I)*(3.8_DP*(N-I+1) + 1.0_DP)
  APOLY(I) = EPS*APOLY(I)*(3.8_DP*(I-1) + 1.0_DP)
end do
if ((APOLY(1)%HI == 0.0_DP) .or. (APOLY(N+1)%HI == 0.0_DP)) then
  write(unit=*,fmt=*)"WARNING: THE COMPUTATION OF SOME INCLUSION RADIUS"
  write(unit=*,fmt=*)"MAY FAIL. THIS IS REPORTED BY RADIUS = 0"
  write(unit=2,fmt=*)"WARNING: THE COMPUTATION OF SOME INCLUSION RADIUS"
  write(unit=2,fmt=*)"MAY FAIL. THIS IS REPORTED BY RADIUS = 0"
end if
do I = 1,N
  ERR(I) = .true.
  if (abs( RADIUS(I)%HI + 1.0_DP) < EPS) then
    ERR(I) = .false.
  end if
end do

! Starts Aberth's iterations
do ITER = 1, NITMAX
  do I = 1, N
    if (ERR(I)) then
      call NEWTON(N, POLY, APOLY, APOLYR, ROOT(I), SMALL, RADIUS(I), CORR, &
                  ERR(I))
      if (ERR(I)) then
        call ABERTH(N, I, ROOT, ABCORR)
        ROOT(I) = ROOT(I) - CORR / (QC(ONE,ZERO) - CORR*ABCORR)
      else
        NZEROS = NZEROS + 1
        if (NZEROS == N) then
          NITER = ITER
          return
        end if
      end if
    end if
  end do
end do

NITER = ITER
return
end subroutine POLZEROS




!************************************************************************
!                             SUBROUTINE NEWTON                         *
!************************************************************************
! Compute  the Newton's correction, the inclusion radius (4) and checks *
! the stop condition (3)                                                *
!************************************************************************
! Input variables:                                                      *
!     N     : degree of the polynomial p(x)                             *
!     POLY  : coefficients of the polynomial p(x)                       *
!     APOLY : upper bounds on the backward perturbations on the         *
!             coefficients of p(x) when applying Ruffini-Horner's rule  *
!     APOLYR: upper bounds on the backward perturbations on the         *
!             coefficients of p(x) when applying (2), (2')              *
!     Z     : value at which the Newton correction is computed          *
!     SMALL : the min positive TYPE (quad) ::, SMALL=2**(-1074) for the IEEE.
!************************************************************************
! Output variables:                                                     *
!     RADIUS: upper bound to the distance of Z from the closest root of *
!             the polynomial computed according to (4).                 *
!     CORR  : Newton's correction                                       *
!     AGAIN : this variable is .true. if the computed value p(z) is     *
!             reliable, i.e., (3) is not satisfied in Z. AGAIN is       *
!             .false., otherwise.                                       *
!************************************************************************

subroutine NEWTON(N, POLY, APOLY, APOLYR, Z, SMALL, RADIUS, CORR, AGAIN)
integer, intent(in)                   :: N
type (QC), intent(in), dimension(:)   :: POLY
type (QC), intent(in)                 :: Z
type (QC), intent(out)                :: CORR
type (QUAD), intent(in), dimension(:) :: APOLY, APOLYR
real (kind=DP), intent(in)            :: SMALL
type (QUAD), intent(out)              :: RADIUS
logical, intent(out)                  :: AGAIN

! Local variables
integer                :: I
type (QC)              :: P, P1, ZI, DEN, PPSP
type (QUAD)            :: AP, AZ, AZI, ABSP
type (QUAD), parameter :: ZERO = QUAD( 0.0_DP, 0.0_DP ),   &
                          ONE = QUAD( 1.0_DP, 0.0_DP )

AZ = abs(Z)
! If |z|<=1 then apply Ruffini-Horner's rule for p(z)/p'(z)
! and for the computation of the inclusion radius
if (AZ%HI <= 1.0_DP) then
  P = POLY(N+1)
  AP = APOLY(N+1)
  P1 = P
  do I = N, 2, -1
    P = P*Z + POLY(I)
    P1 = P1*Z + P
    AP = AP*AZ + APOLY(I)
  end do
  P = P*Z + POLY(1)
  AP = AP*AZ + APOLY(1)
  CORR = P/P1
  ABSP = abs(P)
  AP = abs(AP)
  AGAIN = (ABSP > QUAD(SMALL, 0.0_DP) + AP)
  if (.not.AGAIN) then
    RADIUS = real(N,kind=DP)*(ABSP + AP) / abs(P1)
  end if
  return
else
! If |z| > 1 then apply Ruffini-Horner's rule to the reversed polynomial
! and use formula (2) for p(z)/p'(z). Analogously do for the inclusion radius.
  ZI = QC(ONE,ZERO) / Z
  AZI = ONE / AZ
  P = POLY(1)
  P1 = P
  AP = APOLYR(N+1)
  do I = N, 2, -1
    P = P*ZI + POLY(N-I+2)
    P1 = P1*ZI + P
    AP = AP*AZI + APOLYR(I)
  end do
  P = P*ZI + POLY(N+1)
  AP = AP*AZI + APOLYR(1)
  ABSP = abs(P)
  AGAIN = (ABSP > QUAD(SMALL, 0.0_DP) + AP)
  PPSP = (P*Z)/P1
  DEN = QC(QUAD(real(N,kind=DP), 0.0_DP), ZERO) * PPSP - QC(ONE,ZERO)
  CORR = Z*(PPSP/DEN)
  if (AGAIN) then
    return
  end if
  RADIUS = abs(PPSP) + (AP*AZ) / abs(P1)
  RADIUS = N*RADIUS/abs(DEN)
  RADIUS = RADIUS*AZ
end if
return
end subroutine NEWTON



!************************************************************************
!                             SUBROUTINE ABERTH                         *
!************************************************************************
! Compute  the Aberth correction. To save time, the reciprocation of    *
! ROOT(J)-ROOT(I) could be performed in single precision (complex*8)    *
! In principle this might cause overflow if both ROOT(J) and ROOT(I)    *
! have too small moduli.                                                *
!************************************************************************
! Input variables:                                                      *
!     N     : degree of the polynomial                                  *
!     ROOT  : vector containing the current approximations to the roots *
!     J     : index of the component of ROOT with respect to which the  *
!             Aberth correction is computed                             *
!************************************************************************
! Output variable:                                                      *
!     ABCORR: Aberth's correction (compare (1))                         *
!************************************************************************

subroutine ABERTH(N, J, ROOT, ABCORR)
integer, intent(in)                 :: N, J
type (QC), intent(in), dimension(:) :: ROOT
type (QC), intent(out)              :: ABCORR

! Local variables
integer                :: I
type (QC)              :: Z, ZJ
type (QUAD), parameter :: ZERO = QUAD( 0.0_DP, 0.0_DP ),    &
                          ONE = QUAD( 1.0_DP, 0.0_DP )

ABCORR = cmplx(ZERO, ZERO)
ZJ = ROOT(J)
do I = 1, J-1
  Z = ZJ - ROOT(I)
  ABCORR = ABCORR + QC(ONE,ZERO) / Z
end do
do I = J+1, N
  Z = ZJ - ROOT(I)
  ABCORR = ABCORR + QC(ONE,ZERO) / Z
end do
return
end subroutine ABERTH



!************************************************************************
!                             SUBROUTINE START                          *
!************************************************************************
! Compute  the starting approximations of the roots                     *
!************************************************************************
! Input variables:                                                      *
!     N     :  number of the coefficients of the polynomial             *
!     A     :  moduli of the coefficients of the polynomial             *
!     SMALL : the min positive REAL (dp), SMALL=2**(-1074) for the IEEE. *
!     BIG   : the max REAL (dp), BIG=2**1023 for the IEEE standard.      *
! Output variables:                                                     *
!     Y     :  starting approximations                                  *
!     RADIUS:  if a component is -1 then the corresponding root has a   *
!              too big or too small modulus in order to be represented  *
!              as double float with no overflow/underflow               *
!     NZ    :  number of roots which cannot be represented without      *
!              overflow/underflow                                       *
! Auxiliary variables:                                                  *
!     H     :  needed for the computation of the convex hull            *
!************************************************************************
! This routines selects starting approximations along circles center at *
! 0 and having suitable radii. The computation of the number of circles *
! and of the corresponding radii is performed by computing the upper    *
! convex hull of the set (i,log(A(i))), i=1,...,n+1.                    *
!************************************************************************

subroutine START(N, A, Y, RADIUS, NZ, SMALL, BIG, H)
integer, intent(in)                       :: N
integer, intent(out)                      :: NZ
logical, intent(out), dimension(:)        :: H
type (QC), intent(out), dimension(:)      :: Y
real (kind=DP), intent(in)                :: SMALL, BIG
type (QUAD), intent(in out), dimension(:) :: A
type (QUAD), intent(out), dimension(:)    :: RADIUS

! Local variables
integer                :: I, IOLD, NZEROS, J, JJ
type (QUAD)            :: R, TH, ANG, TEMP
real (kind=DP)         :: XSMALL, XBIG
type (QUAD), parameter :: PI2 = TWOPI, SIGMA = QUAD( 0.7_DP, 0.0_DP )

XSMALL = log(SMALL)
XBIG = log(BIG)
NZ = 0
! Compute the logarithm A(I) of the moduli of the coefficients of
! the polynomial and then the upper convex hull of the set (A(I),I)
do I = 1, N+1
  if (A(I)%HI /= 0.0_DP) then
    A(I) = log(A(I))
  else
    A(I) = QUAD(-1.0e+50_DP, 0.0_DP)
  end if
end do
call CNVEX(N+1, A, H)

! Given the upper convex hull of the set (A(I),I) compute the moduli
! of the starting approximations by means of Rouche's theorem
IOLD = 1
TH = PI2 / real(N,kind=DP)
do I = 2, N+1
  if (H(I)) then
    NZEROS = I - IOLD
    TEMP = (A(IOLD) - A(I)) / real(NZEROS,kind=DP)
! Check if the modulus is too small
    if ((TEMP%HI < -XBIG) .and. (TEMP%HI >= XSMALL)) then
      write(unit=*,fmt=*)"WARNING:", NZEROS, " ZERO(S) ARE TOO SMALL TO"
      write(unit=*,fmt=*)"REPRESENT THEIR INVERSES AS TYPE (qc), THEY"
      write(unit=*,fmt=*)"ARE REPLACED BY SMALL NUMBERS, THE CORRESPONDING"
      write(unit=*,fmt=*)"RADII ARE SET TO -1"
      write(unit=2,fmt=*)"WARNING:", NZEROS, " ZERO(S) ARE TOO SMALL TO "
      write(unit=2,fmt=*)"REPRESENT THEIR INVERSES AS TYPE (qc), THEY"
      write(unit=2,fmt=*)"ARE REPLACED BY SMALL NUMBERS, THE CORRESPONDING"
      write(unit=2,fmt=*)"RADII ARE SET TO -1"
      NZ = NZ + NZEROS
      R = QUAD(1.0_DP, 0.0_DP) / BIG
    end if
    if (TEMP%HI < XSMALL) then
      NZ = NZ + NZEROS
      write(unit=*,fmt=*)"WARNING: ", NZEROS, " ZERO(S) ARE TOO SMALL TO BE"
      write(unit=*,fmt=*)"REPRESENTED AS TYPE (qc), THEY ARE SET TO 0"
      write(unit=*,fmt=*)"THE CORRESPONDING RADII ARE SET TO -1"
      write(unit=2,fmt=*)"WARNING: ", NZEROS, " ZERO(S) ARE TOO SMALL TO BE"
      write(unit=2,fmt=*)"REPRESENTED AS TYPE (qc), THEY ARE SET 0"
      write(unit=2,fmt=*)"THE CORRESPONDING RADII ARE SET TO -1"
    end if
! Check if the modulus is too big
    if (TEMP%HI > XBIG) then
      R = QUAD(BIG, 0.0_DP)
      NZ = NZ + NZEROS
      write(unit=*,fmt=*)"WARNING: ", NZEROS, " ZEROS(S) ARE TOO BIG TO BE"
      write(unit=*,fmt=*)"REPRESENTED AS TYPE (qc),"
      write(unit=*,fmt=*)"THE CORRESPONDING RADII ARE SET TO -1"
      write(unit=2,fmt=*)"WARNING: ",NZEROS, " ZERO(S) ARE TOO BIG TO BE"
      write(unit=2,fmt=*)"REPRESENTED AS TYPE (qc),"
      write(unit=2,fmt=*)"THE CORRESPONDING RADII ARE SET TO -1"
    end if
    if ((TEMP%HI <= XBIG) .and. (TEMP%HI > max(-XBIG, XSMALL))) then
      R = exp(TEMP)
    end if
! Compute NZEROS approximations equally distributed in the disk of radius R
    ANG = PI2 / real(NZEROS,kind=DP)
    do J = IOLD, I-1
      JJ = J - IOLD + 1
      if (R%HI <= 1.0_DP/BIG .or. R == QUAD(BIG, 0.0_DP)) then
        RADIUS(J) = QUAD(-1.0_DP, 0.0_DP)
      end if
      Y(J)%QR = R * cos(ANG*JJ + TH*I + SIGMA)
      Y(J)%QI = R * sin(ANG*JJ + TH*I + SIGMA)
    end do
    IOLD = I
  end if
end do
return
end subroutine START


!************************************************************************
!                             SUBROUTINE CNVEX                          *
!************************************************************************
! Compute  the upper convex hull of the set (i,a(i)), i.e., the set of  *
! vertices (i_k,a(i_k)), k=1,2,...,m, such that the points (i,a(i)) lie *
! below the straight lines passing through two consecutive vertices.    *
! The abscissae of the vertices of the convex hull equal the indices of *
! the TRUE  components of the logical output vector H.                  *
! The used method requires O(nlog n) comparisons and is based on a      *
! divide-and-conquer technique. Once the upper convex hull of two       *
! contiguous sets  (say, {(1,a(1)),(2,a(2)),...,(k,a(k))} and           *
! {(k,a(k)), (k+1,a(k+1)),...,(q,a(q))}) have been computed, then       *
! the upper convex hull of their union is provided by the subroutine    *
! CMERGE. The program starts with sets made up by two consecutive       *
! points, which trivially constitute a convex hull, then obtains sets   *
! of 3,5,9... points,  up to  arrive at the entire set.                 *
! The program uses the subroutine  CMERGE; the subroutine CMERGE uses   *
! the subroutines LEFT, RIGHT and CTEST. The latter tests the convexity *
! of the angle formed by the points (i,a(i)), (j,a(j)), (k,a(k)) in the *
! vertex (j,a(j)) up to within a given tolerance TOLER, where i<j<k.    *
!************************************************************************

subroutine CNVEX(N, A, H)
integer, intent(in)                   :: N
logical, intent(out), dimension(:)    :: H
type (QUAD), intent(in), dimension(:) :: A

! Local variables
integer :: I, J, K, M, NJ, JC

H(1:N) = .true.

! compute K such that N-2 <= 2**K < N-1
K = int(log(N - 2.0_DP) / log(2.0_DP))
if (2**(K+1) <= (N-2)) then
  K = K+1
end if

! For each M=1,2,4,8,...,2**K, consider the NJ pairs of consecutive
! sets made up by M+1 points having the common vertex
! (JC,A(JC)), where JC=M*(2*J+1)+1 and J=0,...,NJ,
! NJ = MAX(0, INT((N-2-M)/(M+M))).
! Compute the upper convex hull of their union by means of subroutine CMERGE
M = 1
do I = 0, K
  NJ = max(0, int((N-2-M)/(M+M)))
  do J = 0, NJ
    JC = (J+J+1)*M + 1
    call CMERGE(N, A, JC, M, H)
  end do
  M = M + M
end do
return
end subroutine CNVEX



!************************************************************************
!                             SUBROUTINE LEFT                           *
!************************************************************************
! Given as input the integer I and the vector H of logical, compute the *
! the maximum integer IL such that IL < I and H(IL) is TRUE.            *
!************************************************************************
! Input variables:                                                      *
!     H   : vector of logical                                           *
!     I   : integer                                                     *
!************************************************************************
! Output variable:                                                      *
!     IL  : maximum integer such that IL < I, H(IL)=.TRUE.              *
!************************************************************************

subroutine LEFT(H, I, IL)
integer, intent(in)               :: I
integer, intent(out)              :: IL
logical, intent(in), dimension(:) :: H

integer :: L

do L = I-1, 0, -1
  if (H(L)) then
    IL = L
    return
  end if
end do
return
end subroutine LEFT



!************************************************************************
!                             SUBROUTINE RIGHT                          *
!************************************************************************
!************************************************************************
! Given as input the integer I and the vector H of logical, compute the *
! the minimum integer IR such that IR>I and H(IL) is TRUE.              *
!************************************************************************
!************************************************************************
! Input variables:                                                      *
!     N   : length of the vector H                                      *
!     H   : vector of logical                                           *
!     I   : integer                                                     *
!************************************************************************
! Output variable:                                                      *
!     IR  : minimum integer such that IR>I, H(IR)=.TRUE.                *
!************************************************************************

subroutine RIGHT(N, H, I, IR)
integer, intent(in)               :: N, I
integer, intent(out)              :: IR
logical, intent(in), dimension(:) :: H

integer :: R

do R = I+1, N
  if (H(R)) then
    IR = R
    return
  end if
end do
return
end subroutine RIGHT



!************************************************************************
!                             SUBROUTINE CMERGE                         *
!************************************************************************
! Given the upper convex hulls of two consecutive sets of pairs         *
! (j,A(j)), compute the upper convex hull of their union                *
!************************************************************************
! Input variables:                                                      *
!     N    : length of the vector A                                     *
!     A    : vector defining the points (j,A(j))                        *
!     I    : abscissa of the common vertex of the two sets              *
!     M    : the number of elements of each set is M+1                  *
!************************************************************************
! Input/Output variable:                                                *
!     H    : vector defining the vertices of the convex hull, i.e.,     *
!            H(j) is .TRUE. if (j,A(j)) is a vertex of the convex hull  *
!            This vector is used also as output.                        *
!************************************************************************

subroutine CMERGE(N, A, I, M, H)
integer, intent(in)                   :: N, M, I
logical, intent(in out), dimension(:) :: H
type (QUAD), intent(in), dimension(:) :: A

! Local variables
integer :: IR, IL, IRR, ILL
logical :: TSTL, TSTR

! at the left and the right of the common vertex (I,A(I)) determine
! the abscissae IL,IR, of the closest vertices of the upper convex
! hull of the left and right sets, respectively
call LEFT(H, I, IL)
call RIGHT(N, H, I, IR)

! check the convexity of the angle formed by IL,I,IR
if (CTEST(A, IL, I, IR)) then
  return
else
! continue the search of a pair of vertices in the left and right
! sets which yield the upper convex hull
  H(I) = .false.
  do
    if (IL == (I-M)) then
      TSTL = .true.
    else
      call LEFT(H, IL, ILL)
      TSTL = CTEST(A, ILL, IL, IR)
    end if
    if (IR == min(N, I+M)) then
      TSTR = .true.
    else
      call RIGHT(N, H, IR, IRR)
      TSTR = CTEST(A, IL, IR, IRR)
    end if
    H(IL) = TSTL
    H(IR) = TSTR
    if (TSTL .and. TSTR) then
      return
    end if
    if (.not.TSTL) then
      IL = ILL
    end if
    if (.not.TSTR) then
      IR = IRR
    end if
  end do
end if

return
end subroutine CMERGE



!************************************************************************
!                             FUNCTION CTEST                            *
!************************************************************************
! Test the convexity of the angle formed by (IL,A(IL)), (I,A(I)),       *
! (IR,A(IR)) at the vertex (I,A(I)), up to within the tolerance         *
! TOLER. If convexity holds then the function is set to .TRUE.,         *
! otherwise CTEST=.FALSE. The parameter TOLER is set to 0.4 by default. *
!************************************************************************
! Input variables:                                                      *
!     A       : vector of double                                        *
!     IL,I,IR : integers such that IL < I < IR                          *
!************************************************************************
! Output:                                                               *
!     .TRUE. if the angle formed by (IL,A(IL)), (I,A(I)), (IR,A(IR)) at *
!            the vertex (I,A(I)), is convex up to within the tolerance  *
!            TOLER, i.e., if                                            *
!            (A(I)-A(IL))*(IR-I) - (A(IR)-A(I))*(I-IL) > TOLER.         *
!     .FALSE.,  otherwise.                                              *
!************************************************************************

function CTEST(A, IL, I, IR) result(ok)
integer, intent(in)                   :: I, IL, IR
type (QUAD), intent(in), dimension(:) :: A
logical                               :: ok

! Local variables
type (QUAD)               :: S1, S2
real (kind=DP), parameter :: TOLER = 0.4_DP

S1 = A(I) - A(IL)
S2 = A(IR) - A(I)
S1 = S1*(IR-I)
S2 = S2*(I-IL)
ok = .false.
if (S1 > (S2 + TOLER)) then
  ok = .true.
end if
return
end function CTEST


end module POLY_ZEROES