MODULE Solve_Complex_Poly ! ALGORITHM 419 COLLECTED ALGORITHMS FROM ACM. ! Code converted using TO_F90 by Alan Miller ! Date: 2000-01-08 Time: 16:02:44 ! Latest revision - 10 February 2000 ! amiller @ bigpond.net.au ! ALGORITHM APPEARED IN COMM. ACM, VOL. 15, NO. 02, P. 097. USE quadruple_precision IMPLICIT NONE ! COMMON AREA ! COMMON /global/ pr, pi, hr, hi, qpr, qpi, qhr, qhi, shr, shi, sr, & ! si, tr, ti, pvr, pvi, are, mre, eta, infin, nn ! TYPE (quad) :: sr, si, tr, ti, pvr, pvi, are, mre, eta, infin, & ! pr(50), pi(50), hr(50), hi(50), qpr(50), qpi(50), qhr(50), & ! qhi(50), shr(50), shi(50) ! N.B. pi has been changed to ci as the constant pi = 3.14159... ! is available from the quad.precision module. TYPE (quad), ALLOCATABLE, SAVE :: pr(:), ci(:), hr(:), hi(:), qpr(:), qpi(:), & qhr(:), qhi(:), shr(:), shi(:) TYPE (quad), SAVE :: sr, si, tr, ti, pvr, pvi, are, mre, eta, & infin INTEGER, SAVE :: nn PRIVATE PUBLIC :: cpoly CONTAINS SUBROUTINE cpoly(opr, opi, degree, zeror, zeroi, fail) ! FINDS THE ZEROS OF A COMPLEX POLYNOMIAL. ! OPR, OPI - TYPE (quad) VECTORS OF REAL AND IMAGINARY PARTS OF THE ! 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 LEADING COEFFICIENT ! IS ZERO OR IF CPOLY HAS FOUND FEWER THAN DEGREE ZEROS. ! THE PROGRAM HAS BEEN WRITTEN TO REDUCE THE CHANCE OF OVERFLOW OCCURRING. ! IF IT DOES OCCUR, THERE IS STILL A POSSIBILITY THAT THE ZEROFINDER WILL WORK ! PROVIDED THE OVERFLOWED QUANTITY IS REPLACED BY A LARGE NUMBER. TYPE (quad), INTENT(IN) :: opr(:) TYPE (quad), INTENT(IN) :: opi(:) INTEGER, INTENT(IN) :: degree TYPE (quad), INTENT(OUT) :: zeror(:) TYPE (quad), INTENT(OUT) :: zeroi(:) LOGICAL, INTENT(OUT) :: fail TYPE (quad) :: xx, yy, cosr, sinr, smalno, base, xxx, zr, zi, bnd LOGICAL :: conv INTEGER :: cnt1, cnt2, i, idnn2 ! INITIALIZATION OF CONSTANTS CALL mcon(eta, infin, smalno, base) are = eta mre = 2.0_dp * SQRT(2.0_dp) * eta xx = .70710678 yy = -xx cosr = -.060756474 sinr = .99756405 fail = .false. nn = degree + 1 ! ALGORITHM FAILS IF THE LEADING COEFFICIENT IS ZERO. IF (opr(1)%hi == 0.0_dp .AND. opi(1)%hi == 0.0_dp) THEN fail = .true. RETURN END IF ! Allocate arrays IF (ALLOCATED(pr)) DEALLOCATE(pr, ci, hr, hi, qpr, qpi, qhr, qhi, shr, shi) ALLOCATE(pr(nn), ci(nn), hr(nn), hi(nn), qpr(nn), qpi(nn), qhr(nn), & qhi(nn), shr(nn), shi(nn)) ! REMOVE THE ZEROS AT THE ORIGIN IF ANY. 10 IF (opr(nn)%hi == 0.0_dp .AND. opi(nn)%hi == 0.0_dp) THEN idnn2 = degree - nn + 2 zeror(idnn2) = 0.0_dp zeroi(idnn2) = 0.0_dp nn = nn - 1 GO TO 10 END IF ! MAKE A COPY OF THE COEFFICIENTS. DO i = 1, nn pr(i) = opr(i) ci(i) = opi(i) shr(i) = cmod(pr(i), ci(i)) END DO ! SCALE THE POLYNOMIAL. bnd = scalepoly(nn, shr, eta, infin, smalno, base) IF (bnd%hi /= 1.0_dp) THEN DO i = 1, nn pr(i) = bnd * pr(i) ci(i) = bnd * ci(i) END DO END IF ! START THE ALGORITHM FOR ONE ZERO . 40 IF (nn <= 2) THEN ! CALCULATE THE FINAL ZERO AND RETURN. CALL cdivid(-pr(2), -ci(2), pr(1), ci(1), zeror(degree), zeroi(degree)) RETURN END IF ! CALCULATE BND, A LOWER BOUND ON THE MODULUS OF THE ZEROS. DO i = 1, nn shr(i) = cmod(pr(i), ci(i)) END DO CALL cauchy(nn, shr, shi, bnd) ! 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) ci(i) = qpi(i) END DO GO TO 40 END IF ! IF THE ITERATION IS UNSUCCESSFUL ANOTHER SHIFT IS CHOSEN. END DO ! IF 9 SHIFTS FAIL, THE OUTER LOOP IS REPEATED WITH ANOTHER SEQUENCE OF SHIFTS. END DO ! THE ZEROFINDER HAS FAILED ON TWO MAJOR PASSES. ! RETURN EMPTY HANDED. fail = .true. RETURN END SUBROUTINE cpoly SUBROUTINE noshft(l1) ! COMPUTES THE DERIVATIVE POLYNOMIAL AS THE INITIAL H ! POLYNOMIAL AND COMPUTES L1 NO-SHIFT H POLYNOMIALS. INTEGER, INTENT(IN) :: l1 TYPE (quad) :: xni, t1, t2 INTEGER :: i, j, jj, n, nm1 n = nn - 1 nm1 = n - 1 DO i = 1, n xni = nn - i hr(i) = xni * pr(i) / n hi(i) = xni * ci(i) / n END DO DO jj = 1, l1 IF (cmod(hr(n), hi(n)) > eta*10.0_dp*cmod(pr(n), ci(n))) THEN CALL cdivid(-pr(nn), -ci(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 + ci(j) END DO hr(1) = pr(1) hi(1) = ci(1) ELSE ! 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) END DO hr(1) = 0.0_dp hi(1) = 0.0_dp END IF END DO RETURN 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. ! L2 - LIMIT OF FIXED SHIFT STEPS ! ZR,ZI - APPROXIMATE ZERO IF CONV IS .TRUE. ! CONV - LOGICAL INDICATING CONVERGENCE OF STAGE 3 ITERATION INTEGER, INTENT(IN) :: l2 TYPE (quad), INTENT(OUT) :: zr TYPE (quad), INTENT(OUT) :: zi LOGICAL, INTENT(OUT) :: conv TYPE (quad) :: otr, oti, svsr, svsi LOGICAL :: test, pasd, bool INTEGER :: i, j, n n = nn - 1 ! EVALUATE P AT S. CALL polyev(nn, sr, si, pr, ci, 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) < .5_dp*cmod(zr,zi)) THEN IF (pasd) THEN ! 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) END DO 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) END DO sr = svsr si = svsi CALL polyev(nn, sr, si, pr, ci, qpr, qpi, pvr, pvi) CALL calct(bool) CYCLE END IF pasd = .true. ELSE pasd = .false. END IF END IF END DO ! ATTEMPT AN ITERATION WITH FINAL H POLYNOMIAL FROM SECOND STAGE. CALL vrshft(10, zr, zi, conv) RETURN END SUBROUTINE fxshft SUBROUTINE vrshft(l3, zr, zi, conv) ! CARRIES OUT THE THIRD STAGE ITERATION. ! L3 - LIMIT OF STEPS IN STAGE 3. ! ZR,ZI - ON ENTRY CONTAINS THE INITIAL ITERATE, IF THE ! ITERATION CONVERGES IT CONTAINS THE FINAL ITERATE ON EXIT. ! CONV - .TRUE. IF ITERATION CONVERGES INTEGER, INTENT(IN) :: l3 TYPE (quad), INTENT(IN OUT) :: zr TYPE (quad), INTENT(IN OUT) :: zi LOGICAL, INTENT(OUT) :: conv TYPE (quad) :: 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, ci, qpr, qpi, pvr, pvi) mp = cmod(pvr, pvi) ms = cmod(sr,si) IF (mp <= 20.0_dp*errev(nn, qpr, qpi, ms, mp, are, mre)) THEN ! 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 END IF IF (i /= 1) THEN IF (.NOT.(b .OR. mp < omp .OR. relstp >= quad(.05_dp,0.0_dp))) THEN ! 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_dp+r1) - si * r1 si = sr * r1 + si * (1.0_dp+r1) sr = r2 CALL polyev(nn, sr, si, pr, ci, qpr, qpi, pvr, pvi) DO j = 1, 5 CALL calct(bool) CALL nexth(bool) END DO omp = infin GO TO 20 END IF ! EXIT IF POLYNOMIAL VALUE INCREASES SIGNIFICANTLY. IF (mp*.1_dp > omp) RETURN END IF omp = mp ! CALCULATE NEXT ITERATE. 20 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 END IF END DO RETURN END SUBROUTINE vrshft SUBROUTINE calct(bool) ! COMPUTES T = -P(S)/H(S). ! BOOL - LOGICAL, SET TRUE IF H(S) IS ESSENTIALLY ZERO. LOGICAL, INTENT(OUT) :: bool TYPE (quad) :: 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_dp * cmod(hr(n), hi(n)) IF (.NOT.bool) THEN CALL cdivid(-pvr, -pvi, hvr, hvi, tr, ti) RETURN END IF tr = 0.0_dp ti = 0.0_dp RETURN END SUBROUTINE calct SUBROUTINE nexth(bool) ! CALCULATES THE NEXT SHIFTED H POLYNOMIAL. ! BOOL - LOGICAL, IF .TRUE. H(S) IS ESSENTIALLY ZERO LOGICAL, INTENT(IN) :: bool TYPE (quad) :: t1, t2 INTEGER :: j, n n = nn - 1 IF (.NOT.bool) THEN 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) END DO hr(1) = qpr(1) hi(1) = qpi(1) RETURN END IF ! IF H(S) IS ZERO REPLACE H WITH QH. DO j = 2, n hr(j) = qhr(j-1) hi(j) = qhi(j-1) END DO hr(1) = 0.0_dp hi(1) = 0.0_dp RETURN END SUBROUTINE nexth SUBROUTINE polyev(nn, sr, si, pr, ci, 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. INTEGER, INTENT(IN) :: nn TYPE (quad), INTENT(IN) :: sr TYPE (quad), INTENT(IN) :: si TYPE (quad), INTENT(IN) :: pr(:) TYPE (quad), INTENT(IN) :: ci(:) TYPE (quad), INTENT(OUT) :: qr(:) TYPE (quad), INTENT(OUT) :: qi(:) TYPE (quad), INTENT(OUT) :: pvr TYPE (quad), INTENT(OUT) :: pvi TYPE (quad) :: t INTEGER :: i qr(1) = pr(1) qi(1) = ci(1) pvr = qr(1) pvi = qi(1) DO i = 2, nn t = pvr * sr - pvi * si + pr(i) pvi = pvr * si + pvi * sr + ci(i) pvr = t qr(i) = pvr qi(i) = pvi END DO RETURN END SUBROUTINE polyev FUNCTION errev(nn, qr, qi, ms, mp, are, mre) RESULT(fn_val) ! BOUNDS THE ERROR IN EVALUATING THE POLYNOMIAL BY THE HORNER RECURRENCE. ! QR,QI - THE PARTIAL SUMS ! MS -MODULUS OF THE POINT ! MP -MODULUS OF POLYNOMIAL VALUE ! ARE, MRE -ERROR BOUNDS ON COMPLEX ADDITION AND MULTIPLICATION INTEGER, INTENT(IN) :: nn TYPE (quad), INTENT(IN) :: qr(:) TYPE (quad), INTENT(IN) :: qi(:) TYPE (quad), INTENT(IN) :: ms TYPE (quad), INTENT(IN) :: mp TYPE (quad), INTENT(IN) :: are TYPE (quad), INTENT(IN) :: mre TYPE (quad) :: fn_val TYPE (quad) :: e INTEGER :: i e = cmod(qr(1), qi(1)) * mre / (are+mre) DO i = 1, nn e = e * ms + cmod(qr(i), qi(i)) END DO fn_val = e * (are+mre) - mp * mre RETURN END FUNCTION errev SUBROUTINE cauchy(nn, pt, q, fn_val) ! CAUCHY COMPUTES A LOWER BOUND ON THE MODULI OF THE ZEROS OF A ! POLYNOMIAL - PT IS THE MODULUS OF THE COEFFICIENTS. INTEGER, INTENT(IN) :: nn TYPE (quad), INTENT(OUT) :: pt(:), q(:), fn_val TYPE (quad) :: x, xm, f, dx, df INTEGER :: i, n pt(nn) = -pt(nn) ! COMPUTE UPPER ESTIMATE OF BOUND. n = nn - 1 x = EXP((LOG(-pt(nn)) - LOG(pt(1))) / n) IF (pt(n)%hi /= 0.0_dp) 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) UNITL F LE 0. 10 xm = x * .1_dp f = pt(1) DO i = 2, nn f = f * xm + pt(i) END DO IF (f%hi > 0.0_dp) THEN x = xm GO TO 10 END IF dx = x ! DO NEWTON ITERATION UNTIL X CONVERGES TO TWO DECIMAL PLACES. 30 IF (ABS(dx/x) > quad(.005_dp,0.0_dp)) THEN q(1) = pt(1) DO i = 2, nn q(i) = q(i-1) * x + pt(i) END DO f = q(nn) df = q(1) DO i = 2, n df = df * x + q(i) END DO dx = f / df x = x - dx GO TO 30 END IF fn_val = x RETURN END SUBROUTINE cauchy FUNCTION scalepoly(nn, pt, eta, infin, smalno, base) RESULT(fn_val) ! 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. ! PT - MODULUS OF COEFFICIENTS OF P ! ETA, INFIN, SMALNO, BASE - CONSTANTS DESCRIBING THE FLOATING POINT ARITHMETIC. INTEGER, INTENT(IN) :: nn TYPE (quad), INTENT(IN) :: pt(:) TYPE (quad), INTENT(IN) :: eta TYPE (quad), INTENT(IN) :: infin TYPE (quad), INTENT(IN) :: smalno TYPE (quad), INTENT(IN) :: base TYPE (quad) :: fn_val TYPE (quad) :: high, low, MAX, MIN, x, sc INTEGER :: i, l ! FIND LARGEST AND SMALLEST MODULI OF COEFFICIENTS. high = SQRT(infin) low = smalno / eta MAX = 0.0_dp MIN = infin DO i = 1, nn x = pt(i) IF (x > MAX) MAX = x IF (x%hi /= 0.0_dp .AND. x < MIN) MIN = x END DO ! SCALE ONLY IF THERE ARE VERY LARGE OR VERY SMALL COMPONENTS. fn_val = 1.0_dp IF (MIN >= low .AND. MAX <= high) RETURN x = low / MIN IF (x%hi <= 1.0_dp) THEN sc = 1.0_dp / (SQRT(MAX)*SQRT(MIN)) ELSE sc = x IF (infin/sc > MAX) sc = 1.0_dp END IF l = LOG(sc) / LOG(base) + .5 fn_val = base ** l RETURN END FUNCTION scalepoly SUBROUTINE cdivid(ar, ai, br, bi, cr, ci) ! COMPLEX DIVISION C = A/B, AVOIDING OVERFLOW. TYPE (quad), INTENT(IN) :: ar TYPE (quad), INTENT(IN) :: ai TYPE (quad), INTENT(IN) :: br TYPE (quad), INTENT(IN) :: bi TYPE (quad), INTENT(OUT) :: cr TYPE (quad), INTENT(OUT) :: ci TYPE (quad) :: r, d, t, infin IF (br%hi == 0.0_dp .AND. bi%hi == 0.0_dp) THEN ! DIVISION BY ZERO, C = INFINITY. CALL mcon(t, infin, t, t) cr = infin ci = infin RETURN END IF IF (ABS(br) < ABS(bi)) THEN r = br / bi d = bi + r * br cr = (ar*r+ai) / d ci = (ai*r-ar) / d RETURN END IF r = bi / br d = br + r * bi cr = (ar+ai*r) / d ci = (ai-ar*r) / d RETURN END SUBROUTINE cdivid FUNCTION cmod(r, i) RESULT(fn_val) ! MODULUS OF A COMPLEX NUMBER AVOIDING OVERFLOW. TYPE (quad), INTENT(IN) :: r TYPE (quad), INTENT(IN) :: i TYPE (quad) :: fn_val TYPE (quad) :: ar, ai ar = ABS(r) ai = ABS(i) IF (ar < ai) THEN fn_val = ai * SQRT(1.0_dp + (ar/ai)**2) RETURN END IF IF (ar > ai) THEN fn_val = ar * SQRT(1.0_dp + (ai/ar)**2) RETURN END IF fn_val = ar * SQRT(2.0_dp) RETURN END FUNCTION cmod SUBROUTINE mcon(eta, infiny, smalno, base) ! MCON PROVIDES MACHINE CONSTANTS USED IN VARIOUS PARTS OF THE PROGRAM. ! THE USER MAY EITHER SET THEM DIRECTLY OR USE THE STATEMENTS BELOW TO ! COMPUTE THEM. 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.0_dp + ETA > 1.0. ! INFINY THE LARGEST FLOATING-POINT NUMBER ! SMALNO THE SMALLEST POSITIVE FLOATING-POINT NUMBER ! BASE THE BASE OF THE FLOATING-POINT NUMBER SYSTEM USED TYPE (quad), INTENT(OUT) :: eta, infiny, smalno, base base = RADIX(0.0_dp) eta = EPSILON(0.0_dp) * 1.0E-15_dp infiny = HUGE(0.0_dp) smalno = TINY(0.0_dp) RETURN END SUBROUTINE mcon END MODULE Solve_Complex_Poly PROGRAM cpolydr ! DRIVER TO TEST CPOLY USE quadruple_precision USE Solve_Complex_Poly IMPLICIT NONE LOGICAL :: fail TYPE (quad) :: p(50), ci(50), zr(50), zi(50) INTEGER :: i WRITE (6,5000) 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 ci(1:11) = quad(0.0_dp,0.0_dp) CALL prtc(11, p, ci) CALL cpoly(p, ci, 10, zr, zi, fail) IF (fail) GO TO 60 CALL prtz(10, zr, zi) 20 WRITE (6,5100) p(1) = 1 p(2) = 0 p(3) = -10001.0001_dp p(4) = 0 ci(1) = 0 ci(2) = -10001.0001_dp ci(3) = 0 ci(4) = 1 CALL prtc(4, p, ci) CALL cpoly(p, ci, 3, zr, zi, fail) IF (fail) GO TO 70 CALL prtz(3,zr,zi) 30 WRITE (6,5200) p(1) = 1.0 p(2) = -1.998046875 p(3) = 0.0 p(4) = .7567065954208374_dp p(5) = -.2002119533717632_dp p(6) = 1.271507365163416D-2 p(7) = 0 p(8) = -1.154642632172909D-5 p(9) = 1.584803612786345D-7 p(10) = -4.652065399568528D-10 p(11) = 0 ci(1) = 0 ci(2) = p(2) ci(3) = 2.658859252929688_dp ci(4) = -7.567065954208374D-1 ci(5) = 0 ci(6) = p(6) ci(7) = -7.820779428584501D-4 ci(8) = -p(8) ci(9) = 0 ci(10) = p(10) ci(11) = 9.094947017729282D-13 CALL prtc(11, p, ci) CALL cpoly(p, ci, 10, zr, zi, fail) IF (fail) GO TO 80 CALL prtz(10, zr, zi) 40 WRITE (6,5300) p(1) = 1 p(2) = -10 p(3) = 3 p(4) = 284 p(5) = -1293 p(6) = 2374 p(7) = -1587 p(8) = -920 p(9) = 2204 p(10) = -1344 p(11) = 288 ci(1) = 0 ci(2) = -10 ci(3) = 100 ci(4) = -334 ci(5) = 200 ci(6) = 1394 ci(7) = -3836 ci(8) = 4334 ci(9) = -2352 ci(10) = 504 ci(11) = 0 CALL prtc(11, p, ci) CALL cpoly(p, ci, 10, zr, zi, fail) IF (fail) GO TO 90 CALL prtz(10, zr, zi) 50 WRITE (6,5400) p(1) = 1 p(2) = 0 p(3) = -264 p(4) = 0 p(5) = 7920 p(6) = 0 p(7) = -59136 p(8) = 0 p(9) = 126720 p(10) = 0 p(11) = -67584 p(12) = 0 p(13) = 4095 ci(1) = 0 ci(2) = -24 ci(3) = 0 ci(4) = 1760 ci(5) = 0 ci(6) = -25344 ci(7) = 0 ci(8) = 101376 ci(9) = 0 ci(10) = -112640 ci(11) = 0 ci(12) = 24576 ci(13) = 0 CALL prtc(13, p, ci) CALL cpoly(p, ci, 12, zr, zi, fail) IF (fail) GO TO 100 CALL prtz(12, zr, zi) STOP 60 WRITE (6,5500) GO TO 20 70 WRITE (6,5500) GO TO 30 80 WRITE (6,5500) GO TO 40 90 WRITE (6,5500) GO TO 50 100 WRITE (6,5500) STOP 5000 FORMAT ('1EXAMPLE 1. POLYNOMIAL WITH ZEROS 1,2,...,10.') 5100 FORMAT ('1EXAMPLE 2. ZEROS ON IMAGINARY AXIS DEGREE 3.') 5200 FORMAT ('1EXAMPLE 3. ZEROS AT 1+I,1/2*(1+I)....1/(2**-9)*(1+I)') 5300 FORMAT ('1EXAMPLE 4. MULTIPLE ZEROS') 5400 FORMAT ('1EXAMPLE 5. 12 ZEROS EVENLY DISTRIBUTED ON A CIRCLE OF', & ' RADIUS 1. CENTERED AT 0+2I.') 5500 FORMAT (//' CPOLY HAS FAILED ON THIS EXAMPLE') CONTAINS SUBROUTINE prtc(n, p, q) INTEGER, INTENT(IN) :: n TYPE (quad), INTENT(IN) :: p(50) TYPE (quad), INTENT(IN) :: q(50) WRITE (6,5000) (p(i)%hi,q(i)%hi,i = 1,n) RETURN 5000 FORMAT (//' COEFFICIENTS'/ 50(2G26.16/)) END SUBROUTINE prtc SUBROUTINE prtz(n, zr, zi) INTEGER, INTENT(IN) :: n TYPE (quad), INTENT(IN) :: zr(50) TYPE (quad), INTENT(IN) :: zi(50) WRITE (6,5000) (zr(i)%hi,zi(i)%hi,i = 1,n) RETURN 5000 FORMAT (//' ZEROS'/50(2G26.16/)) END SUBROUTINE prtz END PROGRAM cpolydr