MODULE Solve_Real_Poly
! CACM Algorithm 493 by Jenkins & Traub

! Code converted using TO_F90 by Alan Miller
! Date: 2000-01-08  Time: 16:02:50

! Latest revision - 10 February 2000
! amiller @ bigpond.net.au

USE quadruple_precision
IMPLICIT NONE

! COMMON /global/ p, qp, k, qk, svk, sr, si, u, v, a, b, c, d, a1,  &
!    a2, a3, a6, a7, e, f, g, h, szr, szi, lzr, lzi, eta, are, mre, n, nn

TYPE (quad), ALLOCATABLE, SAVE  :: p(:), qp(:), k(:), qk(:), svk(:)
TYPE (quad), SAVE               :: sr, si, u, v, a, b, c, d, a1, a2, a3, a6,  &
                                   a7, ee, f, g, h, szr, szi, lzr, lzi
REAL (r10), SAVE                :: eta, are, mre
INTEGER, SAVE                   :: n, nn

PRIVATE
PUBLIC  :: rpoly


CONTAINS


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

! FINDS THE ZEROS OF A REAL POLYNOMIAL
! OP  - TYPE (quad) VECTOR OF COEFFICIENTS IN ORDER OF DECREASING POWERS.
! DEGREE   - INTEGER DEGREE OF POLYNOMIAL.
! ZEROR, ZEROI - OUTPUT TYPE (quad) VECTORS OF
!                REAL AND IMAGINARY PARTS OF THE ZEROS.
! FAIL  - OUTPUT LOGICAL PARAMETER, TRUE ONLY IF THE LEADING COEFFICIENT
!         IS ZERO, OR IF RPOLY HAS FOUND FEWER THAN DEGREE ZEROS.
!         IN THE LATTER CASE DEGREE IS RESET TO THE NUMBER OF ZEROS FOUND.
! TO CHANGE THE SIZE OF POLYNOMIALS WHICH CAN BE SOLVED, RESET THE DIMENSIONS
! OF THE ARRAYS IN THE COMMON AREA AND IN THE FOLLOWING DECLARATIONS.
! THE SUBROUTINE USES SINGLE PRECISION CALCULATIONS FOR SCALING, BOUNDS AND
! ERROR CALCULATIONS.  ALL CALCULATIONS FOR THE ITERATIONS ARE DONE IN DOUBLE
! PRECISION.

TYPE (quad), INTENT(IN)   :: op(:)
INTEGER, INTENT(IN OUT)   :: degree
TYPE (quad), INTENT(OUT)  :: zeror(:)
TYPE (quad), INTENT(OUT)  :: zeroi(:)
LOGICAL, INTENT(OUT)      :: fail

TYPE (quad), ALLOCATABLE  :: temp(:)
REAL (r10), ALLOCATABLE   :: pt(:)

TYPE (quad)  :: t, aa, bb, cc, factor
REAL (r10)   :: lo, MAX, MIN, xx, yy, cosr, sinr, xxx, x, sc, bnd,  &
                xm, ff, df, dx, infin, smalno, base
INTEGER      :: cnt, nz, i, j, jj, l, nm1
LOGICAL      :: zerok

! THE FOLLOWING STATEMENTS SET MACHINE CONSTANTS USED IN VARIOUS PARTS OF THE
! PROGRAM.  THE MEANING OF THE FOUR CONSTANTS ARE...
! ETA     THE MAXIMUM RELATIVE REPRESENTATION ERROR WHICH CAN BE DESCRIBED AS
!         THE SMALLEST POSITIVE FLOATING POINT NUMBER SUCH THAT 1._r10+ETA > 1.
! INFINY  THE LARGEST FLOATING-POINT NUMBER.
! SMALNO  THE SMALLEST POSITIVE FLOATING-POINT NUMBER IF THE EXPONENT RANGE
!         DIFFERS IN SINGLE AND TYPE (quad) THEN SMALNO AND INFIN
!         SHOULD INDICATE THE SMALLER RANGE.
! BASE    THE BASE OF THE FLOATING-POINT NUMBER SYSTEM USED.

base = RADIX(0.0)
eta = EPSILON(1.0_r10) * 1.0E-18
infin = HUGE(0.0)
smalno = TINY(0.0)

