module SORTX
use QUADRUPLE_PRECISION
implicit none
private
public :: SORT

contains

!************************************************************************
!                             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)

integer, intent(in)                       :: N
type (QC), intent(in out), dimension(:)   :: X
type (QUAD), intent(in out), dimension(:) :: Y
logical, intent(in out), dimension(:)     :: E

! Local variables
integer     :: K, I, IMAX
type (QC)   :: TEMP
type (QUAD) :: YT, AMAX, RXI
logical     :: ET

do K = 1, N-1
  AMAX = X(K)%QR
  IMAX = K
  do I = K+1,N
    RXI = X(I)%QR
    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

end module SORTX



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

use QUADRUPLE_PRECISION
use POLY_ZEROES
use SORTX
implicit none

integer                      :: I, J, N, NITMAX, ITER, IER
real (kind=DP)               :: EPS, BIG, SMALL
type (QUAD)                  :: NEW_ROOT, PR
type (QUAD), dimension(0:20) :: P
type (QUAD), dimension(20)   :: RADIUS
type (QC), dimension(20)     :: ROOT
type (QC), dimension(0:20)   :: POLY
logical, dimension(21)       :: ERR
character (len=36)           :: STRING
character (len=*), parameter :: fmt300 =  &
                      "(a, i3, a, l1, a, e21.15, a, e21.15, a, e7.1, a, e7.1)"

EPS    = 1.0e-29_DP
SMALL  = tiny(1.0_DP) * 1.0e+16_DP
BIG    = huge(1.0_DP)
NITMAX = 50

P(0) = QUAD(-1.0_DP, 0.0_DP)
P(1) = QUAD( 1.0_DP, 0.0_DP)
N = 20
do I = 2, N
  NEW_ROOT = QUAD( real(I,kind=DP), 0.0_DP )
  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

open(unit=2, file="poly20.out", status="UNKNOWN", action="WRITE")
write(unit=2, fmt="(a)") "Coefficients of input polynomial"
do I = 0, N
  call QUAD_STRING(P(I), STRING, IER)
  if (IER == 0) then
    write(unit=2, fmt="(i4, a, a36)") I, "  ", STRING
  end if
end do
write(unit=2, fmt=*)

do I = 0, N
  POLY(I) = cmplx(P(I), QUAD(0.0_DP, 0.0_DP) )
end do

call POLZEROS (N, POLY, EPS, BIG, SMALL, NITMAX, ROOT, RADIUS, ERR, ITER)
call SORT(N, ROOT, RADIUS, ERR)
write(unit=2, fmt="(a)") "Solution:"
write(unit=*, fmt="(t16, a, t40, a)") " Real part", "Imaginary part"
write(unit=2, fmt="(t16, a, t40, a)") " Real part", "Imaginary part"
do I = 1, N
  PR = abs(ROOT(I))
  if(PR%HI /= 0.0_DP) then
    PR = RADIUS(I) / PR
  end if
  if(RADIUS(I)%HI <= 0.0_DP) then
    PR = QUAD( -1.0_DP, 0.0_DP )
  end if
  write(unit=*, fmt=fmt300) " ", I, " ERR=", ERR(I), " X=", ROOT(I)%QR%HI, &
                            ", ", ROOT(I)%QI%HI, "  R=", RADIUS(I)%HI,  &
                            " RR=", PR%HI
  write(unit=2, fmt=fmt300) " ", I, " ERR=", ERR(I), " X=", ROOT(I)%QR%HI, &
                            ", ", ROOT(I)%QI%HI, "  R=", RADIUS(I)%HI,  &
                            " RR=", PR%HI

end do
write(unit=*, fmt=*)" ITERATIONS =", ITER
write(unit=2, fmt=*)" ITERATIONS =", ITER

stop

end program POLY20