C ALGORITHM 419 COLLECTED ALGORITHMS FROM ACM. C ALGORITHM APPEARED IN COMM. ACM, VOL. 15, NO. 02, C P. 097. C C Added changes from Remark on Algorithm 419 by David H. Withers C CACM (March 1974) Vol 17 No 3 p. 157 PROGRAM CPOLYDR C C DRIVER TO TEST CPOLY C C DIMENSION MULT(50) LOGICAL FLAG(50) REAL*8 FZR(50),FZI(50),BND(50) LOGICAL FAIL DOUBLE PRECISION P(50),PI(50),ZR(50),ZI(50) WRITE(6,100) 100 FORMAT('1EXAMPLE 1. POLYNOMIAL WITH ZEROS 1,2,...,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 DO 10 I=1,11 10 PI(I)=0 CALL PRTC(11,P,PI) CALL CPOLY(P,PI,10,ZR,ZI,FAIL) IF(FAIL) GO TO 95 CALL PRTZ (10,ZR,ZI) 2 WRITE(6,101) 101 FORMAT('1EXAMPLE 2. ZEROS ON IMAGINARY AXIS DEGREE 3.') P(1)=1 P(2)=0 P(3)=-10001.0001D0 P(4)=0 PI(1)=0 PI(2)=-10001.0001D0 PI(3)=0 PI(4)=1 CALL PRTC(4,P,PI) CALL CPOLY(P,PI,3,ZR,ZI,FAIL) IF (FAIL) GO TO 96 CALL PRTZ (3,ZR,ZI) 3 WRITE(6,102) 102 FORMAT('1EXAMPLE 3. ZEROS AT 1+I,1/2*(1+I)....1/(2**-9)*(1+I)') P(1)=1.0 P(2)=-1.998046875 P(3)=0.0 P(4)=.7567065954208374D0 P(5)=-.2002119533717632D0 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 PI(1)=0 PI(2)=P(2) PI(3)=2.658859252929688D0 PI(4)=-7.567065954208374D-1 PI(5)=0 PI(6)=P(6) PI(7)=-7.820779428584501D-4 PI(8)=-P(8) PI(9)=0 PI(10)=P(10) PI(11)=9.094947017729282D-13 CALL PRTC(11,P,PI) CALL CPOLY(P,PI,10,ZR,ZI,FAIL) IF (FAIL) GO TO 97 CALL PRTZ(10,ZR,ZI) 4 WRITE(6,103) 103 FORMAT('1EXAMPLE 4. MULTIPLE ZEROS') 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 PI(1)=0 PI(2)=-10 PI(3)=100 PI(4)=-334 PI(5)=200 PI(6)=1394 PI(7) =-3836 PI(8)=4334 PI(9)=-2352 PI(10)=504 PI(11)=0 CALL PRTC(11,P,PI) CALL CPOLY(P,PI,10,ZR,ZI,FAIL) IF (FAIL) GO TO 98 CALL PRTZ(10,ZR,ZI) 5 WRITE(6,104) 104 FORMAT('1EXAMPLE 5. 12 ZEROS EVENLY DISTRIBUTE ON A CIRCLE OF RADI -US 1. CENTERED AT 0+2I.') 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 PI(1)=0 PI(2)=-24 PI(3)=0 PI(4)=1760 PI(5)=0 PI(6)=-25344 PI(7)=0 PI(8)=101376 PI(9)=0 PI(10)=-112640 PI(11)=0 PI(12)=24576 PI(13)=0 CALL PRTC(13,P,PI) CALL CPOLY(P,PI,12,ZR,ZI,FAIL) IF(FAIL) GO TO 99 CALL PRTZ(12,ZR,ZI) RETURN 95 WRITE(6,105) GO TO 2 96 WRITE(6,105) GO TO 3 97 WRITE(6,105) GO TO 4 98 WRITE(6,105) GO TO 5 99 WRITE(6,105) RETURN 105 FORMAT(//' CPOLY HAS FAILED ON THIS EXAMPLE') END SUBROUTINE PRTC(N,P,Q) DOUBLE PRECISION P(50),Q(50) WRITE(6,10) (P(I),Q(I) ,I=1,N) 10 FORMAT(//' COEFFICIENTS' /50(2D26.16/)) RETURN END SUBROUTINE PRTZ(N,ZR,ZI) DOUBLE PRECISION ZR(50),ZI(50) WRITE(6,10) (ZR(I),ZI(I) ,I=1,N) 10 FORMAT(//' ZEROS'/ 50(2D26.16/)) RETURN END SUBROUTINE CPOLY(OPR,OPI,DEGREE,ZEROR,ZEROI,FAIL) CPOL 10 C C Added changes from Remark on Algorithm 419 by David H. Withers C CACM (March 1974) Vol 17 No 3 p. 157 C C FINDS THE ZEROS OF A COMPLEX POLYNOMIAL. C OPR, OPI - DOUBLE PRECISION VECTORS OF REAL AND C IMAGINARY PARTS OF THE COEFFICIENTS IN C ORDER OF DECREASING POWERS. C DEGREE - INTEGER DEGREE OF POLYNOMIAL. C ZEROR, ZEROI - OUTPUT DOUBLE PRECISION VECTORS OF C REAL AND IMAGINARY PARTS OF THE ZEROS. C FAIL - OUTPUT LOGICAL PARAMETER, TRUE ONLY IF C LEADING COEFFICIENT IS ZERO OR IF CPOLY C HAS FOUND FEWER THAN DEGREE ZEROS. C THE PROGRAM HAS BEEN WRITTEN TO REDUCE THE CHANCE OF OVERFLOW C OCCURRING. IF IT DOES OCCUR, THERE IS STILL A POSSIBILITY THAT C THE ZEROFINDER WILL WORK PROVIDED THE OVERFLOWED QUANTITY IS C REPLACED BY A LARGE NUMBER. C 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) C TO CHANGE THE SIZE OF POLYNOMIALS WHICH CAN BE SOLVED, REPLACE C THE DIMENSION OF THE ARRAYS IN THE COMMON AREA. 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 C INITIALIZATION OF CONSTANTS CALL MCON(ETA,INFIN,SMALNO,BASE) ARE = ETA MRE = 2.0D0*DSQRT(2.0D0)*ETA XX = .70710678 YY = -XX COSR = -.069756474 SINR = .99756405 FAIL = .FALSE. NN = DEGREE+1 C ALGORITHM FAILS IF THE LEADING COEFFICIENT IS ZERO. IF (OPR(1) .NE. 0.0D0 .OR. OPI(1) .NE. 0.0D0) GO TO 10 FAIL = .TRUE. RETURN C REMOVE THE ZEROS AT THE ORIGIN IF ANY. 10 IF (OPR(NN) .NE. 0.0D0 .OR. OPI(NN) .NE. 0.0D0) GO TO 20 IDNN2 = DEGREE-NN+2 ZEROR(IDNN2) = 0.0D0 ZEROI(IDNN2) = 0.0D0 NN = NN-1 GO TO 10 C MAKE A COPY OF THE COEFFICIENTS. 20 DO 30 I = 1,NN PR(I) = OPR(I) PI(I) = OPI(I) SHR(I) = CMOD(PR(I),PI(I)) 30 CONTINUE C SCALE THE POLYNOMIAL. BND = SCALE (NN,SHR,ETA,INFIN,SMALNO,BASE) IF (BND .EQ. 1.0D0) GO TO 40 DO 35 I = 1,NN PR(I) = BND*PR(I) PI(I) = BND*PI(I) 35 CONTINUE C START THE ALGORITHM FOR ONE ZERO . 40 IF (NN.GT. 2) GO TO 50 C CALCULATE THE FINAL ZERO AND RETURN. CALL CDIVID(-PR(2),-PI(2),PR(1),PI(1),ZEROR(DEGREE), * ZEROI(DEGREE)) RETURN C CALCULATE BND, A LOWER BOUND ON THE MODULUS OF THE ZEROS. 50 DO 60 I = 1,NN SHR(I) = CMOD(PR(I),PI(I)) 60 CONTINUE BND = CAUCHY(NN,SHR,SHI) C OUTER LOOP TO CONTROL 2 MAJOR PASSES WITH DIFFERENT SEQUENCES C OF SHIFTS. DO 100 CNT1 = 1,2 C FIRST STAGE CALCULATION, NO SHIFT. CALL NOSHFT(5) C INNER LOOP TO SELECT A SHIFT. DO 90 CNT2 = 1,9 C SHIFT IS CHOSEN WITH MODULUS BND AND AMPLITUDE ROTATED BY C 94 DEGREES FROM THE PREVIOUS SHIFT XXX = COSR*XX-SINR*YY YY = SINR*XX+COSR*YY XX = XXX SR = BND*XX SI = BND*YY C SECOND STAGE CALCULATION, FIXED SHIFT. CALL FXSHFT(10*CNT2,ZR,ZI,CONV) IF (.NOT. CONV) GO TO 80 C THE SECOND STAGE JUMPS DIRECTLY TO THE THIRD STAGE ITERATION. C IF SUCCESSFUL THE ZERO IS STORED AND THE POLYNOMIAL DEFLATED. IDNN2 = DEGREE-NN+2 ZEROR(IDNN2) = ZR ZEROI(IDNN2) = ZI NN = NN-1 DO 70 I = 1,NN PR(I) = QPR(I) PI(I) = QPI(I) 70 CONTINUE GO TO 40 80 CONTINUE C IF THE ITERATION IS UNSUCCESSFUL ANOTHER SHIFT IS CHOSEN. 90 CONTINUE C IF 9 SHIFTS FAIL, THE OUTER LOOP IS REPEATED WITH ANOTHER C SEQUENCE OF SHIFTS. 100 CONTINUE C THE ZEROFINDER HAS FAILED ON TWO MAJOR PASSES. C RETURN EMPTY HANDED. FAIL = .TRUE. RETURN END SUBROUTINE NOSHFT(L1) NOSH1130 C COMPUTES THE DERIVATIVE POLYNOMIAL AS THE INITIAL H C POLYNOMIAL AND COMPUTES L1 NO-SHIFT H POLYNOMIALS. C 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 N = NN-1 NM1 = N-1 DO 10 I = 1,N XNI = NN-I HR(I) = XNI*PR(I)/FLOAT(N) HI(I) = XNI*PI(I)/FLOAT(N) 10 CONTINUE DO 50 JJ = 1,L1 IF (CMOD(HR(N),HI(N)) .LE. ETA*10.0D0*CMOD(PR(N),PI(N))) * GO TO 30 CALL CDIVID(-PR(NN),-PI(NN),HR(N),HI(N),TR,TI) DO 20 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+PI(J) 20 CONTINUE HR(1) = PR(1) HI(1) = PI(1) GO TO 50 C IF THE CONSTANT TERM IS ESSENTIALLY ZERO, SHIFT H COEFFICIENTS. 30 DO 40 I = 1,NM1 J = NN-I HR(J) = HR(J-1) HI(J) = HI(J-1) 40 CONTINUE HR(1) = 0.0D0 HI(1) = 0.0D0 50 CONTINUE RETURN END SUBROUTINE FXSHFT(L2,ZR,ZI,CONV) FXSH1550 C COMPUTES L2 FIXED-SHIFT H POLYNOMIALS AND TESTS FOR C CONVERGENCE. C INITIATES A VARIABLE-SHIFT ITERATION AND RETURNS WITH THE C APPROXIMATE ZERO IF SUCCESSFUL. C L2 - LIMIT OF FIXED SHIFT STEPS C ZR,ZI - APPROXIMATE ZERO IF CONV IS .TRUE. C CONV - LOGICAL INDICATING CONVERGENCE OF STAGE 3 ITERATION C 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 N = NN-1 C EVALUATE P AT S. CALL POLYEV(NN,SR,SI,PR,PI,QPR,QPI,PVR,PVI) TEST = .TRUE. PASD = .FALSE. C CALCULATE FIRST T = -P(S)/H(S). CALL CALCT(BOOL) C MAIN LOOP FOR ONE SECOND STAGE STEP. DO 50 J = 1,L2 OTR = TR OTI = TI C COMPUTE NEXT H POLYNOMIAL AND NEW T. CALL NEXTH(BOOL) CALL CALCT(BOOL) ZR = SR+TR ZI = SI+TI C TEST FOR CONVERGENCE UNLESS STAGE 3 HAS FAILED ONCE OR THIS C IS THE LAST H POLYNOMIAL . IF ( BOOL .OR. .NOT. TEST .OR. J .EQ. L2) GO TO 50 IF (CMOD(TR-OTR,TI-OTI) .GE. .5D0*CMOD(ZR,ZI)) GO TO 40 IF (.NOT. PASD) GO TO 30 C THE WEAK CONVERGENCE TEST HAS BEEN PASSED TWICE, START THE C THIRD STAGE ITERATION, AFTER SAVING THE CURRENT H POLYNOMIAL C AND SHIFT. DO 10 I = 1,N SHR(I) = HR(I) SHI(I) = HI(I) 10 CONTINUE SVSR = SR SVSI = SI CALL VRSHFT(10,ZR,ZI,CONV) IF (CONV) RETURN C THE ITERATION FAILED TO CONVERGE. TURN OFF TESTING AND RESTORE C H,S,PV AND T. TEST = .FALSE. DO 20 I = 1,N HR(I) = SHR(I) HI(I) = SHI(I) 20 CONTINUE SR = SVSR SI = SVSI CALL POLYEV(NN,SR,SI,PR,PI,QPR,QPI,PVR,PVI) CALL CALCT(BOOL) GO TO 50 30 PASD = .TRUE. GO TO 50 40 PASD = .FALSE. 50 CONTINUE C ATTEMPT AN ITERATION WITH FINAL H POLYNOMIAL FROM SECOND STAGE. CALL VRSHFT(10,ZR,ZI,CONV) RETURN END SUBROUTINE VRSHFT(L3,ZR,ZI,CONV) VRSH2230 C CARRIES OUT THE THIRD STAGE ITERATION. C L3 - LIMIT OF STEPS IN STAGE 3. C ZR,ZI - ON ENTRY CONTAINS THE INITIAL ITERATE, IF THE C ITERATION CONVERGES IT CONTAINS THE FINAL ITERATE C ON EXIT. C CONV - .TRUE. IF ITERATION CONVERGES C 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 CONV = .FALSE. B = .FALSE. SR = ZR SI = ZI C MAIN LOOP FOR STAGE THREE DO 60 I = 1,L3 C EVALUATE P AT S AND TEST FOR CONVERGENCE. CALL POLYEV(NN,SR,SI,PR,PI,QPR,QPI,PVR,PVI) MP = CMOD(PVR,PVI) MS = CMOD(SR,SI) IF (MP .GT. 20.0D0*ERREV(NN,QPR,QPI,MS,MP,ARE,MRE)) * GO TO 10 C POLYNOMIAL VALUE IS SMALLER IN VALUE THAN A BOUND ON THE ERROR C IN EVALUATING P, TERMINATE THE ITERATION. CONV = .TRUE. ZR = SR ZI = SI RETURN 10 IF (I .EQ. 1) GO TO 40 IF (B .OR. MP .LT.OMP .OR. RELSTP .GE. .05D0) * GO TO 30 C ITERATION HAS STALLED. PROBABLY A CLUSTER OF ZEROS. DO 5 FIXED C SHIFT STEPS INTO THE CLUSTER TO FORCE ONE ZERO TO DOMINATE. TP = RELSTP B = .TRUE. IF (RELSTP .LT. ETA) TP = ETA R1 = DSQRT(TP) R2 = SR*(1.0D0+R1)-SI*R1 SI = SR*R1+SI*(1.0D0+R1) SR = R2 CALL POLYEV(NN,SR,SI,PR,PI,QPR,QPI,PVR,PVI) DO 20 J = 1,5 CALL CALCT(BOOL) CALL NEXTH(BOOL) 20 CONTINUE OMP = INFIN GO TO 50 C EXIT IF POLYNOMIAL VALUE INCREASES SIGNIFICANTLY. 30 IF (MP*.1D0 .GT. OMP) RETURN 40 OMP = MP C CALCULATE NEXT ITERATE. 50 CALL CALCT(BOOL) CALL NEXTH(BOOL) CALL CALCT(BOOL) IF (BOOL) GO TO 60 RELSTP = CMOD(TR,TI)/CMOD(SR,SI) SR = SR+TR SI = SI+TI 60 CONTINUE RETURN END SUBROUTINE CALCT(BOOL) CALC2890 C COMPUTES T = -P(S)/H(S). C BOOL - LOGICAL, SET TRUE IF H(S) IS ESSENTIALLY ZERO. C 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 N = NN-1 C EVALUATE H(S). CALL POLYEV(N,SR,SI,HR,HI,QHR,QHI,HVR,HVI) BOOL = CMOD(HVR,HVI) .LE. ARE*10.0D0*CMOD(HR(N),HI(N)) IF (BOOL) GO TO 10 CALL CDIVID(-PVR,-PVI,HVR,HVI,TR,TI) RETURN 10 TR = 0.0D0 TI = 0.0D0 RETURN END SUBROUTINE NEXTH(BOOL) NEXT3110 C CALCULATES THE NEXT SHIFTED H POLYNOMIAL. C BOOL - LOGICAL, IF .TRUE. H(S) IS ESSENTIALLY ZERO C 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 N = NN-1 NM1 = N-1 IF (BOOL) GO TO 20 DO 10 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) 10 CONTINUE HR(1) = QPR(1) HI(1) = QPI(1) RETURN C IF H(S) IS ZERO REPLACE H WITH QH. 20 DO 30 J = 2,N HR(J) = QHR(J-1) HI(J) = QHI(J-1) 30 CONTINUE HR(1) = 0.0D0 HI(1) = 0.0D0 RETURN END SUBROUTINE POLYEV(NN,SR,SI,PR,PI,QR,QI,PVR,PVI) POLY3430 C EVALUATES A POLYNOMIAL P AT S BY THE HORNER RECURRENCE C 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 QR(1) = PR(1) QI(1) = PI(1) PVR = QR(1) PVI = QI(1) DO 10 I = 2,NN T = PVR*SR-PVI*SI+PR(I) PVI = PVR*SI+PVI*SR+PI(I) PVR = T QR(I) = PVR QI(I) = PVI 10 CONTINUE RETURN END DOUBLE PRECISION FUNCTION ERREV(NN,QR,QI,MS,MP,ARE,MRE) ERRE3610 C BOUNDS THE ERROR IN EVALUATING THE POLYNOMIAL BY THE HORNER C RECURRENCE. C QR,QI - THE PARTIAL SUMS C MS -MODULUS OF THE POINT C MP -MODULUS OF POLYNOMIAL VALUE C ARE, MRE -ERROR BOUNDS ON COMPLEX ADDITION AND MULTIPLICATION DOUBLE PRECISION QR(NN),QI(NN),MS,MP,ARE,MRE,E,CMOD E = CMOD(QR(1),QI(1))*MRE/(ARE+MRE) DO 10 I = 1,NN E = E*MS+CMOD(QR(I),QI(I)) 10 CONTINUE ERREV = E*(ARE+MRE)-MP*MRE RETURN END DOUBLE PRECISION FUNCTION CAUCHY(NN,PT,Q) CAUC3760 C CAUCHY COMPUTES A LOWER BOUND ON THE MODULI OF THE ZEROS OF A C POLYNOMIAL - PT IS THE MODULUS OF THE COEFFICIENTS. DOUBLE PRECISION Q(NN),PT(NN),X,XM,F,DX,DF, * DABS,DEXP,DLOG PT(NN) = -PT(NN) C COMPUTE UPPER ESTIMATE OF BOUND. N = NN-1 X = DEXP( (DLOG(-PT(NN)) - DLOG(PT(1)))/FLOAT(N) ) IF (PT(N).EQ.0.0D0) GO TO 20 C IF NEWTON STEP AT THE ORIGIN IS BETTER, USE IT. XM = -PT(NN)/PT(N) IF (XM.LT.X) X=XM C CHOP THE INTERVAL (0,X) UNITL F LE 0. 20 XM = X*.1D0 F = PT(1) DO 30 I = 2,NN F = F*XM+PT(I) 30 CONTINUE IF (F.LE. 0.0D0) GO TO 40 X = XM GO TO 20 40 DX = X C DO NEWTON ITERATION UNTIL X CONVERGES TO TWO DECIMAL PLACES. 50 IF (DABS(DX/X) .LE. .005D0) GO TO 70 Q(1) = PT(1) DO 60 I = 2,NN Q(I) = Q(I-1)*X+PT(I) 60 CONTINUE F = Q(NN) DF = Q(1) DO 65 I = 2,N DF = DF*X+Q(I) 65 CONTINUE DX = F/DF X = X-DX GO TO 50 70 CAUCHY = X RETURN END DOUBLE PRECISION FUNCTION SCALE(NN,PT,ETA,INFIN,SMALNO,BASE) SCAL4160 C RETURNS A SCALE FACTOR TO MULTIPLY THE COEFFICIENTS OF THE C POLYNOMIAL. THE SCALING IS DONE TO AVOID OVERFLOW AND TO AVOID C UNDETECTED UNDERFLOW INTERFERING WITH THE CONVERGENCE C CRITERION. THE FACTOR IS A POWER OF THE BASE. C PT - MODULUS OF COEFFICIENTS OF P C ETA,INFIN,SMALNO,BASE - CONSTANTS DESCRIBING THE C FLOATING POINT ARITHMETIC. DOUBLE PRECISION PT(NN),ETA,INFIN,SMALNO,BASE,HI,LO, * MAX,MIN,X,SC,DSQRT,DLOG C FIND LARGEST AND SMALLEST MODULI OF COEFFICIENTS. HI = DSQRT(INFIN) LO = SMALNO/ETA MAX = 0.0D0 MIN = INFIN DO 10 I = 1,NN X = PT(I) IF (X .GT. MAX) MAX = X IF (X .NE. 0.0D0 .AND. X.LT.MIN) MIN = X 10 CONTINUE C SCALE ONLY IF THERE ARE VERY LARGE OR VERY SMALL COMPONENTS. SCALE = 1.0D0 IF (MIN .GE. LO .AND. MAX .LE. HI) RETURN X = LO/MIN IF (X .GT. 1.0D0) GO TO 20 SC = 1.0D0/(DSQRT(MAX)*DSQRT(MIN)) GO TO 30 20 SC = X IF (INFIN/SC .GT. MAX) SC = 1.0D0 30 L = DLOG(SC)/DLOG(BASE) + .500 SCALE = BASE**L RETURN END SUBROUTINE CDIVID(AR,AI,BR,BI,CR,CI) CDIV4490 C COMPLEX DIVISION C = A/B, AVOIDING OVERFLOW. DOUBLE PRECISION AR,AI,BR,BI,CR,CI,R,D,T,INFIN,DABS IF (BR .NE. 0.0D0 .OR. BI .NE. 0.0D0) GO TO 10 C DIVISION BY ZERO, C = INFINITY. CALL MCON (T,INFIN,T,T) CR = INFIN CI = INFIN RETURN 10 IF (DABS(BR) .GE. DABS(BI)) GO TO 20 R = BR/BI D = BI+R*BR CR = (AR*R+AI)/D CI = (AI*R-AR)/D RETURN 20 R = BI/BR D = BR+R*BI CR = (AR+AI*R)/D CI = (AI-AR*R)/D RETURN END DOUBLE PRECISION FUNCTION CMOD(R,I) CMOD4700 C MODULUS OF A COMPLEX NUMBER AVOIDING OVERFLOW. DOUBLE PRECISION R,I,AR,AI,DABS,DSQRT AR = DABS(R) AI = DABS(I) IF (AR .GE. AI) GO TO 10 CMOD = AI*DSQRT(1.0D0+(AR/AI)**2) RETURN 10 IF (AR .LE. AI) GO TO 20 CMOD = AR*DSQRT(1.0D0+(AI/AR)**2) RETURN 20 CMOD = AR*DSQRT(2.0D0) RETURN END SUBROUTINE MCON(ETA,INFINY,SMALNO,BASE) MCON4840 C MCON PROVIDES MACHINE CONSTANTS USED IN VARIOUS PARTS OF THE C PROGRAM. THE USER MAY EITHER SET THEM DIRECTLY OR USE THE C STATEMENTS BELOW TO COMPUTE THEM. THE MEANING OF THE FOUR C CONSTANTS ARE - C ETA THE MAXIMUM RELATIVE REPRESENTATION ERROR C WHICH CAN BE DESCRIBED AS THE SMALLEST POSITIVE C FLOATING-POINT NUMBER SUCH THAT 1.0D0 + ETA IS C GREATER THAN 1.0D0. C INFINY THE LARGEST FLOATING-POINT NUMBER C SMALNO THE SMALLEST POSITIVE FLOATING-POINT NUMBER C BASE THE BASE OF THE FLOATING-POINT NUMBER SYSTEM USED C LET T BE THE NUMBER OF BASE-DIGITS IN EACH FLOATING-POINT C NUMBER(DOUBLE PRECISION). THEN ETA IS EITHER .5*B**(1-T) C OR B**(1-T) DEPENDING ON WHETHER ROUNDING OR TRUNCATION C IS USED. C LET M BE THE LARGEST EXPONENT AND N THE SMALLEST EXPONENT C IN THE NUMBER SYSTEM. THEN INFINY IS (1-BASE**(-T))*BASE**M C AND SMALNO IS BASE**N. C THE VALUES FOR BASE,T,M,N BELOW CORRESPOND TO THE IBM/360. DOUBLE PRECISION ETA,INFINY,SMALNO,BASE INTEGER M,N,T BASE = 16.0D0 T = 14 M = 63 N = -65 ETA = BASE**(1-T) INFINY = BASE*(1.0D0-BASE**(-T))*BASE**(M-1) SMALNO = (BASE**(N+3))/BASE**3 RETURN END