! ARE AND MRE REFER TO THE UNIT ERROR IN + AND * RESPECTIVELY.
! THEY ARE ASSUMED TO BE THE SAME AS ETA.
are = eta
mre = eta
lo = smalno / eta

! INITIALIZATION OF CONSTANTS FOR SHIFT ROTATION
xx = .70710678
yy = -xx
cosr = -.069756474
sinr = .99756405
fail = .false.
n = degree
nn = n + 1

! ALGORITHM FAILS IF THE LEADING COEFFICIENT IS ZERO.
IF (op(1)%hi == 0._r10) THEN
  fail = .true.
  degree = 0
  RETURN
END IF

! REMOVE THE ZEROS AT THE ORIGIN IF ANY
10 IF (op(nn)%hi == 0.0_r10) THEN
  j = degree - n + 1
  zeror(j) = 0._r10
  zeroi(j) = 0._r10
  nn = nn - 1
  n = n - 1
  GO TO 10
END IF

! Allocate various arrays

IF (ALLOCATED(p)) DEALLOCATE(p, qp, k, qk, svk, temp, pt)
ALLOCATE( p(nn), qp(nn), k(nn), qk(nn), svk(nn), temp(nn), pt(nn) )

! MAKE A COPY OF THE COEFFICIENTS
p(1:nn) = op(1:nn)

! START THE ALGORITHM FOR ONE ZERO
30 IF (n <= 2) THEN
  IF (n < 1) RETURN

! CALCULATE THE FINAL ZERO OR PAIR OF ZEROS
  IF (n /= 2) THEN
    zeror(degree) = -p(2) / p(1)
    zeroi(degree) = 0.0_r10
    RETURN
  END IF
  CALL quadratic(p(1), p(2), p(3), zeror(degree-1), zeroi(degree-1),  &
                 zeror(degree), zeroi(degree))
  RETURN
END IF

! FIND LARGEST AND SMALLEST MODULI OF COEFFICIENTS.
MAX = 0.
MIN = infin
DO  i = 1, nn
  x = ABS(p(i))
  IF (x > MAX) MAX = x
  IF (x /= 0. .AND. x < MIN) MIN = x
END DO

! SCALE IF THERE ARE LARGE OR VERY SMALL COEFFICIENTS COMPUTES 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.

sc = lo / MIN
IF (sc <= 1.0) THEN
  IF (MAX < 10.) GO TO 60
  IF (sc == 0.) sc = smalno
ELSE
  IF (infin/sc < MAX) GO TO 60
END IF
l = LOG(sc) / LOG(base) + .5
factor = (base) ** l
IF (factor%hi /= 1._r10) THEN
  DO i = 1, nn
    p(i) = factor * p(i)
  END DO
END IF

! COMPUTE LOWER BOUND ON MODULI OF ZEROS.
60 DO i = 1, nn
  pt(i) = ABS(p(i))
END DO
pt(nn) = -pt(nn)

! COMPUTE UPPER ESTIMATE OF BOUND
x = EXP( (LOG(-pt(nn)) - LOG(pt(1)) ) / n )
IF (pt(n) /= 0.) THEN

! IF NEWTON STEP AT THE ORIGIN IS BETTER, USE IT.
  xm = -pt(nn) / pt(n)
  IF (xm < x) x = xm
END IF

! CHOP THE INTERVAL (0,X) UNTIL FF <= 0
80 xm = x * .1
ff = pt(1)
DO  i = 2, nn
  ff = ff * xm + pt(i)
END DO
IF (ff > 0.) THEN
  x = xm
  GO TO 80
END IF
dx = x

! DO NEWTON ITERATION UNTIL X CONVERGES TO TWO DECIMAL PLACES
100 IF (ABS(dx/x) > .005) THEN
  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
  GO TO 100
END IF
bnd = x

! COMPUTE THE DERIVATIVE AS THE INTIAL K POLYNOMIAL
! AND DO 5 STEPS WITH NO SHIFT
nm1 = n - 1
DO  i = 2, n
  k(i) = (nn-i) * p(i) / REAL(n)
END DO
k(1) = p(1)
aa = p(nn)
bb = p(n)
zerok = k(n)%hi == 0._r10
DO  jj = 1, 5
  cc = k(n)
  IF (.NOT.zerok) THEN
