C C ________________________________________________________ C | | C | USE JENKENS-TRAUB ALGORITHM TO COMPUTE ZEROS OF A MONIC| C | POLYNOMIAL P(1) + P(2)X +...+ P(ND)X**(ND-1) + X**ND | C | | C | INPUT: | C | | C | P --THE FIRST ND COEFFICIENTS OF POLYNOMIAL| C | (COMPLEX DOUBLE PRECISION ARRAY) | C | | C | ND --DEGREE OF POLYNOMIAL | C | | C | W --WORK ARRAY (REQUIRED LENGTH DEPENDS ON | C | THE DIFFICULTY IN COMPUTING ZEROS. IN | C | MANY CASES, BETWEEN 4ND+32 AND 6ND | C | COMPLEX DOUBLE PRECISION ENTRIES ARE | C | SUFFICIENT) | C | | C | OUTPUT: | C | | C | Z --ZEROS (COMPLEX DOUBLE PRECISION ARRAY) | C | | C | W(1) --NUMBER OF COMPUTED ZEROS | C | | C | NOTE: | C | | C | THIS ALGORITHM CAN GENERATE UNDERFLOWS SINCE THE | C | FUNCTION VALUES TEND TO ZERO AS THE ITERATIONS | C | APPROACH A ROOT. THUS THE APPROPRIATE COMPILER | C | OPTIONS MUST BE INVOKED SO THAT AN UNDERFLOW IS | C | REPLACED BY ZERO AND THE EXECUTION CONTINUES. FOR | C | IBM FORTRAN, TRY THE FOLLOWING STATEMENT: | C | CALL ERRSET(208,256,-1,-1,0) | C | | C | BUILTIN FUNCTIONS: CDABS,CDEXP,DACOS,DCMPLX,DLOG | C | PACKAGE FUNCTIONS: DMAG,DREAL,DIMAG | C | PACKAGE SUBROUTINES: CC,CEF,ERR,FF,JT,QAD,RAD,RFN | C |________________________________________________________| C SUBROUTINE CZERO(Z,P,ND,W) COMPLEX*16 P(1),W(1),Z(1),PI,R,S,U,V,X DOUBLE PRECISION A,B,C,D,E,DMAG,DREAL,S1,T,T1,T2,V1,Y,Y0 INTEGER F,G,H,I,J,K,L,M,N,NC,ND,NZ,N1,O,O1,O2,O3,O4,Q N = ND W(1) = N DO 10 I = 1,N 10 Z(I) = P(I) 20 IF ( N .GT. 1 ) GOTO 25 Z(1) = -P(1) RETURN 25 IF ( N .GT. 2 ) GOTO 30 CALL QAD(Z(2),Z(1)) RETURN C -------------------------------- C |*** REMOVE ZEROS AT ORIGIN ***| C -------------------------------- 30 IF ( DMAG(Z(1)) .NE. 0.D0 ) GOTO 50 DO 40 I = 2,N 40 Z(I-1) = Z(I) Z(N) = (0.D0,0.D0) N = N - 1 GOTO 20 50 NZ = N E = 1.D0 C --------------------------------- C |*** COMPUTE MACHINE EPSILON ***| C --------------------------------- 60 E = .5D0*E A = 1.D0 + E IF( A .GT. 1.D0 ) GOTO 60 E = E + E M = 5 PI = DCMPLX(0.D0,2.D0*DACOS(-1.D0)) C ---------------------------------------------------- C |*** SOME COMPILERS USE DARCOS INSTEAD OF DACOS ***| C ---------------------------------------------------- X = (1.D0,0.D0) Y = 65536. Y0 = 1./Y 70 N1 = N + 1 O1 = N1 + N O2 = O1 + N O3 = O2 + N1/2 O4 = O3 + N1/2 O = O4 - 1 CALL CC(Z,W(O1),W(O2),N,Y,Y0) G = O4 + N1 NC = 32*N K = 10 L = 8 Q = 0 J = 0 F = 0 80 F = F + 1 Q = Q + J IF ( F .LT. 6 ) GOTO 90 IF ( L .GT. NC ) GOTO 370 90 L = L + L C -------------------------------------------------------------- C |*** ESTIMATE RADIUS OF SMALLEST CIRCLE CONTAINING A ZERO ***| C -------------------------------------------------------------- CALL RAD(Z,W,W(N1),K,N,A) CALL CEF(W(O4),Z,W(O1),W(O2),W(O3),A,X,N,N1,T,T1,Y,Y0) T = N*T*E J = 1 I = N1 H = I 100 I = I/2 J = J + J IF ( I .GT. 0 ) GOTO 100 IF ( H+H .EQ. J ) GOTO 120 H = J J = H + O DO 110 I = G,J 110 W(I) = (0.D0,0.D0) 120 J = L + O4 CALL FF(W(O4),W(J),H,L) C = DMAG(W(O4)) J = 0 H = L - 1 C ---------------------------------- C |*** DETERMINE STARTING GUESS ***| C ---------------------------------- DO 130 I=1,H D = DMAG(W(O4+I)) IF ( D .GT. C ) GOTO 130 J = I C = D 130 CONTINUE S = A*X*CDEXP(J*PI/L) S1 = CDABS(S) IF ( S1 .LT. 1. ) S1 = 1. X = X*CDEXP(PI/(4*L)) X = X/CDABS(X) C ---------------------------------- C |*** JENKINS-TRAUB ITERATIONS ***| C ---------------------------------- DO 140 I=1,M 140 CALL JT(Z,W,N,V,S,S1) C = 0.D0 J = 0 150 CALL JT(Z,W,N,V,S,S1) V1 = DMAG(V) IF ( V1 .EQ. 0. ) GOTO 240 IF ( S1 .GT. 1. ) GOTO 170 U = (1.D0,0.D0) I = N 160 I = I - 1 U = W(I) + U*S IF ( I .GT. 1 ) GOTO 160 IF ( DMAG(U) .EQ. 0.D0 ) GOTO 80 GOTO 190 170 U = W(1) R = 1./S DO 180 I = 2,N 180 U = W(I) + U*R IF ( DMAG(U) .EQ. 0.D0 ) GOTO 80 U = U*R 190 R = S - V/U J = J + 1 D = DMAG(R-S) S = R IF ( S1 .LE. A ) GOTO 200 CALL ERR(Z,N,S1,T2) T2 = N*T2*E IF ( V1 .LE. T2 ) GOTO 240 GOTO 210 200 IF ( N*DLOG(T1/S1)+DLOG(T/V1) .GE. 0. ) GOTO 230 210 S1 = CDABS(S) IF ( S1 .LT. 1. ) S1 = 1. IF ( J .LT. 3 ) GOTO 220 IF ( .5D0*D .GT. B+C ) GOTO 80 IF ( J .LT. M ) GOTO 220 IF ( 4.D0*D .GT. B+C ) GOTO 80 220 B = C C = D GOTO 150 230 CALL ERR(Z,N,S1,T2) T2 = N*E*T2 IF ( V1 .GT. T2 ) GOTO 210 240 R = Z(N) + S I = N 250 I = I - 1 V = Z(I) + R*S Z(I) = R R = V IF ( I .GT. 1 ) GOTO 250 Z(N) = S Q = Q + J N = N - 1 IF ( N .GT. 2 ) GOTO 70 CALL QAD(Z(2),Z(1)) NC = ND N = 1 260 L = N N = ND M = N + N J = 1 K = N - 1 DO 270 I = 2,K W(J) = J*P(I) W(J+N) = J*I*P(I+1) W(J+M) = DMAG(P(J)) 270 J = I W(K) = K*P(N) W(N) = N W(M-1) = N*K W(M) = 0.D0 W(M+K) = DMAG(P(K)) W(M+N) = DMAG(P(N)) O1 = M + 1 N1 = N + 1 C -------------------------- C |*** REFINE THE ZEROS ***| C -------------------------- DO 340 K = L,NZ S = Z(K) DO 280 J = 1,3 CALL RFN(P,W,W(N1),W(O1),N,S,V,C,T) T = N*T*E 280 IF ( C .LE. T ) GOTO 330 A = CDABS(S) IF ( A .GT. 1.D0 ) GOTO 300 V = (1.D0,0.D0) T = 1. I = N 290 V = V*S + P(I) T = T*A + DREAL(W(I+M)) I = I - 1 IF ( I .GT. 0 ) GOTO 290 GOTO 320 300 R = 1./S A = 1./A V = (0.D0,0.D0) T = 0. DO 310 I = 1,N V = V*R + P(I) 310 T = T*A + DREAL(W(I+M)) V = V*R + DCMPLX(1.D0,0.D0) T = T*A + 1. 320 A = DMAG(V) IF ( A .GT. N*T*E ) GOTO 340 330 Z(K) = S 340 CONTINUE DO 350 I = L,N 350 W(N1-I) = Z(I) L = N1 - L DO 360 I = 1,L 360 Z(I) = W(I) W(1) = NC RETURN 370 NC = ND - N + 1 WRITE(6,*) 'SINCE THE STOPPING CRITERION HAS NOT BEEN SATISFIED' WRITE(6,*) 'AFTER',Q,'ITERATIONS, WE STOP WHILE COMPUTING' WRITE(6,*) 'ROOT NUMBER',NC Z(N) = S GOTO 260 END C% C ------------------------------ C |*** ROOTS OF A QUADRATIC ***| C ------------------------------ SUBROUTINE QAD(B,C) COMPLEX*16 B,C,R,S DOUBLE PRECISION DMAG,T T = DMAG(B) R = DCMPLX(4.D0,0.D0)*C IF ( T .GT. DSQRT(DMAG(R)) ) GOTO 10 R = CDSQRT(B*B-R) GOTO 20 10 S = 1.D0/T R = T*CDSQRT((S*B)**2-R*S*S) 20 S = R + B R = R - B T = DMAG(R) IF ( T .GE. DMAG(S) ) GOTO 30 C = -(C+C)/S B = DCMPLX(-.5D0,0.D0)*S RETURN 30 IF ( T .EQ. 0.D0 ) RETURN C = (C+C)/R B = DCMPLX(.5D0,0.D0)*R RETURN END C% C --------------------------------- C |*** JENKINS-TRAUB ITERATION ***| C --------------------------------- SUBROUTINE JT(P,Q,N,V,S,S1) COMPLEX*16 P(1),Q(1),R,S,T,U,V,W DOUBLE PRECISION DMAG,S1 INTEGER I,M,N IF ( S1 .GT. 1. ) GOTO 70 U = (0.D0,0.D0) V = (1.D0,0.D0) I = N 10 V = P(I) + V*S U = Q(I) + U*S I = I - 1 IF ( I .GT. 0 ) GOTO 10 IF ( DMAG(U) .NE. 0.D0 ) GOTO 50 M = N 20 R = Q(M) Q(M) = (0.D0,0.D0) I = M 30 I = I - 1 U = Q(I) + R*S Q(I) = R R = U IF ( I .GT. 1 ) GOTO 30 U = (0.D0,0.D0) I = M 40 U = Q(I) + U*S I = I - 1 IF ( I .GT. 0 ) GOTO 40 M = M - 1 IF ( DMAG(U) .EQ. 0.D0 ) GOTO 20 50 R = V/U T = P(N) - R*Q(N) + S Q(N) = (1.D0,0.D0) I = N 60 I = I - 1 U = P(I) - R*Q(I) + T*S Q(I) = T T = U IF ( I .GT. 1 ) GOTO 60 RETURN 70 W = 1./S U = Q(1) V = P(1) DO 80 I = 2,N U = U*W + Q(I) 80 V = V*W + P(I) U = U*W V = V*W + 1. IF ( DMAG(U) .NE. 0. ) GOTO 120 M = N 90 R = Q(M) Q(M) = (0.D0,0.D0) I = M 100 I = I - 1 U = Q(I) + R*S Q(I) = R R = U IF ( I .GT. 1 ) GOTO 100 U = (0.D0,0.D0) M = M - 1 DO 110 I = 1,N 110 U = U*W + Q(I) U = U*W IF ( DMAG(U) .EQ. 0. ) GOTO 90 120 R = V/U U = P(1) - R*Q(1) Q(1) = -W*U DO 130 I = 2,N U = P(I) - R*Q(I) + U*W 130 Q(I) = -W*U Q(N) = (1.D0,0.D0) RETURN END C% C ----------------------- C |*** A NEWTON-STEP ***| C ----------------------- SUBROUTINE RFN(P,Q,R,Z,N,S,V,B,C) COMPLEX*16 P(1),Q(1),R(1),Z(1),M,S,T,U,V,W DOUBLE PRECISION A,B,C,DMAG,DREAL INTEGER I,N A = CDABS(S) IF ( A .GT. 1.D0 ) GOTO 20 V = (1.D0,0.D0) U = (0.D0,0.D0) W = (0.D0,0.D0) C = 1. I = N 10 V = V*S + P(I) U = U*S + Q(I) W = W*S + R(I) C = C*A + DREAL(Z(I)) I = I - 1 IF ( I .GT. 0 ) GOTO 10 GOTO 40 20 T = DCMPLX(1.D0,0.D0)/S A = 1./A V = P(1) U = Q(1) W = R(1) C = DREAL(Z(1)) DO 30 I = 2,N V = V*T + P(I) U = U*T + Q(I) W = W*T + R(I) 30 C = C*A + DREAL(Z(I)) V = V*T + DCMPLX(1.D0,0.D0) U = U*T W = W*T C = C*A + 1. 40 B = DMAG(V) IF ( B .EQ. 0.D0 ) RETURN A = DMAG(U) IF ( A .EQ. 0.D0 ) RETURN T = V/U A = DMAG(W) IF ( A .EQ. 0.D0 ) GOTO 50 U = U/W M = U - T A = DMAG(M) IF ( A .EQ. 0.D0 ) RETURN M = U/M S = S - M*T RETURN 50 S = S - T RETURN END C% C -------------------------------------------------------------- C |*** ESTIMATE RADIUS OF SMALLEST CIRCLE CONTAINING A ZERO ***| C -------------------------------------------------------------- SUBROUTINE RAD(P,Q,R,K,N,A) COMPLEX*16 P(1),Q(1),R(1),T DOUBLE PRECISION A,C,D,F,DMAG,S INTEGER I,J,K,L,M,N SAVE J,C,F DATA J/1/ IF ( J .EQ. K ) GOTO 20 J = 1 L = 1 DO 10 I = 2,N Q(L) = L*P(I)/N 10 L = I Q(N) = 1 C = N*CDABS(P(1)) F = 1.D0/N A = CDABS(P(1))**(1.D0/DFLOAT(N)) GOTO 40 20 DO 30 I = 1,N 30 Q(I) = R(I) 40 L = J + K 50 IF ( DMAG(Q(1)) .EQ. 0.D0 ) GOTO 90 J = J+1 T = P(1)/Q(1) DO 60 I = 2,N 60 Q(I-1) = P(I) - T*Q(I) S = 1.D0/DFLOAT(J) F = F**(1.D0-S)*CDABS(T)**S D = F*(C/CDABS(Q(1)))**S IF ( D .LT. A ) A = D 70 IF ( J .LT. L ) GOTO 50 K = J DO 80 I = 1,N 80 R(I) = Q(I) RETURN 90 M = J 100 J = J + 1 DO 110 I = 2,N 110 Q(I-1) = Q(I) Q(N) = (0.D0,0.D0) IF ( DMAG(Q(1)) .EQ. 0.D0 ) GOTO 100 S = 1.D0/J F = F**(DFLOAT(M)*S) D = F*(C/CDABS(Q(1)))**S IF ( D .LT. A ) A = D T = P(1)/Q(1) DO 120 I = 2,N 120 Q(I-1) = P(I) - T*Q(I) Q(N) = (1.D0,0.D0) F = F*CDABS(T)**S GOTO 70 END C% C -------------------------------- C |*** FAST FOURIER TRANSFORM ***| C -------------------------------- SUBROUTINE FF(A,B,N,L) COMPLEX*16 A(1),B(1),S,T,W DOUBLE PRECISION P INTEGER I,IL,I1,I2,J,K,L,L1,L2,M,M0,M1,N M = N IF ( L .GE. M ) GOTO 20 M = M - 1 DO 10 J = L,M,L DO 10 I = 1,L 10 A(I) = A(I) + A(I+J) M = L 20 P = DACOS(-1.D0) C ---------------------------------------------------- C |*** SOME COMPILERS USE DARCOS INSTEAD OF DACOS ***| C ---------------------------------------------------- W = CDEXP(M*DCMPLX(0.D0,P)/L) L1 = L - 1 L2 = L/2 IF ( M .EQ. L ) GOTO 40 DO 30 J = M,L1,M DO 30 I = 1,M 30 A(I+J) = A(I) 40 M0 = M M = M/2 M1 = M - 1 IF ( M .EQ. 0 ) RETURN S = 1 T = -1 I1 = 0 DO 60 J = 1,L,M0 IL = J + M1 I2 = I1 + L2 DO 50 I = J,IL K = I + M B(I+I1) = A(I) + S*A(K) 50 B(I+I2) = A(I) + T*A(K) S = W*S T = W*T I1 = I1 - M 60 CONTINUE DO 70 I = 1,L 70 A(I) = B(I) W = CDSQRT(W) GOTO 40 END C% DOUBLE PRECISION FUNCTION DMAG(Z) COMPLEX*16 Z DOUBLE PRECISION DIMAG,DREAL DMAG = DABS(DREAL(Z)) + DABS(DIMAG(Z)) RETURN END DOUBLE PRECISION FUNCTION DIMAG(A) DOUBLE PRECISION C(2) COMPLEX*16 A,B EQUIVALENCE (B,C) B = A DIMAG = C(2) RETURN END DOUBLE PRECISION FUNCTION DREAL(A) DOUBLE PRECISION C(2) COMPLEX*16 A,B EQUIVALENCE (B,C) B = A DREAL = C(1) RETURN END C% SUBROUTINE CC(Z,F,E,N,Y,Y0) COMPLEX*16 F(1),Z(1),X DOUBLE PRECISION DMAG,E(1),S,T,Y,Y0 INTEGER I,N DO 40 I = 1,N S = 0. X = Z(I) T = DMAG(X) IF ( T .GT. 1. ) GOTO 20 IF ( T .GT. Y0 ) GOTO 30 IF ( T .EQ. 0. ) GOTO 30 10 T = T*Y X = X*Y S = S - 1. IF ( T .LE. Y0 ) GOTO 10 GOTO 30 20 T = T*Y0 X = X*Y0 S = S + 1. IF ( T .GT. 1. ) GOTO 20 30 E(I) = S F(I) = X 40 CONTINUE RETURN END C% SUBROUTINE CEF(W,Z,F,E,D,R,X,N,N1,T,T1,Y,Y0) COMPLEX*16 F(1),W(1),Z(1),C,X,P DOUBLE PRECISION D(1),DMAG,E(1),A,B,Q,R,S,T,T1,Y,Y0 INTEGER I,N,N1 P = 1. Q = 0. S = E(1) C = X*R IF ( R .GT. 1. ) GOTO 50 T = 1. T1 = 1. DO 20 I = 1,N A = E(I) + Q IF ( A .GT. S ) S = A W(I) = P*F(I) D(I) = A T = R*T + DMAG(Z(N1-I)) P = P*C IF ( DMAG(P) .GT. Y0 ) GOTO 20 10 P = P*Y Q = Q - 1. IF ( DMAG(P) .LE. Y0 ) GOTO 10 20 CONTINUE 30 IF ( Q .GT. S ) S = Q DO 40 I = 1,N A = D(I) - S 40 W(I) = W(I)*Y**A A = Q - S W(N1) = P*Y**A RETURN 50 T1 = R B = 1./R T = 0. DO 70 I = 1,N A = E(I) + Q IF ( A .GT. S ) S = A W(I) = P*F(I) D(I) = A T = B*T + DMAG(Z(I)) P = P*C IF ( DMAG(P) .LE. 1. ) GOTO 70 60 P = P*Y0 Q = Q + 1. IF ( DMAG(P) .GT. 1. ) GOTO 60 70 CONTINUE T = B*T + 1. GOTO 30 END C% SUBROUTINE ERR(Z,N,S,T) COMPLEX*16 Z(1) DOUBLE PRECISION DMAG,R,S,T INTEGER I,N IF ( S .GT. 1. ) GOTO 20 T = 1. I = N 10 T = T*S + DMAG(Z(I)) I = I - 1 IF ( I .GT. 0 ) GOTO 10 RETURN 20 R = 1./S T = DMAG(Z(1)) DO 30 I = 2,N 30 T = T*R + DMAG(Z(I)) T = T*R + 1. RETURN END