/* 419.f -- translated by f2c (version of 23 April 1993 18:34:30). You must link the resulting object file with the libraries: -lf2c -lm (in that order) */ #include "f2c.h" /* Common Block Declarations */ struct { doublereal pr[50], pi[50], hr[50], hi[50], qpr[50], qpi[50], qhr[50], qhi[ 50], shr[50], shi[50], sr, si, tr, ti, pvr, pvi, are, mre, eta, infin; integer nn; } global_; #define global_1 global_ /* Table of constant values */ static integer c__9 = 9; static integer c__1 = 1; static integer c__2 = 2; static integer c__3 = 3; static integer c__4 = 4; static integer c__5 = 5; static integer c__11 = 11; static integer c__10 = 10; static integer c__13 = 13; static integer c__12 = 12; static integer c__17 = 17; static integer c__16 = 16; /* ALGORITHM 419 COLLECTED ALGORITHMS FROM ACM. */ /* ALGORITHM APPEARED IN COMM. ACM, VOL. 15, NO. 02, */ /* P. 097. */ /* Main program */ MAIN__(void) { /* Builtin functions */ integer s_wsle(cilist *), do_lio(integer *, integer *, char *, ftnlen), e_wsle(void); /* Local variables */ static logical fail; extern /* Subroutine */ int pchk_(doublereal *, doublereal *, integer *, doublereal *, doublereal *), prtc_(integer *, doublereal *, doublereal *), prtz_(integer *, doublereal *, doublereal *); static integer i; static doublereal p[50]; extern /* Subroutine */ int cpoly_(doublereal *, doublereal *, integer *, doublereal *, doublereal *, logical *); static doublereal pi[50], zi[50], zr[50]; /* Fortran I/O blocks */ static cilist io___1 = { 0, 6, 0, 0, 0 }; static cilist io___7 = { 0, 6, 0, 0, 0 }; static cilist io___8 = { 0, 6, 0, 0, 0 }; static cilist io___9 = { 0, 6, 0, 0, 0 }; static cilist io___10 = { 0, 6, 0, 0, 0 }; static cilist io___11 = { 0, 6, 0, 0, 0 }; static cilist io___12 = { 0, 6, 0, 0, 0 }; static cilist io___13 = { 0, 6, 0, 0, 0 }; static cilist io___14 = { 0, 6, 0, 0, 0 }; static cilist io___15 = { 0, 6, 0, 0, 0 }; static cilist io___16 = { 0, 6, 0, 0, 0 }; static cilist io___18 = { 0, 6, 0, 0, 0 }; static cilist io___19 = { 0, 6, 0, 0, 0 }; static cilist io___20 = { 0, 6, 0, 0, 0 }; static cilist io___21 = { 0, 6, 0, 0, 0 }; static cilist io___22 = { 0, 6, 0, 0, 0 }; static cilist io___23 = { 0, 6, 0, 0, 0 }; static cilist io___24 = { 0, 6, 0, 0, 0 }; static cilist io___25 = { 0, 6, 0, 0, 0 }; static cilist io___26 = { 0, 6, 0, 0, 0 }; static cilist io___27 = { 0, 6, 0, 0, 0 }; static cilist io___28 = { 0, 6, 0, 0, 0 }; static cilist io___29 = { 0, 6, 0, 0, 0 }; static cilist io___30 = { 0, 6, 0, 0, 0 }; /* DRIVER TO TEST CPOLY */ s_wsle(&io___1); do_lio(&c__9, &c__1, "EXAMPLE 0. POLYNOMIAL WITH ZERO 0.5", 36L); e_wsle(); p[0] = 4.; p[1] = -2.; pi[0] = 4.; pi[1] = -2.; prtc_(&c__2, p, pi); cpoly_(p, pi, &c__1, zr, zi, &fail); if (fail) { s_wsle(&io___7); do_lio(&c__9, &c__1, "CPOLY HAS FAILED ON THIS EXAMPLE", 32L); e_wsle(); } prtz_(&c__1, zr, zi); pchk_(p, pi, &c__1, zr, zi); s_wsle(&io___8); do_lio(&c__9, &c__1, "EXAMPLE 0a. POLYNOMIAL WITH ZEROS 0.0 -2.0", 43L); e_wsle(); p[0] = 1.; p[1] = 2.; p[2] = 0.; pi[0] = 1.; pi[1] = 2.; pi[2] = 0.; prtc_(&c__3, p, pi); cpoly_(p, pi, &c__2, zr, zi, &fail); if (fail) { s_wsle(&io___9); do_lio(&c__9, &c__1, "CPOLY HAS FAILED ON THIS EXAMPLE", 32L); e_wsle(); } prtz_(&c__2, zr, zi); pchk_(p, pi, &c__2, zr, zi); s_wsle(&io___10); do_lio(&c__9, &c__1, "EXAMPLE 0b. POLYNOMIAL WITH ZEROS (1+2i) (3-4i)", 48L); e_wsle(); p[0] = 1.; p[1] = -4.; p[2] = 11.; pi[0] = 0.; pi[1] = 2.; pi[2] = 2.; prtc_(&c__3, p, pi); cpoly_(p, pi, &c__2, zr, zi, &fail); if (fail) { s_wsle(&io___11); do_lio(&c__9, &c__1, "CPOLY HAS FAILED ON THIS EXAMPLE", 32L); e_wsle(); } prtz_(&c__2, zr, zi); pchk_(p, pi, &c__2, zr, zi); s_wsle(&io___12); do_lio(&c__9, &c__1, "EXAMPLE 0c. POLYNOMIAL WITH ZEROS 10**15 and", 45L) ; do_lio(&c__9, &c__1, " 10**(-15)", 10L); e_wsle(); p[0] = 1.; p[1] = -1e15; p[2] = 1.; pi[0] = 0.; pi[1] = 0.; pi[2] = 0.; prtc_(&c__3, p, pi); cpoly_(p, pi, &c__2, zr, zi, &fail); if (fail) { s_wsle(&io___13); do_lio(&c__9, &c__1, "CPOLY HAS FAILED ON THIS EXAMPLE", 32L); e_wsle(); } prtz_(&c__2, zr, zi); pchk_(p, pi, &c__2, zr, zi); s_wsle(&io___14); do_lio(&c__9, &c__1, "EXAMPLE 0d. POLYNOMIAL WITH ZEROS", 34L); do_lio(&c__9, &c__1, " (1+2i) (3-4i) (-6+5i)", 22L); e_wsle(); p[0] = 1.; p[1] = 2.; p[2] = -3.; p[3] = 76.; pi[0] = 0.; pi[1] = -3.; pi[2] = 34.; pi[3] = -43.; prtc_(&c__4, p, pi); cpoly_(p, pi, &c__3, zr, zi, &fail); if (fail) { s_wsle(&io___15); do_lio(&c__9, &c__1, "CPOLY HAS FAILED ON THIS EXAMPLE", 32L); e_wsle(); } prtz_(&c__3, zr, zi); pchk_(p, pi, &c__3, zr, zi); s_wsle(&io___16); do_lio(&c__9, &c__1, "EXAMPLE 1. POLYNOMIAL WITH ZEROS 1,2,3,4", 41L); e_wsle(); p[0] = 1.; p[1] = -10.; p[2] = 35.; p[3] = -50.; p[4] = 24.; for (i = 1; i <= 5; ++i) { /* L10: */ pi[i - 1] = 0.; } prtc_(&c__5, p, pi); cpoly_(p, pi, &c__4, zr, zi, &fail); if (fail) { s_wsle(&io___18); do_lio(&c__9, &c__1, "CPOLY HAS FAILED ON THIS EXAMPLE", 32L); e_wsle(); } prtz_(&c__4, zr, zi); pchk_(p, pi, &c__4, zr, zi); s_wsle(&io___19); do_lio(&c__9, &c__1, "EXAMPLE 1a. POLYNOMIAL WITH ZEROS 1,2,...,10.", 46L); e_wsle(); p[0] = 1.; p[1] = -55.; p[2] = 1320.; p[3] = -18150.; p[4] = 157773.; p[5] = -902055.; p[6] = 3416930.; p[7] = -8409500.; p[8] = 12753576.; p[9] = -10628640.; p[10] = 3628800.; for (i = 1; i <= 11; ++i) { /* L11: */ pi[i - 1] = 0.; } prtc_(&c__11, p, pi); cpoly_(p, pi, &c__10, zr, zi, &fail); if (fail) { s_wsle(&io___20); do_lio(&c__9, &c__1, "CPOLY HAS FAILED ON THIS EXAMPLE", 32L); e_wsle(); } prtz_(&c__10, zr, zi); pchk_(p, pi, &c__10, zr, zi); s_wsle(&io___21); do_lio(&c__9, &c__1, "EXAMPLE 2. ZEROS ON IMAGINARY AXIS DEGREE 3.", 44L); e_wsle(); p[0] = 1.; p[1] = 0.; p[2] = -10001.0001; p[3] = 0.; pi[0] = 0.; pi[1] = -10001.0001; pi[2] = 0.; pi[3] = 1.; prtc_(&c__4, p, pi); cpoly_(p, pi, &c__3, zr, zi, &fail); if (fail) { s_wsle(&io___22); do_lio(&c__9, &c__1, "CPOLY HAS FAILED ON THIS EXAMPLE", 32L); e_wsle(); } prtz_(&c__3, zr, zi); pchk_(p, pi, &c__3, zr, zi); s_wsle(&io___23); do_lio(&c__9, &c__1, "EXAMPLE 3. ZEROS AT 1+I,1/2*(1+I) . 1/(2**-9)*(1+I)" , 51L); e_wsle(); p[0] = 1.f; p[1] = -1.998046875f; p[2] = 0.f; p[3] = .7567065954208374; p[4] = -.2002119533717632; p[5] = .01271507365163416; p[6] = 0.; p[7] = -1.154642632172909e-5; p[8] = 1.584803612786345e-7; p[9] = -4.652065399568528e-10; p[10] = 0.; pi[0] = 0.; pi[1] = p[1]; pi[2] = 2.658859252929688; pi[3] = -.7567065954208374; pi[4] = 0.; pi[5] = p[5]; pi[6] = -7.820779428584501e-4; pi[7] = -p[7]; pi[8] = 0.; pi[9] = p[9]; pi[10] = 9.094947017729282e-13; prtc_(&c__11, p, pi); cpoly_(p, pi, &c__10, zr, zi, &fail); if (fail) { s_wsle(&io___24); do_lio(&c__9, &c__1, "CPOLY HAS FAILED ON THIS EXAMPLE", 32L); e_wsle(); } prtz_(&c__10, zr, zi); pchk_(p, pi, &c__10, zr, zi); s_wsle(&io___25); do_lio(&c__9, &c__1, "EXAMPLE 4. MULTIPLE ZEROS", 25L); e_wsle(); p[0] = 1.; p[1] = -10.; p[2] = 3.; p[3] = 284.; p[4] = -1293.; p[5] = 2374.; p[6] = -1587.; p[7] = -920.; p[8] = 2204.; p[9] = -1344.; p[10] = 288.; pi[0] = 0.; pi[1] = -10.; pi[2] = 100.; pi[3] = -334.; pi[4] = 200.; pi[5] = 1394.; pi[6] = -3836.; pi[7] = 4334.; pi[8] = -2352.; pi[9] = 504.; pi[10] = 0.; prtc_(&c__11, p, pi); cpoly_(p, pi, &c__10, zr, zi, &fail); if (fail) { s_wsle(&io___26); do_lio(&c__9, &c__1, "CPOLY HAS FAILED ON THIS EXAMPLE", 32L); e_wsle(); } prtz_(&c__10, zr, zi); pchk_(p, pi, &c__10, zr, zi); s_wsle(&io___27); do_lio(&c__9, &c__1, "EXAMPLE 5. 12 ZEROS EVENLY DISTRIBUTE ON A CIRCLE OF", 52L); do_lio(&c__9, &c__1, " RADIUS 1. CENTERED AT 0+2I.", 28L); e_wsle(); p[0] = 1.; p[1] = 0.; p[2] = -264.; p[3] = 0.; p[4] = 7920.; p[5] = 0.; p[6] = -59136.; p[7] = 0.; p[8] = 126720.; p[9] = 0.; p[10] = -67584.; p[11] = 0.; p[12] = 4095.; pi[0] = 0.; pi[1] = -24.; pi[2] = 0.; pi[3] = 1760.; pi[4] = 0.; pi[5] = -25344.; pi[6] = 0.; pi[7] = 101376.; pi[8] = 0.; pi[9] = -112640.; pi[10] = 0.; pi[11] = 24576.; pi[12] = 0.; prtc_(&c__13, p, pi); cpoly_(p, pi, &c__12, zr, zi, &fail); if (fail) { s_wsle(&io___28); do_lio(&c__9, &c__1, "CPOLY HAS FAILED ON THIS EXAMPLE", 32L); e_wsle(); } prtz_(&c__12, zr, zi); pchk_(p, pi, &c__12, zr, zi); s_wsle(&io___29); do_lio(&c__9, &c__1, "EXAMPLE 6. 12 ZEROS EVENLY DISTRIBUTE ON A CIRCLE OF", 52L); do_lio(&c__9, &c__1, " RADIUS 5. CENTERED AT 0,0", 26L); e_wsle(); p[0] = 1.; p[1] = 0.; p[2] = 0.; p[3] = 0.; p[4] = 0.; p[5] = 0.; p[6] = 0.; p[7] = 0.; p[8] = 0.; p[9] = 0.; p[10] = 0.; p[11] = 0.; p[12] = 0.; p[13] = 0.; p[14] = 0.; p[15] = 0.; p[16] = -152587890625.f; pi[0] = 0.; pi[1] = 0.; pi[2] = 0.; pi[3] = 0.; pi[4] = 0.; pi[5] = 0.; pi[6] = 0.; pi[7] = 0.; pi[8] = 0.; pi[9] = 0.; pi[10] = 0.; pi[11] = 0.; pi[12] = 0.; pi[13] = 0.; pi[14] = 0.; pi[15] = 0.; pi[16] = 0.; prtc_(&c__17, p, pi); cpoly_(p, pi, &c__16, zr, zi, &fail); if (fail) { s_wsle(&io___30); do_lio(&c__9, &c__1, "CPOLY HAS FAILED ON THIS EXAMPLE", 32L); e_wsle(); } prtz_(&c__16, zr, zi); pchk_(p, pi, &c__16, zr, zi); return 0; } /* MAIN__ */ /* Subroutine */ int calct_(logical *bool) { /* System generated locals */ doublereal d__1, d__2; /* Local variables */ extern doublereal cmod_(doublereal *, doublereal *); static integer n; extern /* Subroutine */ int cdivid_(doublereal *, doublereal *, doublereal *, doublereal *, doublereal *, doublereal *), polyev_( integer *, doublereal *, doublereal *, doublereal *, doublereal *, doublereal *, doublereal *, doublereal *, doublereal *); static doublereal hvi, hvr; /* COMPUTES T = -P(S)/H(S). */ /* BOOL - LOGICAL, SET TRUE IF H(S) IS ESSENTIALLY ZERO. */ /* COMMON AREA */ n = global_1.nn - 1; /* EVALUATE H(S). */ polyev_(&n, &global_1.sr, &global_1.si, global_1.hr, global_1.hi, global_1.qhr, global_1.qhi, &hvr, &hvi); *bool = cmod_(&hvr, &hvi) <= global_1.are * 10. * cmod_(&global_1.hr[n - 1], &global_1.hi[n - 1]); if (*bool) { goto L10; } d__1 = -global_1.pvr; d__2 = -global_1.pvi; cdivid_(&d__1, &d__2, &hvr, &hvi, &global_1.tr, &global_1.ti); return 0; L10: global_1.tr = 0.; global_1.ti = 0.; return 0; } /* calct_ */ doublereal cauchy_(integer *nn, doublereal *pt, doublereal *q) { /* System generated locals */ integer i__1; doublereal ret_val, d__1; /* Builtin functions */ double log(doublereal), exp(doublereal); /* Local variables */ static doublereal f; static integer i, n; static doublereal x, df, dx, xm; /* CAUCHY COMPUTES A LOWER BOUND ON THE MODULI OF THE ZEROS OF A */ /* POLYNOMIAL - PT IS THE MODULUS OF THE COEFFICIENTS. */ /* Parameter adjustments */ --q; --pt; /* Function Body */ pt[*nn] = -pt[*nn]; /* COMPUTE UPPER ESTIMATE OF BOUND. */ n = *nn - 1; x = exp((log(-pt[*nn]) - log(pt[1])) / (real) n); if (pt[n] == 0.) { goto L20; } /* IF NEWTON STEP AT THE ORIGIN IS BETTER, USE IT. */ xm = -pt[*nn] / pt[n]; if (xm < x) { x = xm; } /* CHOP THE INTERVAL (0,X) UNITL F LE 0. */ L20: xm = x * .1; f = pt[1]; i__1 = *nn; for (i = 2; i <= i__1; ++i) { f = f * xm + pt[i]; /* L30: */ } if (f <= 0.) { goto L40; } x = xm; goto L20; L40: dx = x; /* DO NEWTON ITERATION UNTIL X CONVERGES TO TWO DECIMAL PLACES. */ L50: if ((d__1 = dx / x, abs(d__1)) <= .005) { goto L70; } q[1] = pt[1]; i__1 = *nn; for (i = 2; i <= i__1; ++i) { q[i] = q[i - 1] * x + pt[i]; /* L60: */ } f = q[*nn]; df = q[1]; i__1 = n; for (i = 2; i <= i__1; ++i) { df = df * x + q[i]; /* L65: */ } dx = f / df; x -= dx; goto L50; L70: ret_val = x; return ret_val; } /* cauchy_ */ /* Subroutine */ int cdivid_(doublereal *ar, doublereal *ai, doublereal *br, doublereal *bi, doublereal *cr, doublereal *ci) { extern /* Subroutine */ int mcon_(doublereal *, doublereal *, doublereal * , doublereal *); static doublereal d, r, t, infin; /* COMPLEX DIVISION C = A/B, AVOIDING OVERFLOW. */ if (*br != 0. || *bi != 0.) { goto L10; } /* DIVISION BY ZERO, C = INFINITY. */ mcon_(&t, &infin, &t, &t); *cr = infin; *ci = infin; return 0; L10: if (abs(*br) >= abs(*bi)) { goto L20; } r = *br / *bi; d = *bi + r * *br; *cr = (*ar * r + *ai) / d; *ci = (*ai * r - *ar) / d; return 0; L20: r = *bi / *br; d = *br + r * *bi; *cr = (*ar + *ai * r) / d; *ci = (*ai - *ar * r) / d; return 0; } /* cdivid_ */ doublereal cmod_(doublereal *r, doublereal *i) { /* System generated locals */ doublereal ret_val, d__1; /* Builtin functions */ double sqrt(doublereal); /* Local variables */ static doublereal ai, ar; /* MODULUS OF A COMPLEX NUMBER AVOIDING OVERFLOW. */ ar = abs(*r); ai = abs(*i); if (ar >= ai) { goto L10; } /* Computing 2nd power */ d__1 = ar / ai; ret_val = ai * sqrt(d__1 * d__1 + 1.); return ret_val; L10: if (ar <= ai) { goto L20; } /* Computing 2nd power */ d__1 = ai / ar; ret_val = ar * sqrt(d__1 * d__1 + 1.); return ret_val; L20: ret_val = ar * sqrt(2.); return ret_val; } /* cmod_ */ /* Subroutine */ int cpoly_(doublereal *opr, doublereal *opi, integer *degree, doublereal *zeror, doublereal *zeroi, logical *fail) { /* System generated locals */ integer i__1; doublereal d__1, d__2; /* Builtin functions */ double sqrt(doublereal); /* Local variables */ static doublereal base; extern doublereal cmod_(doublereal *, doublereal *); extern /* Subroutine */ int mcon_(doublereal *, doublereal *, doublereal * , doublereal *); static logical conv; static doublereal cosr, sinr; static integer idnn2, i; extern doublereal scale_(integer *, doublereal *, doublereal *, doublereal *, doublereal *, doublereal *); static doublereal zi; extern /* Subroutine */ int cdivid_(doublereal *, doublereal *, doublereal *, doublereal *, doublereal *, doublereal *); static doublereal zr; extern doublereal cauchy_(integer *, doublereal *, doublereal *); static doublereal xx, yy, smalno; extern /* Subroutine */ int noshft_(integer *), fxshft_(integer *, doublereal *, doublereal *, logical *); static doublereal bnd, xxx; static integer cnt1, cnt2; /* FINDS THE ZEROS OF A COMPLEX POLYNOMIAL. */ /* OPR, OPI - DOUBLE PRECISION VECTORS OF REAL AND */ /* IMAGINARY PARTS OF THE COEFFICIENTS IN */ /* ORDER OF DECREASING POWERS. */ /* DEGREE - INTEGER DEGREE OF POLYNOMIAL. */ /* ZEROR, ZEROI - OUTPUT DOUBLE PRECISION 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. */ /* COMMON AREA */ /* TO CHANGE THE SIZE OF POLYNOMIALS WHICH CAN BE SOLVED, REPLACE */ /* THE DIMENSION OF THE ARRAYS IN THE COMMON AREA. */ /* INITIALIZATION OF CONSTANTS */ /* Parameter adjustments */ --zeroi; --zeror; --opi; --opr; /* Function Body */ mcon_(&global_1.eta, &global_1.infin, &smalno, &base); global_1.are = global_1.eta; global_1.mre = sqrt(2.) * 2. * global_1.eta; xx = .70710678f; yy = -xx; cosr = -.060756474f; sinr = .99756405f; *fail = FALSE_; global_1.tr = 0.; global_1.ti = 0.; global_1.sr = 0.; global_1.si = 0.; global_1.nn = *degree + 1; /* ALGORITHM FAILS IF THE LEADING COEFFICIENT IS ZERO. */ if (opr[1] != 0. || opi[1] != 0.) { goto L10; } *fail = TRUE_; return 0; /* REMOVE THE ZEROS AT THE ORIGIN IF ANY. */ L10: if (opr[global_1.nn] != 0. || opi[global_1.nn] != 0.) { goto L20; } idnn2 = *degree - global_1.nn + 2; zeror[idnn2] = 0.; zeroi[idnn2] = 0.; --global_1.nn; goto L10; /* MAKE A COPY OF THE COEFFICIENTS. */ L20: i__1 = global_1.nn; for (i = 1; i <= i__1; ++i) { global_1.pr[i - 1] = opr[i]; global_1.pi[i - 1] = opi[i]; global_1.shr[i - 1] = cmod_(&global_1.pr[i - 1], &global_1.pi[i - 1]); /* L30: */ } /* SCALE THE POLYNOMIAL. */ bnd = scale_(&global_1.nn, global_1.shr, &global_1.eta, &global_1.infin, & smalno, &base); if (bnd == 1.) { goto L40; } i__1 = global_1.nn; for (i = 1; i <= i__1; ++i) { global_1.pr[i - 1] = bnd * global_1.pr[i - 1]; global_1.pi[i - 1] = bnd * global_1.pi[i - 1]; /* L35: */ } /* START THE ALGORITHM FOR ONE ZERO . */ L40: if (global_1.nn > 2) { goto L50; } /* CALCULATE THE FINAL ZERO AND RETURN. */ d__1 = -global_1.pr[1]; d__2 = -global_1.pi[1]; cdivid_(&d__1, &d__2, global_1.pr, global_1.pi, &zeror[*degree], &zeroi[* degree]); return 0; /* CALCULATE BND, A LOWER BOUND ON THE MODULUS OF THE ZEROS. */ L50: i__1 = global_1.nn; for (i = 1; i <= i__1; ++i) { global_1.shr[i - 1] = cmod_(&global_1.pr[i - 1], &global_1.pi[i - 1]); /* L60: */ } bnd = cauchy_(&global_1.nn, global_1.shr, global_1.shi); /* OUTER LOOP TO CONTROL 2 MAJOR PASSES WITH DIFFERENT SEQUENCES */ /* OF SHIFTS. */ for (cnt1 = 1; cnt1 <= 2; ++cnt1) { /* FIRST STAGE CALCULATION, NO SHIFT. */ noshft_(&c__5); /* INNER LOOP TO SELECT A SHIFT. */ for (cnt2 = 1; cnt2 <= 9; ++cnt2) { /* 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; global_1.sr = bnd * xx; global_1.si = bnd * yy; /* SECOND STAGE CALCULATION, FIXED SHIFT. */ i__1 = cnt2 * 10; fxshft_(&i__1, &zr, &zi, &conv); if (! conv) { goto L80; } /* THE SECOND STAGE JUMPS DIRECTLY TO THE THIRD STAGE ITERATION. */ /* IF SUCCESSFUL THE ZERO IS STORED AND THE POLYNOMIAL DEFLATED. */ idnn2 = *degree - global_1.nn + 2; zeror[idnn2] = zr; zeroi[idnn2] = zi; --global_1.nn; i__1 = global_1.nn; for (i = 1; i <= i__1; ++i) { global_1.pr[i - 1] = global_1.qpr[i - 1]; global_1.pi[i - 1] = global_1.qpi[i - 1]; /* L70: */ } goto L40; L80: /* IF THE ITERATION IS UNSUCCESSFUL ANOTHER SHIFT IS CHOSEN. */ /* L90: */ ; } /* IF 9 SHIFTS FAIL, THE OUTER LOOP IS REPEATED WITH ANOTHER */ /* SEQUENCE OF SHIFTS. */ /* L100: */ } /* THE ZEROFINDER HAS FAILED ON TWO MAJOR PASSES. */ /* RETURN EMPTY HANDED. */ *fail = TRUE_; return 0; } /* cpoly_ */ doublereal errev_(integer *nn, doublereal *qr, doublereal *qi, doublereal *ms, doublereal *mp, doublereal *are, doublereal *mre) { /* System generated locals */ integer i__1; doublereal ret_val; /* Local variables */ extern doublereal cmod_(doublereal *, doublereal *); static doublereal e; static integer i; /* 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 */ /* Parameter adjustments */ --qi; --qr; /* Function Body */ e = cmod_(&qr[1], &qi[1]) * *mre / (*are + *mre); /* changed 1,NN to 2,NN */ i__1 = *nn; for (i = 2; i <= i__1; ++i) { e = e * *ms + cmod_(&qr[i], &qi[i]); /* L10: */ } ret_val = e * (*are + *mre) - *mp * *mre; return ret_val; } /* errev_ */ /* Subroutine */ int fxshft_(integer *l2, doublereal *zr, doublereal *zi, logical *conv) { /* System generated locals */ integer i__1, i__2; doublereal d__1, d__2; /* Local variables */ extern doublereal cmod_(doublereal *, doublereal *); static logical pasd, bool, test; static doublereal svsi, svsr; static integer i, j, n; extern /* Subroutine */ int calct_(logical *), nexth_(logical *), vrshft_( integer *, doublereal *, doublereal *, logical *), polyev_( integer *, doublereal *, doublereal *, doublereal *, doublereal *, doublereal *, doublereal *, doublereal *, doublereal *); static doublereal oti, otr; /* 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 */ /* COMMON AREA */ n = global_1.nn - 1; /* EVALUATE P AT S. */ polyev_(&global_1.nn, &global_1.sr, &global_1.si, global_1.pr, global_1.pi, global_1.qpr, global_1.qpi, &global_1.pvr, & global_1.pvi); test = TRUE_; pasd = FALSE_; /* CALCULATE FIRST T = -P(S)/H(S). */ calct_(&bool); /* MAIN LOOP FOR ONE SECOND STAGE STEP. */ i__1 = *l2; for (j = 1; j <= i__1; ++j) { otr = global_1.tr; oti = global_1.ti; /* COMPUTE NEXT H POLYNOMIAL AND NEW T. */ nexth_(&bool); calct_(&bool); *zr = global_1.sr + global_1.tr; *zi = global_1.si + global_1.ti; /* TEST FOR CONVERGENCE UNLESS STAGE 3 HAS FAILED ONCE OR THIS */ /* IS THE LAST H POLYNOMIAL . */ if (bool || ! test || j == *l2) { goto L50; } d__1 = global_1.tr - otr; d__2 = global_1.ti - oti; if (cmod_(&d__1, &d__2) >= cmod_(zr, zi) * .5) { goto L40; } if (! pasd) { goto L30; } /* THE WEAK CONVERGENCE TEST HAS BEEN PASSED TWICE, START THE */ /* THIRD STAGE ITERATION, AFTER SAVING THE CURRENT H POLYNOMIAL */ /* AND SHIFT. */ i__2 = n; for (i = 1; i <= i__2; ++i) { global_1.shr[i - 1] = global_1.hr[i - 1]; global_1.shi[i - 1] = global_1.hi[i - 1]; /* L10: */ } svsr = global_1.sr; svsi = global_1.si; vrshft_(&c__10, zr, zi, conv); if (*conv) { return 0; } /* THE ITERATION FAILED TO CONVERGE. TURN OFF TESTING AND RESTORE */ /* H,S,PV AND T. */ test = FALSE_; i__2 = n; for (i = 1; i <= i__2; ++i) { global_1.hr[i - 1] = global_1.shr[i - 1]; global_1.hi[i - 1] = global_1.shi[i - 1]; /* L20: */ } global_1.sr = svsr; global_1.si = svsi; polyev_(&global_1.nn, &global_1.sr, &global_1.si, global_1.pr, global_1.pi, global_1.qpr, global_1.qpi, &global_1.pvr, & global_1.pvi); calct_(&bool); goto L50; L30: pasd = TRUE_; goto L50; L40: pasd = FALSE_; L50: ; } /* ATTEMPT AN ITERATION WITH FINAL H POLYNOMIAL FROM SECOND STAGE. */ vrshft_(&c__10, zr, zi, conv); return 0; } /* fxshft_ */ /* Subroutine */ int mcon_(doublereal *eta, doublereal *infiny, doublereal * smalno, doublereal *base) { /* System generated locals */ integer i__1, i__2; doublereal d__1, d__2; /* Builtin functions */ double pow_di(doublereal *, integer *); /* Local variables */ static integer m, n, t; /* 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.0D0 + ETA IS */ /* GREATER THAN 1.0D0. */ /* INFINY THE LARGEST FLOATING-POINT NUMBER */ /* SMALNO THE SMALLEST POSITIVE FLOATING-POINT NUMBER */ /* BASE THE BASE OF THE FLOATING-POINT NUMBER SYSTEM USED */ /* LET T BE THE NUMBER OF BASE-DIGITS IN EACH FLOATING-POINT */ /* NUMBER(DOUBLE PRECISION). THEN ETA IS EITHER .5*B**(1-T) */ /* OR B**(1-T) DEPENDING ON WHETHER ROUNDING OR TRUNCATION */ /* IS USED. */ /* LET M BE THE LARGEST EXPONENT AND N THE SMALLEST EXPONENT */ /* IN THE NUMBER SYSTEM. THEN INFINY IS (1-BASE**(-T))*BASE**M */ /* AND SMALNO IS BASE**N. */ /* THE VALUES FOR BASE,T,M,N BELOW CORRESPOND TO THE VAX */ *base = 2.; t = 52; m = 127; n = -127; i__1 = 1 - t; *eta = pow_di(base, &i__1); i__1 = -t; i__2 = m - 1; *infiny = *base * (1. - pow_di(base, &i__1)) * pow_di(base, &i__2); i__1 = n + 3; /* Computing 3rd power */ d__1 = *base, d__2 = d__1; *smalno = pow_di(base, &i__1) / (d__2 * (d__1 * d__1)); return 0; } /* mcon_ */ /* Subroutine */ int nexth_(logical *bool) { /* System generated locals */ integer i__1; /* Local variables */ static integer j, n; static doublereal t1, t2; /* CALCULATES THE NEXT SHIFTED H POLYNOMIAL. */ /* BOOL - LOGICAL, IF .TRUE. H(S) IS ESSENTIALLY ZERO */ /* COMMON AREA */ n = global_1.nn - 1; if (*bool) { goto L20; } i__1 = n; for (j = 2; j <= i__1; ++j) { t1 = global_1.qhr[j - 2]; t2 = global_1.qhi[j - 2]; global_1.hr[j - 1] = global_1.tr * t1 - global_1.ti * t2 + global_1.qpr[j - 1]; global_1.hi[j - 1] = global_1.tr * t2 + global_1.ti * t1 + global_1.qpi[j - 1]; /* L10: */ } global_1.hr[0] = global_1.qpr[0]; global_1.hi[0] = global_1.qpi[0]; return 0; /* IF H(S) IS ZERO REPLACE H WITH QH. */ L20: i__1 = n; for (j = 2; j <= i__1; ++j) { global_1.hr[j - 1] = global_1.qhr[j - 2]; global_1.hi[j - 1] = global_1.qhi[j - 2]; /* L30: */ } global_1.hr[0] = 0.; global_1.hi[0] = 0.; return 0; } /* nexth_ */ /* Subroutine */ int noshft_(integer *l1) { /* System generated locals */ integer i__1, i__2; doublereal d__1, d__2; /* Local variables */ extern doublereal cmod_(doublereal *, doublereal *); static integer i, j, n; static doublereal t1, t2; static integer jj; extern /* Subroutine */ int cdivid_(doublereal *, doublereal *, doublereal *, doublereal *, doublereal *, doublereal *); static integer nm1; static doublereal xni; /* COMPUTES THE DERIVATIVE POLYNOMIAL AS THE INITIAL H */ /* POLYNOMIAL AND COMPUTES L1 NO-SHIFT H POLYNOMIALS. */ /* COMMON AREA */ n = global_1.nn - 1; nm1 = n - 1; i__1 = n; for (i = 1; i <= i__1; ++i) { xni = (doublereal) (global_1.nn - i); global_1.hr[i - 1] = xni * global_1.pr[i - 1] / (real) n; global_1.hi[i - 1] = xni * global_1.pi[i - 1] / (real) n; /* L10: */ } i__1 = *l1; for (jj = 1; jj <= i__1; ++jj) { if (cmod_(&global_1.hr[n - 1], &global_1.hi[n - 1]) <= global_1.eta * 10. * cmod_(&global_1.pr[n - 1], &global_1.pi[n - 1])) { goto L30; } d__1 = -global_1.pr[global_1.nn - 1]; d__2 = -global_1.pi[global_1.nn - 1]; cdivid_(&d__1, &d__2, &global_1.hr[n - 1], &global_1.hi[n - 1], & global_1.tr, &global_1.ti); i__2 = nm1; for (i = 1; i <= i__2; ++i) { j = global_1.nn - i; t1 = global_1.hr[j - 2]; t2 = global_1.hi[j - 2]; global_1.hr[j - 1] = global_1.tr * t1 - global_1.ti * t2 + global_1.pr[j - 1]; global_1.hi[j - 1] = global_1.tr * t2 + global_1.ti * t1 + global_1.pi[j - 1]; /* L20: */ } global_1.hr[0] = global_1.pr[0]; global_1.hi[0] = global_1.pi[0]; goto L50; /* IF THE CONSTANT TERM IS ESSENTIALLY ZERO, SHIFT H COEFFICIENTS. */ L30: i__2 = nm1; for (i = 1; i <= i__2; ++i) { j = global_1.nn - i; global_1.hr[j - 1] = global_1.hr[j - 2]; global_1.hi[j - 1] = global_1.hi[j - 2]; /* L40: */ } global_1.hr[0] = 0.; global_1.hi[0] = 0.; L50: ; } return 0; } /* noshft_ */ /* Subroutine */ int polyev_(integer *nn, doublereal *sr, doublereal *si, doublereal *pr, doublereal *pi, doublereal *qr, doublereal *qi, doublereal *pvr, doublereal *pvi) { /* System generated locals */ integer i__1; /* Local variables */ static integer i; static doublereal t; /* EVALUATES A POLYNOMIAL P AT S BY THE HORNER RECURRENCE */ /* PLACING THE PARTIAL SUMS IN Q AND THE COMPUTED VALUE IN PV. */ /* Parameter adjustments */ --qi; --qr; --pi; --pr; /* Function Body */ qr[1] = pr[1]; qi[1] = pi[1]; *pvr = qr[1]; *pvi = qi[1]; i__1 = *nn; for (i = 2; i <= i__1; ++i) { t = *pvr * *sr - *pvi * *si + pr[i]; *pvi = *pvr * *si + *pvi * *sr + pi[i]; *pvr = t; qr[i] = *pvr; qi[i] = *pvi; /* L10: */ } return 0; } /* polyev_ */ /* Subroutine */ int pchk_(doublereal *opr, doublereal *opi, integer *degree, doublereal *zeror, doublereal *zeroi) { /* Format strings */ static char fmt_100[] = "(\002CHECK\002,/50(2d26.16/))"; /* System generated locals */ integer i__1, i__2; /* Builtin functions */ integer s_wsfe(cilist *), do_fio(integer *, char *, ftnlen), e_wsfe(void); /* Local variables */ static doublereal temp, sumi, sumr; static integer i, j; /* Fortran I/O blocks */ static cilist io___96 = { 0, 6, 0, fmt_100, 0 }; /* CHECKS THE ZEROS OF A COMPLEX POLYNOMIAL. */ /* OPR, OPI - DOUBLE PRECISION VECTORS OF REAL AND */ /* IMAGINARY PARTS OF THE COEFFICIENTS IN */ /* ORDER OF DECREASING POWERS. */ /* DEGREE - INTEGER DEGREE OF POLYNOMIAL. */ /* ZEROR, ZEROI - DOUBLE PRECISION VECTORS OF */ /* REAL AND IMAGINARY PARTS OF THE ZEROS. */ /* REAL*16 SUMR,SUMI,TEMP */ /* Parameter adjustments */ --zeroi; --zeror; --opi; --opr; /* Function Body */ i__1 = *degree; for (i = 1; i <= i__1; ++i) { sumr = 0.; sumi = 0.; i__2 = *degree + 1; for (j = 1; j <= i__2; ++j) { temp = opr[j] + sumr * zeror[i] - sumi * zeroi[i]; sumi = opi[j] + sumr * zeroi[i] + sumi * zeror[i]; sumr = temp; /* L20: */ } zeror[i] = sumr; zeroi[i] = sumi; /* L10: */ } s_wsfe(&io___96); i__1 = *degree; for (i = 1; i <= i__1; ++i) { do_fio(&c__1, (char *)&zeror[i], (ftnlen)sizeof(doublereal)); do_fio(&c__1, (char *)&zeroi[i], (ftnlen)sizeof(doublereal)); } e_wsfe(); return 0; } /* pchk_ */ /* Subroutine */ int prtc_(integer *n, doublereal *p, doublereal *q) { /* Format strings */ static char fmt_10[] = "(\002COEFFICIENTS\002,/50(2d26.16/))"; /* System generated locals */ integer i__1; /* Builtin functions */ integer s_wsfe(cilist *), do_fio(integer *, char *, ftnlen), e_wsfe(void); /* Local variables */ static integer i; /* Fortran I/O blocks */ static cilist io___97 = { 0, 6, 0, fmt_10, 0 }; /* Parameter adjustments */ --q; --p; /* Function Body */ s_wsfe(&io___97); i__1 = *n; for (i = 1; i <= i__1; ++i) { do_fio(&c__1, (char *)&p[i], (ftnlen)sizeof(doublereal)); do_fio(&c__1, (char *)&q[i], (ftnlen)sizeof(doublereal)); } e_wsfe(); return 0; } /* prtc_ */ /* Subroutine */ int prtz_(integer *n, doublereal *zr, doublereal *zi) { /* Format strings */ static char fmt_10[] = "(\002ZEROS\002,/50(2d26.16/))"; /* System generated locals */ integer i__1; /* Builtin functions */ integer s_wsfe(cilist *), do_fio(integer *, char *, ftnlen), e_wsfe(void); /* Local variables */ static integer i; /* Fortran I/O blocks */ static cilist io___99 = { 0, 6, 0, fmt_10, 0 }; /* Parameter adjustments */ --zi; --zr; /* Function Body */ s_wsfe(&io___99); i__1 = *n; for (i = 1; i <= i__1; ++i) { do_fio(&c__1, (char *)&zr[i], (ftnlen)sizeof(doublereal)); do_fio(&c__1, (char *)&zi[i], (ftnlen)sizeof(doublereal)); } e_wsfe(); return 0; } /* prtz_ */ doublereal scale_(integer *nn, doublereal *pt, doublereal *eta, doublereal * infin, doublereal *smalno, doublereal *base) { /* System generated locals */ integer i__1; doublereal ret_val; /* Builtin functions */ double sqrt(doublereal), log(doublereal), pow_di(doublereal *, integer *); /* Local variables */ static integer i, l; static doublereal x, hi, sc, lo, min_, max_; /* 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. */ /* FIND LARGEST AND SMALLEST MODULI OF COEFFICIENTS. */ /* Parameter adjustments */ --pt; /* Function Body */ hi = sqrt(*infin); lo = *smalno / *eta; max_ = 0.; min_ = *infin; i__1 = *nn; for (i = 1; i <= i__1; ++i) { x = pt[i]; if (x > max_) { max_ = x; } if (x != 0. && x < min_) { min_ = x; } /* L10: */ } /* SCALE ONLY IF THERE ARE VERY LARGE OR VERY SMALL COMPONENTS. */ ret_val = 1.; if (min_ >= lo && max_ <= hi) { return ret_val; } x = lo / min_; if (x > 1.) { goto L20; } sc = 1. / (sqrt(max_) * sqrt(min_)); goto L30; L20: sc = x; if (*infin / sc > max_) { sc = 1.; } L30: l = (integer) (log(sc) / log(*base) + .5f); ret_val = pow_di(base, &l); return ret_val; } /* scale_ */ /* Subroutine */ int vrshft_(integer *l3, doublereal *zr, doublereal *zi, logical *conv) { /* System generated locals */ integer i__1; /* Builtin functions */ double sqrt(doublereal); /* Local variables */ extern doublereal cmod_(doublereal *, doublereal *); static logical bool, b; static integer i, j; extern /* Subroutine */ int calct_(logical *); extern doublereal errev_(integer *, doublereal *, doublereal *, doublereal *, doublereal *, doublereal *, doublereal *); extern /* Subroutine */ int nexth_(logical *); static doublereal r1, r2, mp, ms, tp, relstp; extern /* Subroutine */ int polyev_(integer *, doublereal *, doublereal *, doublereal *, doublereal *, doublereal *, doublereal *, doublereal *, doublereal *); static doublereal omp; /* 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 */ /* COMMON AREA */ *conv = FALSE_; b = FALSE_; global_1.sr = *zr; global_1.si = *zi; /* MAIN LOOP FOR STAGE THREE */ i__1 = *l3; for (i = 1; i <= i__1; ++i) { /* EVALUATE P AT S AND TEST FOR CONVERGENCE. */ polyev_(&global_1.nn, &global_1.sr, &global_1.si, global_1.pr, global_1.pi, global_1.qpr, global_1.qpi, &global_1.pvr, & global_1.pvi); mp = cmod_(&global_1.pvr, &global_1.pvi); ms = cmod_(&global_1.sr, &global_1.si); if (mp > errev_(&global_1.nn, global_1.qpr, global_1.qpi, &ms, &mp, & global_1.are, &global_1.mre) * 20.) { goto L10; } /* POLYNOMIAL VALUE IS SMALLER IN VALUE THAN A BOUND ON THE ERROR */ /* IN EVALUATING P, TERMINATE THE ITERATION. */ *conv = TRUE_; *zr = global_1.sr; *zi = global_1.si; return 0; L10: if (i == 1) { goto L40; } if (b || mp < omp || relstp >= .05) { goto L30; } /* 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 < global_1.eta) { tp = global_1.eta; } r1 = sqrt(tp); r2 = global_1.sr * (r1 + 1.) - global_1.si * r1; global_1.si = global_1.sr * r1 + global_1.si * (r1 + 1.); global_1.sr = r2; polyev_(&global_1.nn, &global_1.sr, &global_1.si, global_1.pr, global_1.pi, global_1.qpr, global_1.qpi, &global_1.pvr, & global_1.pvi); for (j = 1; j <= 5; ++j) { calct_(&bool); nexth_(&bool); /* L20: */ } omp = global_1.infin; goto L50; /* EXIT IF POLYNOMIAL VALUE INCREASES SIGNIFICANTLY. */ L30: if (mp * .1 > omp) { return 0; } L40: omp = mp; /* CALCULATE NEXT ITERATE. */ L50: calct_(&bool); nexth_(&bool); calct_(&bool); if (bool) { goto L60; } relstp = cmod_(&global_1.tr, &global_1.ti) / cmod_(&global_1.sr, & global_1.si); global_1.sr += global_1.tr; global_1.si += global_1.ti; L60: ; } return 0; } /* vrshft_ */ /* Main program alias */ int croots_ () { MAIN__ (); return 0; }