! USE SCALED FORM OF RECURRENCE IF VALUE OF K AT 0 IS NONZERO
    t = -aa / cc
    DO  i = 1, nm1
      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.
  ELSE
! USE UNSCALED FORM OF RECURRENCE
    DO  i = 1, nm1
      j = nn - i
      k(j) = k(j-1)
    END DO
    k(1) = 0._r10
    zerok = k(n)%hi == 0._r10
  END IF
END DO

! SAVE K FOR RESTARTS WITH NEW SHIFTS
temp(1:n) = k(1:n)

! LOOP TO SELECT THE QUADRATIC  CORRESPONDING TO EACH NEW SHIFT
DO  cnt = 1, 20
! QUADRATIC CORRESPONDS TO A DOUBLE SHIFT TO A NON-REAL POINT AND ITS COMPLEX
! CONJUGATE.  THE POINT HAS 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
  u = -2.0_r10 * sr
  v = bnd

! SECOND STAGE CALCULATION, FIXED QUADRATIC
  CALL fxshfr(20*cnt,nz)
  IF (nz /= 0) THEN

! THE SECOND STAGE JUMPS DIRECTLY TO ONE OF THE THIRD STAGE ITERATIONS AND
! RETURNS HERE IF SUCCESSFUL.
! DEFLATE THE POLYNOMIAL, STORE THE ZERO OR ZEROS AND RETURN TO THE MAIN
! ALGORITHM.
    j = degree - n + 1
    zeror(j) = szr
    zeroi(j) = szi
    nn = nn - nz
    n = nn - 1
    p(1:nn) = qp(1:nn)
    IF (nz == 1) GO TO 30
    zeror(j+1) = lzr
    zeroi(j+1) = lzi
    GO TO 30
  END IF

! IF THE ITERATION IS UNSUCCESSFUL ANOTHER QUADRATIC
! IS CHOSEN AFTER RESTORING K
  k(1:n) = temp(1:n)
END DO

! RETURN WITH FAILURE IF NO CONVERGENCE WITH 20 SHIFTS
fail = .true.
degree = degree - n
RETURN
END SUBROUTINE rpoly



SUBROUTINE fxshfr(l2, nz)
! COMPUTES UP TO  L2  FIXED SHIFT K-POLYNOMIALS, TESTING FOR CONVERGENCE IN
! THE LINEAR OR QUADRATIC CASE.  INITIATES ONE OF THE VARIABLE SHIFT
! ITERATIONS AND RETURNS WITH THE NUMBER OF ZEROS FOUND.
! L2 - LIMIT OF FIXED SHIFT STEPS
! NZ - NUMBER OF ZEROS FOUND

INTEGER, INTENT(IN)   :: l2
INTEGER, INTENT(OUT)  :: nz

TYPE (quad)  :: svu, svv, ui, vi, s
REAL (r10)   :: betas, betav, oss, ovv, ss, vv, ts, tv, ots, otv, tvv, tss
INTEGER      :: TYPE, j, iflag
LOGICAL      :: vpass, spass, vtry, stry

nz = 0
betav = .25
betas = .25
oss = sr
ovv = v

! EVALUATE POLYNOMIAL BY SYNTHETIC DIVISION
CALL quadsd(nn, u, v, p, qp, a, b)
CALL calcsc(TYPE)
DO  j = 1, l2
! CALCULATE NEXT K POLYNOMIAL AND ESTIMATE V
  CALL nextk(TYPE)
  CALL calcsc(TYPE)
  CALL newest(TYPE,ui,vi)
  vv = vi

! ESTIMATE S
  ss = 0.
  IF (k(n)%hi /= 0._r10) ss = -p(nn) / k(n)
  tv = 1.
  ts = 1.
  IF (j /= 1 .AND. TYPE /= 3) THEN

! COMPUTE RELATIVE MEASURES OF CONVERGENCE OF S AND V SEQUENCES
    IF (vv /= 0.) tv = ABS((vv-ovv)/vv)
    IF (ss /= 0.) ts = ABS((ss-oss)/ss)

! IF DECREASING, MULTIPLY TWO MOST RECENT CONVERGENCE MEASURES
    tvv = 1.
    IF (tv < otv) tvv = tv * otv
    tss = 1.
    IF (ts < ots) tss = ts * ots

