// c419.java basic conversion from Fortran IV, cpoly // test cases and error check appear first // many Fortran comments included for traceability package myjava; import java.text.*; // to get formatting public class c419 { // ALGORITHM 419 COLLECTED ALGORITHMS FROM ACM. // ALGORITHM APPEARED IN COMM. ACM, VOL. 15, NO. 02, // P. 097. // PROGRAM CROOTS // // DRIVER TO TEST CPOLY // // COMMON AREA used by calct // COMMON/GLOBAL/PR,PI,HR,HI,QPR,QPI,QHR,QHI,SHR,SHI, // * SR,SI,TR,TI,PVR,PVI,ARE,MRE,ETA,INFIN,NN // DOUBLE PRECISION 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) static double pr[] = new double[50]; static double pi[] = new double[50]; static double hr[] = new double[50]; static double hi[] = new double[50]; static double qpr[] = new double[50]; static double qpi[] = new double[50]; static double qhr[] = new double[50]; static double qhi[] = new double[50]; static double shr[] = new double[50]; static double shi[] = new double[50]; static double sr; static double si; static double tr; static double ti; static double pvr[]={0.0}; static double pvi[]={0.0}; // typical use as out parameter static double eta=0.0; static double infin=0.0; static double smalno=0.0; static double base=0.0; static double are=0.0; static double mre=0.0; static int nn; static double maxc; public c419() // test program { // LOGICAL FAIL boolean fail[]={false}; // out parameter // DOUBLE PRECISION P(50),PI(50),ZR(50),ZI(50) double p[] = new double[50]; double pi[] = new double[50]; double zr[] = new double[50]; double zi[] = new double[50]; int n; // degree of polynomial, one less than highest subscript System.out.println("c419.java test cpoly to find zeros of polynomials"); mcon(); // initialize COMMON, machine constants System.out.println("machine constants"); System.out.println("eta="+eta); System.out.println("infin="+infin); System.out.println("smalno="+smalno); System.out.println("base="+base); System.out.println("are="+are); System.out.println("mre="+mre); System.out.println(); System.out.println(" real part imaginary part"+ " modulus"); System.out.println(); p[0] = 0.0; // unused, Fortran indexing from 1..n+1 for degree n pi[0] = 0.0; // WRITE(6,*) 'EXAMPLE 0. POLYNOMIAL WITH ZERO 0.5' System.out.println(" EXAMPLE 0. POLYNOMIAL WITH ZERO 0.5"); n = 1; p[1]=4.0; p[2]=-2.0; pi[1]=4.0; pi[2]=-2.0; prtc(n,p,pi); cpoly(p,pi,n,zr,zi,fail); if(fail[0]) System.out.println("CPOLY HAS FAILED ON THIS EXAMPLE"); prtz(n,zr,zi); pchk(p,pi,n,zr,zi); System.out.println(" EXAMPLE 0a. POLYNOMIAL WITH ZEROS 0.0 -2.0"); n = 2; p[1]=1.0; p[2]=2.0; p[3]=0.0; pi[1]=1.0; pi[2]=2.0; pi[3]=0.0; prtc(n,p,pi); cpoly(p,pi,n,zr,zi,fail); if(fail[0]) System.out.println("CPOLY HAS FAILED ON THIS EXAMPLE"); prtz(n,zr,zi); pchk(p,pi,n,zr,zi); System.out.println(" EXAMPLE 0b. POLYNOMIAL WITH ZEROS (1+2i) (3-4i)"); n = 2; p[1]=1.0; p[2]=-4.0; p[3]=11.0; pi[1]=0.0; pi[2]=2.0; pi[3]=2.0; prtc(n,p,pi); cpoly(p,pi,n,zr,zi,fail); if(fail[0]) System.out.println("CPOLY HAS FAILED ON THIS EXAMPLE"); prtz(n,zr,zi); pchk(p,pi,n,zr,zi); System.out.println(" EXAMPLE 0c. POLYNOMIAL WITH ZEROS 10**15 and"+ " 10**(-15)"); n = 2; p[1]=1.0; p[2]=-1.0E15; p[3]=1.0; pi[1]=0.0; pi[2]=0.0; pi[3]=0.0; prtc(n,p,pi); cpoly(p,pi,n,zr,zi,fail); if(fail[0]) System.out.println("CPOLY HAS FAILED ON THIS EXAMPLE"); prtz(n,zr,zi); pchk(p,pi,n,zr,zi); System.out.println(" EXAMPLE 0d. POLYNOMIAL WITH ZEROS"+ " (1+2i) (3-4i) (-6+5i)"); n = 3; p[1]=1.0; p[2]=2.0; p[3]=-3.0; p[4]=76.0; pi[1]=0.0; pi[2]=-3.0; pi[3]=34.0; pi[4]=-43.0; prtc(n,p,pi); cpoly(p,pi,n,zr,zi,fail); if(fail[0]) System.out.println("CPOLY HAS FAILED ON THIS EXAMPLE"); prtz(n,zr,zi); pchk(p,pi,n,zr,zi); System.out.println(" EXAMPLE 1. POLYNOMIAL WITH ZEROS 1,2,3,4"); n = 4; p[1]=1.0; p[2]=-10.0; p[3]=35.0; p[4]=-50.0; p[5]=24.0; for(int i=1; i<=n+1; i++) pi[i]=0.0; prtc(n,p,pi); cpoly(p,pi,n,zr,zi,fail); if(fail[0]) System.out.println("CPOLY HAS FAILED ON THIS EXAMPLE"); prtz(n,zr,zi); pchk(p,pi,n,zr,zi); System.out.println(" EXAMPLE 1a. POLYNOMIAL WITH ZEROS 1,2,...,10."); n = 10; p[1]=1.0; p[2]=-55.0; p[3]=1320.0; p[4]=-18150.0; p[5]=157773.0; p[6]=-902055.0; p[7] = 3416930.0; p[8]=-8409500.0; p[9]=12753576.0; p[10]=-10628640.0; p[11]=3628800.0; for(int i=1; i<=n+1; i++) pi[i]=0.0; prtc(n,p,pi); cpoly(p,pi,n,zr,zi,fail); if(fail[0]) System.out.println("CPOLY HAS FAILED ON THIS EXAMPLE"); prtz(n,zr,zi); pchk(p,pi,n,zr,zi); System.out.println(" EXAMPLE 2. ZEROS ON IMAGINARY AXIS DEGREE 3."); n = 3; p[1]=1.0; p[2]=0.0; p[3]=-10001.0001; p[4]=0.0; pi[1]=0.0; pi[2]=-10001.0001; pi[3]=0.0; pi[4]=1.0; prtc(n,p,pi); cpoly(p,pi,n,zr,zi,fail); if(fail[0]) System.out.println("CPOLY HAS FAILED ON THIS EXAMPLE"); prtz(n,zr,zi); pchk(p,pi,n,zr,zi); System.out.println(" EXAMPLE 3. ZEROS AT 1+I,1/2*(1+I) . 1/(2**-9)*(1+I)"); n = 10; p[1]=1.0; p[2]=-1.998046875; p[3]=0.0; p[4]=0.7567065954208374; p[5]=-0.2002119533717632; p[6]=1.271507365163416E-2; p[7]=0.0; p[8]=-1.154642632172909E-5; p[9]=1.584803612786345E-7; p[10]=-4.652065399568528E-10; p[11]=0.0; pi[1]=0.0; pi[2]=p[2]; pi[3]=2.658859252929688; pi[4]=-7.567065954208374E-1; pi[5]=0.0; pi[6]=p[6]; pi[7]=-7.820779428584501E-4; pi[8]=-p[8]; pi[9]=0.0; pi[10]=p[10]; pi[11]=9.094947017729282E-13; prtc(n,p,pi); cpoly(p,pi,n,zr,zi,fail); if(fail[0]) System.out.println("CPOLY HAS FAILED ON THIS EXAMPLE"); prtz(n,zr,zi); pchk(p,pi,n,zr,zi); System.out.println(" EXAMPLE 4. MULTIPLE ZEROS"); n = 10; p[1]=1.0; p[2]=-10.0; p[3]=3.0; p[4]=284.0; p[5]=-1293.0; p[6]=2374.0; p[7]=-1587.0; p[8]=-920.0; p[9]=2204.0; p[10]=-1344.0; p[11]=288.0; pi[1]=0.0; pi[2]=-10.0; pi[3]=100.0; pi[4]=-334.0; pi[5]=200.0; pi[6]=1394.0; pi[7] =-3836.0; pi[8]=4334.0; pi[9]=-2352.0; pi[10]=504.0; pi[11]=0.0; prtc(n,p,pi); cpoly(p,pi,n,zr,zi,fail); if(fail[0]) System.out.println("CPOLY HAS FAILED ON THIS EXAMPLE"); prtz(n,zr,zi); pchk(p,pi,n,zr,zi); System.out.println(" EXAMPLE 5. 12 ZEROS EVENLY DISTRIBUTE ON A CIRCLE OF"+ " RADIUS 1. CENTERED AT 0+2I."); n = 12; p[1]=1.0; p[2]=0.0; p[3]=-264.0; p[4]=0.0; p[5]=7920.0; p[6]=0.0; p[7]=-59136.0; p[8]=0.0; p[9]=126720.0; p[10]=0.0; p[11]=-67584.0; p[12]=0.0; p[13]=4095.0; pi[1]=0.0; pi[2]=-24.0; pi[3]=0.0; pi[4]=1760.0; pi[5]=0.0; pi[6]=-25344.0; pi[7]=0.0; pi[8]=101376.0; pi[9]=0.0; pi[10]=-112640.0; pi[11]=0.0; pi[12]=24576.0; pi[13]=0.0; prtc(n,p,pi); cpoly(p,pi,n,zr,zi,fail); if(fail[0]) System.out.println("CPOLY HAS FAILED ON THIS EXAMPLE"); prtz(n,zr,zi); pchk(p,pi,n,zr,zi); System.out.println(" EXAMPLE 6. 16 ZEROS EVENLY DISTRIBUTE ON A CIRCLE"+ " OF RADIUS 5. CENTERED AT 0,0"); n = 16; for(int i=1; i<=n+1; i++) {p[i]=0.0; pi[i]=0.0;} p[1]=1.0; p[17]=-Math.pow(5.0,16); prtc(n,p,pi); cpoly(p,pi,n,zr,zi,fail); if(fail[0]) System.out.println("CPOLY HAS FAILED ON THIS EXAMPLE"); prtz(n,zr,zi); pchk(p,pi,n,zr,zi); System.out.println(" EXAMPLE 7a. 12 ZEROS ON A CIRCLE"+ " OF RADIUS 1. CENTERED AT 0,0"); n = 12; for(int i=1; i<=n+1; i++) {p[i]=0.0; pi[i]=0.0;} p[1]=1.0; p[n+1]=1.0; prtc(n,p,pi); cpoly(p,pi,n,zr,zi,fail); if(fail[0]) System.out.println("CPOLY HAS FAILED ON THIS EXAMPLE"); prtz(n,zr,zi); pchk(p,pi,n,zr,zi); System.out.println(" EXAMPLE 7b. 24 ZEROS ON A CIRCLE"+ " OF RADIUS 1. CENTERED AT 0,0"); n = 24; for(int i=1; i<=n+1; i++) {p[i]=0.0; pi[i]=0.0;} p[1]=1.0; p[n+1]=1.0; prtc(n,p,pi); cpoly(p,pi,n,zr,zi,fail); if(fail[0]) System.out.println("CPOLY HAS FAILED ON THIS EXAMPLE"); prtz(n,zr,zi); pchk(p,pi,n,zr,zi); System.out.println(" EXAMPLE 7c. 36 ZEROS ON A CIRCLE"+ " OF RADIUS 1. CENTERED AT 0,0"); n = 36; for(int i=1; i<=n+1; i++) {p[i]=0.0; pi[i]=0.0;} p[1]=1.0; p[n+1]=1.0; prtc(n,p,pi); cpoly(p,pi,n,zr,zi,fail); if(fail[0]) System.out.println("CPOLY HAS FAILED ON THIS EXAMPLE"); prtz(n,zr,zi); pchk(p,pi,n,zr,zi); System.out.println(" EXAMPLE 7d. 48 ZEROS ON A CIRCLE"+ " OF RADIUS 1. CENTERED AT 0,0"); n = 48; for(int i=1; i<=n+1; i++) {p[i]=0.0; pi[i]=0.0;} p[1]=1.0; p[n+1]=1.0; prtc(n,p,pi); cpoly(p,pi,n,zr,zi,fail); if(fail[0]) System.out.println("CPOLY HAS FAILED ON THIS EXAMPLE"); prtz(n,zr,zi); pchk(p,pi,n,zr,zi); System.out.println("end of c419,java test driver for cpoly"); } // c419 test driver for cpoly static void calct(boolean bool[]) { // SUBROUTINE CALCT(BOOL) CALC2890 // COMPUTES T = -P(S)/H(S). // BOOL - LOGICAL, SET TRUE IF H(S) IS ESSENTIALLY ZERO. // 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 // DOUBLE PRECISION 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) // DOUBLE PRECISION HVR,HVI,CMOD // LOGICAL BOOL int n; double hvr[]={0.0}; double hvi[]={0.0}; double atr[]={0.0}; double ati[]={0.0}; n = nn-1; // EVALUATE H(S). polyev(n,sr,si,hr,hi,qhr,qhi,hvr,hvi); bool[0] = cmod(hvr[0],hvi[0]) <= are*5.0*cmod(hr[n],hi[n]); // was 10 if(!bool[0]) { cdivid(-pvr[0],-pvi[0],hvr[0],hvi[0],atr,ati); tr=atr[0]; ti=ati[0]; } else { tr = 0.0; ti = 0.0; } } // end calct static double cauchy(int nn, double pt[], double q[]) { // DOUBLE PRECISION FUNCTION CAUCHY(NN,PT,Q) CAUC3760 // CAUCHY COMPUTES A LOWER BOUND ON THE MODULI OF THE ZEROS OF A // POLYNOMIAL - PT IS THE MODULUS OF THE COEFFICIENTS. // DOUBLE PRECISION Q(NN),PT(NN),X,XM,F,DX,DF, // * DABS,DEXP,DLOG int n; double x, xm, dx, df; double f; pt[nn] = -pt[nn]; // COMPUTE UPPER ESTIMATE OF BOUND. n = nn-1; x = Math.exp( (Math.log(-pt[nn]) - Math.log(pt[1]))/(double)n ); if(pt[n]!=0.0) { // IF NEWTON STEP AT THE ORIGIN IS BETTER, USE IT. xm = -pt[nn]/pt[n]; if(xm0.005) { q[1] = pt[1]; for(int i=2; i<=nn; i++) q[i] = q[i-1]*x+pt[i]; f = q[nn]; df = q[1]; for(int i=2; i<=n; i++) df = df*x+q[i]; dx = f/df; x = x-dx; } // end while return x; } // end cauchy static void cdivid(double ar, double ai, double br, double bi, double cr[], double ci[]) { // SUBROUTINE CDIVID(AR,AI,BR,BI,CR,CI) CDIV4490 // COMPLEX DIVISION C = A/B, AVOIDING OVERFLOW. // DOUBLE PRECISION AR,AI,BR,BI,CR,CI,R,D,T,INFIN,DABS double r, d, t; if(br==0.0 && bi==0.0) { // DIVISION BY ZERO, C = INFINITY. cr[0] = infin; ci[0] = infin; return; } if(Math.abs(br)ai) return ar*Math.sqrt(1.0+(ai/ar)*(ai/ar)); return ar*Math.sqrt(2.0); } // end cmod static void cpoly(double opr[], double opi[], int degree, double zeror[], double zeroi[], boolean fail[]) { // polynomial stored highest coef in opr[1], constant coef in opr[nn] // nn = degree+1 // SUBROUTINE CPOLY(OPR,OPI,DEGREE,ZEROR,ZEROI,FAIL) CPOL 10 // 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 // COMMON/GLOBAL/PR,PI,HR,HI,QPR,QPI,QHR,QHI,SHR,SHI, // * SR,SI,TR,TI,PVR,PVI,ARE,MRE,ETA,INFIN,NN // DOUBLE PRECISION 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) // TO CHANGE THE SIZE OF POLYNOMIALS WHICH CAN BE SOLVED, REPLACE // THE DIMENSION OF THE ARRAYS IN THE COMMON AREA. //EXTERNAL SCALE // DOUBLE PRECISION XX,YY,COSR,SINR,SMALNO,BASE,XXX,ZR,ZI,BND, // * OPR(1),OPI(1),ZEROR(1),ZEROI(1), // * CMOD,SCALE,CAUCHY,DSQRT // LOGICAL FAIL,CONV // INTEGER DEGREE,CNT1,CNT2 boolean conv[]={false}; int idnn2; double xx, yy, cosr, sinr, xxx, bnd; double zr[]={0.0}; double zi[]={0.0}; // INITIALIZATION OF CONSTANTS xx = Math.sqrt(2.0)/2.0; // .70710678; yy = -xx; cosr = Math.cos(Math.PI*94.0/180.0); // -.060756474; sinr = Math.sin(Math.PI*94.0/180.0); // .99756405; fail[0] = false; tr = 0.0; ti = 0.0; sr = 0.0; si = 0.0; nn = degree+1; // ALGORITHM FAILS IF THE LEADING COEFFICIENT IS ZERO. if(opr[1]==0.0 && opi[1]==0.0) { fail[0] = true; return; } // REMOVE THE ZEROS AT THE ORIGIN IF ANY. while(opr[nn]==0.0 && opi[nn]==0.0) // 10 { idnn2 = degree-nn+2; zeror[idnn2] = 0.0; zeroi[idnn2] = 0.0; nn = nn-1; } // end while // MAKE A COPY OF THE COEFFICIENTS. for(int i=1; i<=nn; i++) { pr[i] = opr[i]; pi[i] = opi[i]; shr[i] = cmod(pr[i],pi[i]); } // SCALE THE POLYNOMIAL. bnd = scale(nn,shr,eta,infin,smalno,base); if(bnd!=1.0) { for(int i= 1; i<=nn; i++) { pr[i] = bnd*pr[i]; pi[i] = bnd*pi[i]; } } // START THE ALGORITHM FOR ONE ZERO . l40: while(true) // 40 { if(nn<=2) { // CALCULATE THE FINAL ZERO AND RETURN. double atr[]={0.0}; double ati[]={0.0}; cdivid(-pr[2],-pi[2],pr[1],pi[1],atr,ati); zeror[degree]=atr[0]; zeroi[degree]=ati[0]; return; } // CALCULATE BND, A LOWER BOUND ON THE MODULUS OF THE ZEROS. for(int i=1; i<=nn; i++) shr[i] = cmod(pr[i],pi[i]); bnd = cauchy(nn,shr,shi); // OUTER LOOP TO CONTROL 2 MAJOR PASSES WITH DIFFERENT SEQUENCES // OF SHIFTS. for(int cnt1=1; cnt1<=2; cnt1++) // 100 { // FIRST STAGE CALCULATION, NO SHIFT. noshft(5); // INNER LOOP TO SELECT A SHIFT. for(int cnt2=1; cnt2<=9; cnt2++) // 90 { // 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. fxshft(10*cnt2,zr,zi,conv); if(conv[0]) { // 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[0]; zeroi[idnn2] = zi[0]; nn = nn-1; for(int i=1; i<=nn; i++) { pr[i] = qpr[i]; pi[i] = qpi[i]; } continue l40; // now should go to 40 via end of while } // IF THE ITERATION IS UNSUCCESSFUL ANOTHER SHIFT IS CHOSEN. } // end cnt2 90 // IF 9 SHIFTS FAIL, THE OUTER LOOP IS REPEATED WITH ANOTHER // SEQUENCE OF SHIFTS. } // end cnt1 100 // THE ZEROFINDER HAS FAILED ON TWO MAJOR PASSES. // RETURN EMPTY HANDED. fail[0]=true; return; } // end while 40 } // end cpoly static double errev(int n, double qr[], double qi[], double ms, double mp) { // DOUBLE PRECISION FUNCTION ERREV(NN,QR,QI,MS,MP,ARE,MRE) ERRE3610 // 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 // DOUBLE PRECISION QR(NN),QI(NN),MS,MP,ARE,MRE,E,CMOD double e; e = cmod(qr[1],qi[1])*mre/(are+mre); for(int i=2; i<=nn; i++) e = e*ms+cmod(qr[i],qi[i]); return e*(are+mre)-mp*mre; } // end errev static void fxshft(int l2, double zr[], double zi[], boolean conv[]) { //SUBROUTINE FXSHFT(L2,ZR,ZI,CONV) FXSH1550 // 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 // COMMON/GLOBAL/PR,PI,HR,HI,QPR,QPI,QHR,QHI,SHR,SHI, // * SR,SI,TR,TI,PVR,PVI,ARE,MRE,ETA,INFIN,NN // DOUBLE PRECISION 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) // DOUBLE PRECISION ZR,ZI,OTR,OTI,SVSR,SVSI,CMOD // LOGICAL CONV,TEST,PASD,BOOL int n; boolean bool[] = {true}; boolean test, pasd; double otr, oti, svsr, svsi; n = nn-1; // EVALUATE P AT S. polyev(nn,sr,si,pr,pi,qpr,qpi,pvr,pvi); test = true; pasd = false; // CALCULATE FIRST T = -P(S)/H(S). calct(bool); // MAIN LOOP FOR ONE SECOND STAGE STEP. // DO 50 J = 1,L2 for(int j=1; j<=l2; j++) { otr = tr; oti = ti; // COMPUTE NEXT H POLYNOMIAL AND NEW T. nexth(bool); calct(bool); zr[0] = sr+tr; zi[0] = si+ti; // TEST FOR CONVERGENCE UNLESS STAGE 3 HAS FAILED ONCE OR THIS // IS THE LAST H POLYNOMIAL . if(bool[0] || !test || j== l2) continue; if(cmod(tr-otr,ti-oti) >= 0.5*cmod(zr[0],zi[0])) { pasd = false; continue; } // 40 if(!pasd){ pasd = true; continue;} // 30 // THE WEAK CONVERGENCE TEST HAS BEEN PASSED TWICE, START THE // THIRD STAGE ITERATION, AFTER SAVING THE CURRENT H POLYNOMIAL // AND SHIFT. for(int i=1; i<=n; i++) { shr[i] = hr[i]; shi[i] = hi[i]; } svsr = sr; svsi = si; vrshft(10,zr,zi,conv); if(conv[0]) return; // THE ITERATION FAILED TO CONVERGE. TURN OFF TESTING AND RESTORE // C H,S,PV AND T. test = false; for(int i=1; i<=n; i++) { hr[i] = shr[i]; hi[i] = shi[i]; } sr = svsr; si = svsi; polyev(nn,sr,si,pr,pi,qpr,qpi,pvr,pvi); calct(bool); } // 50 CONTINUE // ATTEMPT AN ITERATION WITH FINAL H POLYNOMIAL FROM SECOND STAGE. vrshft(10,zr,zi,conv); } // end fxshft static void mcon() //(double eta, double infin, double smalno, double base) // stored in COMMON { // SUBROUTINE MCON(ETA,INFINY,SMALNO,BASE) MCON4840 // 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 + ETA IS // GREATER THAN 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 // 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 IEEE // DOUBLE PRECISION ETA,INFINY,SMALNO,BASE int m,n,t; base = 2.0; t = 52; m = 1023; // 127; n = -1023; // -127; eta = Math.pow(base,1-t); infin = base*(1.0-Math.pow(base,-t))*Math.pow(base,m-1); smalno = Math.pow(base,n+3)/Math.pow(base,3); are = eta; mre = 2.0*Math.sqrt(2.0)*eta; } // end mcon static void nexth(boolean bool[]) { // SUBROUTINE NEXTH(BOOL) NEXT3110 // CALCULATES THE NEXT SHIFTED H POLYNOMIAL. // BOOL - LOGICAL, IF .TRUE. H(S) IS ESSENTIALLY ZERO // 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 // DOUBLE PRECISION 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) // DOUBLE PRECISION T1,T2 // LOGICAL BOOL int n; double t1, t2; n = nn-1; if(bool[0]) { // IF H(S) IS ZERO REPLACE H WITH QH. for(int j=2; j<=n; j++) { hr[j] = qhr[j-1]; hi[j] = qhi[j-1]; } hr[1] = 0.0; hi[1] = 0.0; } else { for(int j=2; j<=n; j++) { t1 = qhr[j-1]; t2 = qhi[j-1]; hr[j] = tr*t1-ti*t2+qpr[j]; hi[j] = tr*t2+ti*t1+qpi[j]; } } hr[1] = qpr[1]; hi[1] = qpi[1]; } // end nexth static void noshft(int l1) { // SUBROUTINE NOSHFT(L1) NOSH1130 // COMPUTES THE DERIVATIVE POLYNOMIAL AS THE INITIAL H // POLYNOMIAL AND COMPUTES L1 NO-SHIFT H POLYNOMIALS. // 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 // DOUBLE PRECISION 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) // DOUBLE PRECISION XNI,T1,T2,CMOD int n, nm1, j; double xni, t1, t2; n = nn-1; nm1 = n-1; for(int i=1; i<=n; i++) { xni = nn-i; hr[i] = xni*pr[i]/(double)n; hr[i] = xni*pi[i]/(double)n; } // DO 50 JJ = 1,L1 for(int jj=1; jj<=l1; jj++) { if(cmod(hr[n],hi[n]) <= eta*10.0*cmod(pr[n],pi[n])) // GO TO 30 { // IF THE CONSTANT TERM IS ESSENTIALLY ZERO, SHIFT H COEFFICIENTS. // DO 40 I = 1,NM1 for(int i=1; i<=nm1; i++) { j = nn-i; hr[j] = hr[j-1]; hi[j] = hi[j-1]; } hr[1] = 0.0; hi[1] = 0.0; continue; // GO TO 50 } double atr[]={0.0}; double ati[]={0.0}; cdivid(-pr[nn],-pi[nn],hr[n],hi[n],atr,ati); tr=atr[0]; ti=ati[0]; for(int i=1; i<=nm1; i++) { 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+pi[j]; } hr[1] = pr[1]; hi[1] = pi[1]; } //50 CONTINUE } // end noshft static void polyev(int nn, double sr, double si, double pr[], double pi[], double qr[], double qi[], double pvr[], double pvi[]) { // SUBROUTINE POLYEV(NN,SR,SI,PR,PI,QR,QI,PVR,PVI) POLY3430 // EVALUATES A POLYNOMIAL P AT S BY THE HORNER RECURRENCE // PLACING THE PARTIAL SUMS IN Q AND THE COMPUTED VALUE IN PV. // DOUBLE PRECISION PR(NN),PI(NN),QR(NN),QI(NN), // * SR,SI,PVR,PVI,T double t; qr[1] = pr[1]; qi[1] = pi[1]; pvr[0] = qr[1]; pvi[0] = qi[1]; for(int i=2; i<=nn; i++) { t = pvr[0]*sr-pvi[0]*si+pr[i]; pvi[0] = pvr[0]*si+pvi[0]*sr+pi[i]; pvr[0] = t; qr[i] = pvr[0]; qi[i] = pvi[0]; } } // end polyev static void pchk(double opr[], double opi[], int degree, double zeror[], double zeroi[]) { // SUBROUTINE PCHK(OPR,OPI,DEGREE,ZEROR,ZEROI) // 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. // INTEGER DEGREE // DOUBLE PRECISION OPR(50),OPI(50),ZEROR(50),ZEROI(50) // DOUBLE PRECISION SUMR,SUMI,TEMP DecimalFormat f1 = new DecimalFormat("0.0000000000000000E00"); DecimalFormat f2 = new DecimalFormat("0.00E00"); double sumr, sumi, temp; f1.setPositivePrefix("+"); double val; for(int i=1; i<=degree; i++) { sumr=0.0; sumi=0.0; for(int j=1; j<=degree+1; j++) { temp = opr[j] + sumr*zeror[i]-sumi*zeroi[i]; sumi = opi[j] + sumr*zeroi[i]+sumi*zeror[i]; sumr = temp; } zeror[i] = sumr; zeroi[i] = sumi; } System.out.println("CHECK"); for(int i=1; i<=degree; i++) { val = Math.sqrt(zeror[i]*zeror[i]+zeroi[i]*zeroi[i])/maxc; System.out.println(" "+f1.format(zeror[i])+" "+f1.format(zeroi[i])+ " err="+f2.format(val)); } System.out.println(); } // end pchk static void prtc(int n, double p[], double q[]) { // SUBROUTINE PRTC(N,P,Q) // DOUBLE PRECISION P(50),Q(50) DecimalFormat f1 = new DecimalFormat("0.0000000000000000E00"); f1.setPositivePrefix("+"); double val; System.out.println("COEFFICIENTS degree="+n); maxc = 0.0; for(int i=1; i<=n+1; i++) { val = Math.sqrt(p[i]*p[i]+q[i]*q[i]); if(val>maxc) maxc=val; System.out.println(" "+f1.format(p[i])+" "+f1.format(q[i])); } System.out.println(); } // end prtc static void prtz(int n, double zr[], double zi[]) { // SUBROUTINE PRTZ(N,ZR,ZI) // DOUBLE PRECISION ZR(50),ZI(50) DecimalFormat f1 = new DecimalFormat("0.0000000000000000E00"); DecimalFormat f2 = new DecimalFormat("0.00E00"); f1.setPositivePrefix("+"); System.out.println("ZEROS"); double val; for(int i=1; i<=n; i++) { val = Math.sqrt(zr[i]*zr[i]+zi[i]*zi[i]); System.out.println(" "+f1.format(zr[i])+" "+f1.format(zi[i])+ " abs="+f2.format(val)); } System.out.println(); } // end pchk static double scale(int nn, double pt[], double eta, double infin, double smalno, double base) { // DOUBLE PRECISION FUNCTION SCALE(NN,PT,ETA,INFIN,SMALNO,BASE) SCAL4160 // 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. // DOUBLE PRECISION PT(NN),ETA,INFIN,SMALNO,BASE,HI,LO, // * MAX,MIN,X,SC,DSQRT,DLOG // FIND LARGEST AND SMALLEST MODULI OF COEFFICIENTS. double hi, lo, max, min, x, sc, l; hi = Math.sqrt(infin); lo = smalno/eta; max = 0.0; min = infin; for(int i=1; i<=nn; i++) { x = pt[i]; if(x > max) max = x; if(x != 0.0 && x= lo && max <= hi) return 1.0; x = lo/min; if(x <= 1.0) { sc = 1.0/(Math.sqrt(max)*Math.sqrt(min)); } else { sc = x; if(infin/sc > max) sc = 1.0; } l = Math.log(sc)/Math.log(base) + 0.5; return Math.pow(base,l); } // end scale static void vrshft(int l3, double zr[], double zi[], boolean conv[]) { // SUBROUTINE VRSHFT(L3,ZR,ZI,CONV) VRSH2230 // 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 // COMMON/GLOBAL/PR,PI,HR,HI,QPR,QPI,QHR,QHI,SHR,SHI, // * SR,SI,TR,TI,PVR,PVI,ARE,MRE,ETA,INFIN,NN // DOUBLE PRECISION 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) // DOUBLE PRECISION ZR,ZI,MP,MS,OMP,RELSTP,R1,R2,CMOD,DSQRT,ERREV,TP // LOGICAL CONV,B,BOOL boolean bool[] = {false}; boolean b; double mp, ms, omp, r1, r2, tp; double relstp=.001; // phoney initialization conv[0] = false; b = false; sr = zr[0]; si = zi[0]; // MAIN LOOP FOR STAGE THREE // DO 60 I = 1,L3 for(int i=1; i<=l3; i++) { // EVALUATE P AT S AND TEST FOR CONVERGENCE. polyev(nn,sr,si,pr,pi,qpr,qpi,pvr,pvi); mp = cmod(pvr[0],pvi[0]); omp = mp; ms = cmod(sr,si); if(mp<=10.0*errev(nn,qpr,qpi,ms,mp)) // was 20.0 { // POLYNOMIAL VALUE IS SMALLER IN VALUE THAN A BOUND ON THE ERROR // IN EVALUATING P, TERMINATE THE ITERATION. conv[0] = true; zr[0] = sr; zi[0] = si; return; } if(i == 1) //GO TO 40 { omp = mp; } else if(b || mp=0.05) // GO TO 30 { // EXIT IF POLYNOMIAL VALUE INCREASES SIGNIFICANTLY. if(mp*0.1 > omp) return; omp = mp; } else { // 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 = Math.sqrt(tp); r2 = sr*(1.0+r1)-si*r1; si = sr*r1+si*(1.0+r1); sr = r2; polyev(nn,sr,si,pr,pi,qpr,qpi,pvr,pvi); for(int j=1; j<=5; j++) { calct(bool); nexth(bool); } omp = infin; } // CALCULATE NEXT ITERATE. (and first) calct(bool); nexth(bool); calct(bool); if(!bool[0]) { relstp = cmod(tr,ti)/cmod(sr,si); sr = sr+tr; si = si+ti; } } // 60 CONTINUE main loop on i } // end vrshft public static void main(String args[]) { new c419(); } }