'***************************************************** '* Find all roots of a complex polynomial using * '* Laguerre formulation in complex domain. * '* ------------------------------------------------- * '* SAMPLE RUN: * '* (Find roots of complex polynomial: * '* (2.5 + I) X3 + (7.5 - 12 I) X2 + (-3.75 + 0.4 I * '* ) X + (4.25 - I) ) * '* * '* Error code = 0 * '* * '* Complex roots are: * '* * '* (-1.079156142110099, 4.904068226593107 ) * '* (-.1118457475859128, .6680069609759574 ) * '* ( .2599674028969959,-.3996613904721711 ) * '* * '* ------------------------------------------------- * '* Ref.: From Numath Library By Tuan Dang Trong in * '* Fortran 77 [BIBLI 18]. * '* * '* Basic Release By J-P Moreau, Paris. * '* (www.jpmoreau.fr) * '***************************************************** 'PROGRAM TEST LAGUERRE DEFDBL A-H, O-Z DEFINT I-N OPTION BASE 0 N = 3 DIM A(N + 1, 2), RAC(N + 1, 2), B(N + 1, 2), C(N + 1, 2), D(N + 1, 2) DIM Z(2), ZZ(2), Z1(2), Z2(2) ' Example #1 A(0, 0) = 2.5: A(0, 1) = 1# A(1, 0) = 7.5: A(1, 1) = -12# A(2, 0) = -3.75: A(2, 1) = .4 A(3, 0) = 4.25: A(3, 1) = -1# ' Example #2 (real roots are -3, 1, 2) ' A(0, 0) = 1#: A(0, 1) = 0# ' A(1, 0) = 0#: A(1, 1) = 0# ' A(2, 0) = -7#: A(2, 1) = 0# ' A(3, 0) = 6#: A(3, 1) = 0# ITMAX = 30 'Maximum number of iterations EPS = .00000001#: EPS2 = .000001# 'Minimum, maximum relative error RAC(0, 0) = 0# 'Approximate value of 1st root RAC(0, 1) = 0# N1 = N CLS 'call complex Laguerre procedure GOSUB 2000 'call CLAGUE(N,A,ITMAX,EPS,EPS2,IMP1,RAC,B,C,D) 'print results PRINT PRINT " Error code ="; IMP1 PRINT PRINT " Complex roots are:" PRINT FOR I = 1 TO N1 PRINT "("; RAC(I, 0); ","; RAC(I, 1); ")" NEXT I PRINT END 'of main program '*** Utility procedures for complex numbers *** 'absolute value of z 500 'FUNCTION CABS(ZZ:Complex) xx = ZZ(0): yy = ZZ(1) x1 = ABS(xx): y1 = ABS(yy) IF x1 = 0# THEN W = y1 ELSE IF y1 = 0# THEN W = x1 ELSE IF x1 > y1 THEN W = x1 * SQR(1# + (y1 / x1) * (y1 / x1)) ELSE W = y1 * SQR(1# + (x1 / y1) * (x1 / y1)) END IF END IF END IF CABS = W RETURN 'Z=Z1/Z2 800 'CDIV(Z1,Z2,Z) D = Z2(0) * Z2(0) + Z2(1) * Z2(1) IF D < .000000000001# THEN RETURN Z(0) = (Z1(0) * Z2(0) + Z1(1) * Z2(1)) / D Z(1) = (Z1(1) * Z2(0) - Z1(0) * Z2(1)) / D RETURN 'Z=Z1*Z2 900 'CMUL(Z1,Z2,Z) Z(0) = Z1(0) * Z2(0) - Z1(1) * Z2(1) Z(1) = Z1(0) * Z2(1) + Z1(1) * Z2(0) RETURN 1000 'Z1=CSQRT(Z) ' SQUARE ROOT OF A COMPLEX NUMBER A+I*B = SQRT(X+I*Y) x = Z(0): y = Z(1) IF x = 0# AND y = 0# THEN A = 0# B = 0# ELSE ZZ(0) = Z(0): ZZ(1) = Z(1): GOSUB 500 A = SQR(ABS(x) + CABS * .5#) IF x >= 0# THEN B = y / (A + A) ELSE IF y < 0# THEN B = -A ELSE B = A A = y / (B + B) END IF END IF END IF Z1(0) = A: Z1(1) = B RETURN 2000 'Subroutine CLAGUE(N,A,ITMAX,EPS,EPS2,IMP1,RAC,B,C,D) '================================================================= ' ROOTS OF A COMPLEX COEFFICIENTS POLYNOMIAL ' BY LAGUERRE FORMULA IN COMPLEX DOMAIN '================================================================= ' INPUTS: ' N : ORDER OF POLYNOMIAL ' A : TABLE OF SIZE (0:N) OF COMPLEX COEFFICIENTS STORED IN ' DECREASING ORDER OF X POWERS ' ITMAX: MAXIMUN NUMBER OF ITERATIONS FOR EACH ROOT ' EPS: MINIMAL RELATIVE ERROR ' EPS2: MAXIMAL RELATIVE ERROR ' RAC(0):APPROXIMATE VALUE OF FIRST ROOT (=0. GENERALLY) ' OUTPUTS: ' IMP1: FLAG = 0 CONVERGENCE WITH AT LEAST EPS2 PRECISION ' = 1 NO CONVERGENCE ' RAC: TABLE OF SIZE (0:N) OF FOUND ROOTS ' RAC(I), I=1,N ' WORKING ZONE: ' W: TABLE OF SIZE (0:3*N), here divided into B, C, D. ' ---------------------------------------------------------------- ' REFERENCE: ' E.DURAND. SOLUTIONS NUMERIQUES DES EQUATIONS ALGEBRIQUES ' TOME I, MASSON & CIE PAGES 269-270 '================================================================= DIM XK(2), XR(2), F(2), FP(2), FS(2), H(2), DEN(2), SQ(2), D2(2) DIM TMP1(2), TMP2(2) IMP1 = 0 B(0, 0) = A(0, 0): B(0, 1) = A(0, 1) C(0, 0) = B(0, 0): C(0, 1) = B(0, 1) D(0, 0) = C(0, 0): D(0, 1) = C(0, 1) IF N = 1 THEN 'RAC(1)=-A(1)/A(0) Z1(0) = A(1, 0): Z1(1) = A(1, 1): Z2(0) = A(0, 0): Z2(1) = A(0, 1) GOSUB 800: RAC(1, 0) = -Z(0): RAC(1, 1) = -Z(1) RETURN END IF XK(0) = RAC(0, 0): XK(1) = RAC(0, 1) 1 IK = 0 EPS1 = EPS IT = 0 2 FOR I = 1 TO N 'B(I)=A(I)+XK*B(I-1) 'Z=XK*B(I-1) Z1(0) = XK(0): Z1(1) = XK(1) Z2(0) = B(I - 1, 0): Z2(1) = B(I - 1, 1) GOSUB 900 'B(I)=A(I)+Z B(I, 0) = A(I, 0) + Z(0) B(I, 1) = A(I, 1) + Z(1) IF I <= N - 1 THEN 'C(I)=B(I)+XK*C(I-1) 'Z=XK*C(I-1) Z2(0) = C(I - 1, 0): Z2(1) = C(I - 1, 1) GOSUB 900 'C(I)=B(I)+Z C(I, 0) = B(I, 0) + Z(0) C(I, 1) = B(I, 1) + Z(1) END IF IF I <= N - 2 THEN 'D(I)=C(I)+XK*D(I-1) 'Z=XK*D(I-1) Z2(0) = D(I - 1, 0): Z2(1) = D(I - 1, 1) GOSUB 900 'D(I)=C(I)+Z D(I, 0) = C(I, 0) + Z(0) D(I, 1) = C(I, 1) + Z(1) END IF NEXT I F(0) = B(N, 0): F(1) = B(N, 1) FP(0) = C(N - 1, 0): FP(1) = C(N - 1, 1) FS(0) = D(N - 2, 0): FS(1) = D(N - 2, 1) 'H=((N-1)*FP)^2-N*(N-1)*F*FS Z1(0) = 1# * (N - 1): Z1(1) = 0# 'Z=Z1*FP Z2(0) = FP(0): Z2(1) = FP(1): GOSUB 900 'TMP1=Z*Z Z1(0) = Z(0): Z1(1) = Z(1): Z2(0) = Z(0): Z2(1) = Z(1) GOSUB 900: TMP1(0) = Z(0): TMP1(1) = Z(1) 'Now TMP1=((N-1)*FP)^2 Z1(0) = 1# * N * (N - 1): Z1(1) = 0# 'Z=Z1*F Z2(0) = F(0): Z2(1) = F(1): GOSUB 900 'CMUL(Z,FS,Z) Z1(0) = Z(0): Z1(1) = Z(1): Z2(0) = FS(0): Z2(1) = FS(1) GOSUB 900: TMP2(0) = Z(0): TMP2(1) = Z(1) 'Now TMP2=N*(N-1)*F*FS 'H=TMP1-TMP2 H(0) = TMP1(0) - TMP2(0) H(1) = TMP1(1) - TMP2(1) ZZ(0) = H(0): ZZ(1) = H(1): GOSUB 500 'Now CABS=CABS(H) IF CABS = 0# THEN FOR I = N TO 1 STEP -1 'RAC(I)=-A(1)/(N*A(0)) Z1(0) = 1# * N: Z1(1) = 0# 'Z=Z1*A(0) Z2(0) = A(0, 0): Z2(1) = A(0, 1): GOSUB 900 'RAC(I)=A(1)/Z Z1(0) = A(1, 0): Z1(1) = A(1, 1): Z2(0) = Z(0): Z2(1) = Z(1) GOSUB 800 RAC(I, 0) = -Z(0) RAC(I, 1) = -Z(1) NEXT I RETURN END IF ZZ(0) = FP(0): ZZ(1) = FP(1): GOSUB 500 'Now CABS=CABS(FP) IF CABS = 0# THEN 'DEN=CSQRT(H) Z(0) = H(0): Z(1) = H(1): GOSUB 1000 DEN(0) = Z1(0): DEN(1) = Z1(1) ELSE 'SQ=CSQRT(H) Z(0) = H(0): Z(1) = H(1): GOSUB 1000 SQ(0) = Z1(0): SQ(1) = Z1(1) 'DEN=FP+SQ DEN(0) = FP(0) + SQ(0) DEN(1) = FP(1) + SQ(1) 'D2=FP-SQ D2(0) = FP(0) - SQ(0) D2(1) = FP(1) - SQ(1) ZZ(0) = DEN(0): ZZ(1) = DEN(1): GOSUB 500 TMP = CABS 'Now TMP=CABS(DEN) ZZ(0) = D2(0): ZZ(1) = D2(1): GOSUB 500 'Now CABS=CABS(D2) IF TMP < CABS THEN DEN(0) = D2(0) DEN(1) = D2(1) END IF END IF ZZ(0) = DEN(0): ZZ(1) = DEN(1): GOSUB 500 IF CABS = 0# THEN XK(0) = XK(0) + .1 XK(1) = XK(1) + .1 GOTO 2 END IF IK = IK + 1 'XR=XK-N*F/DEN Z1(0) = 1# * N: Z1(1) = 0# 'Z=Z1*F Z2(0) = F(0): Z2(1) = F(1): GOSUB 900 'Z=Z/DEN Z1(0) = Z(0): Z1(1) = Z(1): Z2(0) = DEN(0): Z2(1) = DEN(1) GOSUB 800 XR(0) = XK(0) - Z(0) XR(1) = XK(1) - Z(1) 'ZZ=XR-XK ZZ(0) = XR(0) - XK(0) ZZ(1) = XR(1) - XK(1) GOSUB 500: TEST = CABS 'Now TEST=CABS(XR-XK) XK(0) = XR(0): XK(1) = XR(1) ZZ(0) = XR(0): ZZ(1) = XR(1): GOSUB 500 'Now CABS=CABS(XR) IF TEST < EPS1 * CABS THEN GOTO 3 ' For debug only ' PRINT IK, H(0), H(1), DEN(0), DEN(1) ' PRINT XK(0), XK(1), TEST IF IK <= ITMAX THEN GOTO 2 EPS1 = EPS1 * 10# IK = 0 IF EPS1 < EPS2 THEN GOTO 2 IMP1 = 1 RETURN 3 RAC(N, 0) = XR(0) RAC(N, 1) = XR(1) N = N - 1 IF N = 1 THEN GOTO 4 IF N <= 0 THEN RETURN FOR I = 0 TO N A(I, 0) = B(I, 0) A(I, 1) = B(I, 1) NEXT I GOTO 1 4 'RAC(N)=-B(1)/B(0) Z1(0) = B(1, 0): Z1(1) = B(1, 1): Z2(0) = B(0, 0): Z2(1) = B(0, 1) GOSUB 800 RAC(N, 0) = -Z(0) RAC(N, 1) = -Z(1) RETURN 'end of file tclague.bas