! COMPARE WITH CONVERGENCE CRITERIA
    vpass = tvv < betav
    spass = tss < betas
    IF ((spass.OR.vpass)) THEN

! AT LEAST ONE SEQUENCE HAS PASSED THE CONVERGENCE TEST.
! STORE VARIABLES BEFORE ITERATING
      svu = u
      svv = v
      svk(1:n) = k(1:n)
      s = ss

! CHOOSE ITERATION ACCORDING TO THE FASTEST CONVERGING SEQUENCE
      vtry = .false.
      stry = .false.
      IF (spass .AND. ((.NOT.vpass) .OR. tss < tvv)) GO TO 40
      20 CALL quadit(ui, vi, nz)
      IF (nz > 0) RETURN

! QUADRATIC ITERATION HAS FAILED.  FLAG THAT IT HAS
! BEEN TRIED AND DECREASE THE CONVERGENCE CRITERION.
      vtry = .true.
      betav = betav * .25

! TRY LINEAR ITERATION IF IT HAS NOT BEEN TRIED AND
! THE S SEQUENCE IS CONVERGING
      IF (stry .OR. (.NOT.spass)) GO TO 50
      k(1:n) = svk(1:n)
      40 CALL realit(s, nz, iflag)
      IF (nz > 0) RETURN

! LINEAR ITERATION HAS FAILED.  FLAG THAT IT HAS BEEN
! TRIED AND DECREASE THE CONVERGENCE CRITERION
      stry = .true.
      betas = betas * .25
      IF (iflag /= 0) THEN

! IF LINEAR ITERATION SIGNALS AN ALMOST DOUBLE REAL
! ZERO ATTEMPT QUADRATIC INTERATION
        ui = -(s+s)
        vi = s * s
        GO TO 20
      END IF

! RESTORE VARIABLES
      50 u = svu
      v = svv
      k(1:n) = svk(1:n)

! TRY QUADRATIC ITERATION IF IT HAS NOT BEEN TRIED
! AND THE V SEQUENCE IS CONVERGING
      IF (vpass .AND. (.NOT.vtry)) GO TO 20

! RECOMPUTE QP AND SCALAR VALUES TO CONTINUE THE SECOND STAGE
      CALL quadsd(nn, u, v, p, qp, a, b)
      CALL calcsc(TYPE)
    END IF
  END IF
  ovv = vv
  oss = ss
  otv = tv
  ots = ts
END DO
RETURN
END SUBROUTINE fxshfr


SUBROUTINE quadit(uu, vv, nz)
! VARIABLE-SHIFT K-POLYNOMIAL ITERATION FOR A QUADRATIC FACTOR, CONVERGES
! ONLY IF THE ZEROS ARE EQUIMODULAR OR NEARLY SO.
! UU,VV - COEFFICIENTS OF STARTING QUADRATIC
! NZ - NUMBER OF ZERO FOUND

TYPE (quad), INTENT(IN)  :: uu
TYPE (quad), INTENT(IN)  :: vv
INTEGER, INTENT(OUT)     :: nz

TYPE (quad)  :: ui, vi
REAL (r10)   :: mp, omp, ee, relstp, t, zm
INTEGER      :: TYPE, i, j
LOGICAL      :: tried

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

! MAIN LOOP
10 CALL quadratic(quad(1._r10,0.0_r10), u, v, szr, szi, lzr, lzi)

! RETURN IF ROOTS OF THE QUADRATIC ARE REAL AND NOT
! CLOSE TO MULTIPLE OR NEARLY EQUAL AND  OF OPPOSITE SIGN
IF (ABS(ABS(szr)-ABS(lzr)) > .01_r10*ABS(lzr)) RETURN

! EVALUATE POLYNOMIAL BY QUADRATIC SYNTHETIC DIVISION
CALL quadsd(nn, u, v, p, qp, a, b)
mp = ABS(a-szr*b) + ABS(szi*b)

! COMPUTE A RIGOROUS  BOUND ON THE ROUNDING ERROR IN EVALUTING P
zm = SQRT(ABS(v))
ee = 2. * 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.*mre + 4.*are) * ee - (5.*mre + 2.*are) * (ABS(a) + t) +  &
     ABS(b)*zm + 2. * are * ABS(t)

