PROGRAM poly20
! Generate the coefficients of the polynomial:
!    (x-1)(x-2)(x-3) ... (x-20)
! and then solve it.

USE poly_zeroes
IMPLICIT NONE

INTEGER           :: i, j, n, nitmax, iter
REAL (KIND=dp)    :: p(0:20), new_root, eps, big, small, radius(20), pr
COMPLEX (KIND=dp) :: root(20), poly(0:20)
LOGICAL           :: err(21)

INTERFACE
  SUBROUTINE sort(n, x, y, e)
    USE poly_zeroes
    IMPLICIT NONE
    INTEGER, INTENT(IN)               :: n
    COMPLEX (KIND=dp), INTENT(IN OUT) :: x(:)
    REAL (KIND=dp), INTENT(IN OUT)    :: y(:)
    LOGICAL, INTENT(IN OUT)           :: e(:)
  END SUBROUTINE sort
END INTERFACE

eps    = EPSILON(1.0D0)
small  = TINY(1.0D0)
big    = HUGE(1.0D0)
nitmax = 30

p(0) = -1.0_dp
p(1) =  1.0_dp
n = 20
DO i = 2, n
  new_root = DBLE(i)
  p(i) = p(i-1)
  DO j = i-1, 1, -1
    p(j) = p(j-1) - new_root * p(j)
  END DO
  p(0) = - new_root * p(0)
END DO

poly = CMPLX(p, KIND=dp)

CALL polzeros (n, poly, eps, big, small, nitmax, root, radius, err, iter)
CALL sort(n, root, radius, err)
OPEN(UNIT=2, FILE='poly20.out', STATUS='UNKNOWN', ACTION='WRITE')
DO i = 1, n
  pr = ABS(root(i))
  IF(pr /= 0.) pr=radius(i)/pr
  IF(radius(i) <= 0) pr = -1
  WRITE(*,300) i, err(i), root(i), radius(i), pr
  WRITE(2,300) i, err(i), root(i), radius(i), pr
END DO
WRITE(*,*)' ITERATIONS =', iter
WRITE(2,*)' ITERATIONS =', iter

STOP
300 FORMAT(' ',i4,' ERR=',l1,' X=',e21.15,',',e21.15,' R=',e7.1,' RR=',e7.1)

END PROGRAM poly20



!************************************************************************
!                             SUBROUTINE SORT                           *
!************************************************************************
!   SORT  the vector X, according to nonincreasing real parts,          *
!   the same permutation is applied to vectors Y and E.                 *
!************************************************************************
SUBROUTINE sort(n, x, y, e)
USE poly_zeroes
IMPLICIT NONE
INTEGER, INTENT(IN)               :: n
COMPLEX (KIND=dp), INTENT(IN OUT) :: x(:)
REAL (KIND=dp), INTENT(IN OUT)    :: y(:)
LOGICAL, INTENT(IN OUT)           :: e(:)

! Local variables
INTEGER           :: k, i, imax
COMPLEX (KIND=dp) :: temp
REAL (KIND=dp)    :: yt, amax, rxi
LOGICAL           :: et

DO k = 1, n-1
  amax = REAL(x(k))
  imax = k
  DO i = k+1,n
    rxi = REAL(x(i))
    IF (amax < rxi) THEN
      amax = rxi
      imax = i
    END IF
  END DO
  temp = x(k)
  x(k) = x(imax)
  x(imax) = temp
  yt = y(k)
  et = e(k)
  y(k) = y(imax)
  y(imax) = yt
  e(k) = e(imax)
  e(imax) = et
END DO
RETURN
END SUBROUTINE sort