! ITERATION HAS CONVERGED SUFFICIENTLY IF THE
! POLYNOMIAL VALUE IS LESS THAN 20 TIMES THIS BOUND
IF (mp <= 20.*ee) THEN
  nz = 2
  RETURN
END IF
j = j + 1

! STOP ITERATION AFTER 20 STEPS
IF (j > 20) RETURN
IF (j >= 2) THEN
  IF (.NOT.(relstp > .01 .OR. mp < omp .OR. tried)) THEN

! A CLUSTER APPEARS TO BE STALLING THE CONVERGENCE.
! FIVE FIXED SHIFT STEPS ARE TAKEN WITH A U,V CLOSE TO THE CLUSTER
    IF (relstp < eta) 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(TYPE)
      CALL nextk(TYPE)
    END DO
    tried = .true.
    j = 0
  END IF
END IF
omp = mp

! CALCULATE NEXT K POLYNOMIAL AND NEW U AND V
CALL calcsc(TYPE)
CALL nextk(TYPE)
CALL calcsc(TYPE)
CALL newest(TYPE, ui, vi)

! IF VI IS ZERO THE ITERATION IS NOT CONVERGING
IF (vi%hi == 0._r10) RETURN
relstp = ABS((vi-v)/vi)
u = ui
v = vi
GO TO 10
END SUBROUTINE quadit


SUBROUTINE realit(sss, nz, iflag)
! VARIABLE-SHIFT H POLYNOMIAL ITERATION FOR A REAL ZERO.
! SSS   - STARTING ITERATE
! NZ    - NUMBER OF ZERO FOUND
! IFLAG - FLAG TO INDICATE A PAIR OF ZEROS NEAR REAL AXIS.

TYPE (quad), INTENT(IN OUT)  :: sss
INTEGER, INTENT(OUT)         :: nz
INTEGER, INTENT(OUT)         :: iflag

TYPE (quad)  :: pv, kv, t, s
REAL (r10)   :: ms, mp, omp, ee
INTEGER      :: i, j

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

! MAIN LOOP
10 pv = p(1)

! EVALUATE P AT S
qp(1) = pv
DO  i = 2, nn
  pv = pv * s + p(i)
  qp(i) = pv
END DO
mp = ABS(pv)

! COMPUTE A RIGOROUS BOUND ON THE ERROR IN EVALUATING P
ms = ABS(s)
ee = (mre/(are+mre)) * ABS(qp(1))
DO  i = 2, nn
  ee = ee * ms + ABS(qp(i))
END DO

! ITERATION HAS CONVERGED SUFFICIENTLY IF THE
! POLYNOMIAL VALUE IS LESS THAN 20 TIMES THIS BOUND
IF (mp <= 20.*((are+mre)*ee - mre*mp)) THEN
  nz = 1
  szr = s
  szi = 0._r10
  RETURN
END IF
j = j + 1

! STOP ITERATION AFTER 10 STEPS
IF (j > 10) RETURN
IF (j >= 2) THEN
  IF (ABS(t) <= .001*ABS(s-t) .AND. mp > omp) THEN

! A CLUSTER OF ZEROS NEAR THE REAL AXIS HAS BEEN ENCOUNTERED RETURN WITH
! IFLAG SET TO INITIATE A QUADRATIC ITERATION
    iflag = 1
    sss = s
    RETURN
  END IF
END IF

! RETURN IF THE POLYNOMIAL VALUE HAS INCREASED SIGNIFICANTLY
omp = mp

! COMPUTE T, THE NEXT POLYNOMIAL, AND THE NEW ITERATE
kv = k(1)
qk(1) = kv
DO  i = 2, n
  kv = kv * s + k(i)
  qk(i) = kv
END DO
IF (ABS(kv) > ABS(k(n))*10.*eta) THEN

! USE THE SCALED FORM OF THE RECURRENCE IF THE VALUE OF K AT S IS NONZERO
  t = -pv / kv
  k(1) = qp(1)
  DO  i = 2, n
    k(i) = t * qk(i-1) + qp(i)
  END DO
ELSE
! USE UNSCALED FORM
  k(1) = 0.0_r10
  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
t = 0._r10
IF (ABS(kv) > ABS(k(n))*10.*eta) t = -pv / kv
s = s + t
GO TO 10
END SUBROUTINE realit


SUBROUTINE calcsc(TYPE)
! THIS ROUTINE CALCULATES SCALAR QUANTITIES USED TO COMPUTE THE NEXT K
! POLYNOMIAL AND NEW ESTIMATES OF THE QUADRATIC COEFFICIENTS.
! TYPE - INTEGER VARIABLE SET HERE INDICATING HOW THE CALCULATIONS ARE
! NORMALIZED TO AVOID OVERFLOW

INTEGER, INTENT(OUT)  :: TYPE

! SYNTHETIC DIVISION OF K BY THE QUADRATIC 1,U,V
CALL quadsd(n, u, v, k, qk, c, d)
IF (ABS(c) <= ABS(k(n))*100.*eta) THEN
  IF (ABS(d) <= ABS(k(n-1))*100.*eta) THEN
    TYPE = 3

! TYPE=3 INDICATES THE QUADRATIC IS ALMOST A FACTOR OF K
    RETURN
  END IF
END IF
IF (ABS(d) >= ABS(c)) THEN
  TYPE = 2

! TYPE=2 INDICATES THAT ALL FORMULAS ARE DIVIDED BY D
  ee = a / d
  f = c / d
  g = u * b
  h = v * b
  a3 = (a+g) * ee + h * (b/d)
  a1 = b * f - a
  a7 = (f+u) * a + h
  RETURN
END IF
TYPE = 1

! TYPE=1 INDICATES THAT ALL FORMULAS ARE DIVIDED BY C
ee = a / c
f = d / c
g = u * ee
h = v * b
a3 = a * ee + (h/c+g) * b
a1 = b - a * (d/c)
a7 = a + g * d + h * f
RETURN
END SUBROUTINE calcsc


SUBROUTINE nextk(TYPE)
! COMPUTES THE NEXT K POLYNOMIALS USING SCALARS COMPUTED IN CALCSC

INTEGER, INTENT(IN)  :: TYPE

TYPE (quad) :: temp
INTEGER     :: i

IF (TYPE /= 3) THEN
  temp = a
  IF (TYPE == 1) temp = b
  IF (ABS(a1) <= ABS(temp)*eta*10.) THEN

! IF A1 IS NEARLY ZERO THEN USE A SPECIAL FORM OF THE RECURRENCE
    k(1) = 0._r10
    k(2) = -a7 * qp(1)
    DO  i = 3, n
      k(i) = a3 * qk(i-2) - a7 * qp(i-1)
    END DO
    RETURN
  END IF

! USE SCALED FORM OF THE RECURRENCE
  a7 = a7 / a1
  a3 = a3 / a1
  k(1) = qp(1)
  k(2) = qp(2) - a7 * qp(1)
  DO  i = 3, n
    k(i) = a3 * qk(i-2) - a7 * qp(i-1) + qp(i)
  END DO
  RETURN
END IF

! USE UNSCALED FORM OF THE RECURRENCE IF TYPE IS 3
k(1) = 0._r10
k(2) = 0._r10
DO  i = 3, n
  k(i) = qk(i-2)
END DO
RETURN
END SUBROUTINE nextk


SUBROUTINE newest(TYPE, uu, vv)
! COMPUTE NEW ESTIMATES OF THE QUADRATIC COEFFICIENTS
! USING THE SCALARS COMPUTED IN CALCSC.

INTEGER, INTENT(IN)       :: TYPE
TYPE (quad), INTENT(OUT)  :: uu
TYPE (quad), INTENT(OUT)  :: vv

TYPE (quad) :: a4, a5, b1, b2, c1, c2, c3, c4, temp

! USE FORMULAS APPROPRIATE TO SETTING OF TYPE.
IF (TYPE /= 3) THEN
  IF (TYPE /= 2) 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

! EVALUATE NEW QUADRATIC COEFFICIENTS.
  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%hi /= 0._r10) THEN
    uu = u - (u*(c3 + c2) + v*(b1*a1 + b2*a7)) / temp
    vv = v * (1. + c4/temp)
    RETURN
  END IF
END IF

! IF TYPE=3 THE QUADRATIC IS ZEROED
uu = 0._r10
vv = 0._r10
RETURN
END SUBROUTINE newest


SUBROUTINE quadsd(nn, u, v, p, q, a, b)
! DIVIDES P BY THE QUADRATIC  1,U,V  PLACING THE
! QUOTIENT IN Q AND THE REMAINDER IN A,B

INTEGER, INTENT(IN)       :: nn
TYPE (quad), INTENT(IN)   :: u
TYPE (quad), INTENT(IN)   :: v
TYPE (quad), INTENT(IN)   :: p(:)
TYPE (quad), INTENT(OUT)  :: q(:)
TYPE (quad), INTENT(OUT)  :: a
TYPE (quad), INTENT(OUT)  :: b

TYPE (quad) :: c
INTEGER     :: i

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 quadsd


SUBROUTINE quadratic(a, b1, c, sr, si, lr, li)
! CALCULATE THE ZEROS OF THE QUADRATIC A*Z**2 + B1*Z + C.
! THE QUADRATIC FORMULA, MODIFIED TO AVOID OVERFLOW, IS USED TO FIND THE
! LARGER ZERO IF THE ZEROS ARE REAL AND BOTH ZEROS ARE COMPLEX.
! THE SMALLER REAL ZERO IS FOUND DIRECTLY FROM THE PRODUCT OF THE ZEROS C/A.

TYPE (quad), INTENT(IN)   :: a
TYPE (quad), INTENT(IN)   :: b1
TYPE (quad), INTENT(IN)   :: c
TYPE (quad), INTENT(OUT)  :: sr
TYPE (quad), INTENT(OUT)  :: si
TYPE (quad), INTENT(OUT)  :: lr
TYPE (quad), INTENT(OUT)  :: li

TYPE (quad) :: b, d, ee

IF (a%hi /= 0._r10) GO TO 20
sr = 0._r10
IF (b1%hi /= 0._r10) sr = -c / b1
lr = 0._r10
10 si = 0._r10
li = 0._r10
RETURN

20 IF (c%hi == 0._r10) THEN
  sr = 0._r10
  lr = -b1 / a
  GO TO 10
END IF

! COMPUTE DISCRIMINANT AVOIDING OVERFLOW
b = b1 / 2._r10
IF (ABS(b) >= ABS(c)) THEN
  ee = 1._r10 - (a/b) * (c/b)
  d = SQRT(ABS(ee)) * ABS(b)
ELSE
  ee = a
  IF (c%hi < 0._r10) ee = -a
  ee = b * (b/ABS(c)) - ee
  d = SQRT(ABS(ee)) * SQRT(ABS(c))
END IF

IF (ee%hi >= 0._r10) THEN

! REAL ZEROS
  IF (b%hi >= 0._r10) d = -d
  lr = (-b+d) / a
  sr = 0._r10
  IF (lr%hi /= 0._r10) sr = (c/lr) / a
  GO TO 10
END IF

! COMPLEX CONJUGATE ZEROS
sr = -b / a
lr = sr
si = ABS(d/a)
li = -si
RETURN
END SUBROUTINE quadratic

END MODULE Solve_Real_Poly



PROGRAM test_rpoly
USE quadruple_precision
USE Solve_Real_Poly
IMPLICIT NONE

TYPE (quad)  :: p(50), zr(50), zi(50)
INTEGER      :: degree, i
LOGICAL      :: fail

WRITE(*, 5000)
degree = 10
p(1) = 1.
p(2) = -55.
p(3) = 1320.
p(4) = -18150.
p(5) = 157773.
p(6) = -902055.
p(7) = 3416930.
p(8) = -8409500.
p(9) = 12753576.
p(10) = -10628640.
p(11) = 3628800.

CALL rpoly(p, degree, zr, zi, fail)
IF (fail) THEN
  WRITE(*, *) ' ** Failure by RPOLY **'
ELSE
  WRITE(*, '(a/ (2g26.18))') ' Real part              Imaginary part',  &
                             (zr(i)%hi, zi(i)%hi, i=1,degree)
END IF

STOP

5000 FORMAT (' EXAMPLE 1. POLYNOMIAL WITH ZEROS 1,2,...,10.')

END PROGRAM test_rpoly