'****************************************************
'*   Find all roots of a complex polynomial using   *
'*   Newton's iterative formulation.                *
'* ------------------------------------------------ *
'* 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:                              *
'*  0.25996740  -0.39966139                         *
'* -0.11184574   0.66800696                         *
'* -1.07915614   4.90406822                         *
'*                                                  *
'* ------------------------------------------------ *
'* Ref.: From Numath Library By Tuan Dang Trong in  *
'*       Fortran 77 [BIBLI 18].                     *
'*                                                  *
'*              Basic Release By J-P Moreau, Paris. *
'*                       (www.jpmoreau.fr)          *
'****************************************************
DEFDBL A-H, O-Z
DEFINT I-N

  N = 3  'order of polynomial

  DIM A(N + 1), B(N + 1), X(N), Y(N)
' A: real parts of complex coefficients
' B: imaginary parts of complex coefficients
'    (stored in decreasing powers)
' X: real parts of roots
' Y: imaginary parts of roots */
  DIM W1(N + 1), W2(N + 1), W3(N + 1), W4(N + 1)
'     working zones

'define vectors A, B of coefficients
  A(1) = 2.5: B(1) = 1#
  A(2) = 7.5: B(2) = -12#
  A(3) = -3.75: B(3) = .4
  A(4) = 4.25: B(4) = -1#

' Example #2 (real roots are -3, 1, 2)
' A(1)= 1.   B(1)=0.0
' A(2)= 0.   B(2)=0.0
' A(3)=-7.   B(3)=0.0
' A(4)= 6.   B(4)=0.0  */

' call main function
  GOSUB 1000 'call CNEWTON(N,A,B,RR,RI,IER,W1,W2,W3,W4)

  F$ = "##.########"

' print results
  CLS
  PRINT
  PRINT "  Error code = "; IER
  PRINT
  PRINT "  Complex roots are:"
  
  FOR I = 1 TO N
    PRINT " ";
    PRINT USING F$; X(I);
    PRINT "  ";
    PRINT USING F$; Y(I)
  NEXT I
  PRINT

END

1000 'CNEWTON (N, A, B, X, Y, IER, W1, W2, W3, W4)
'---------------------------------------------------------------------
'     CALCULATE ALL THE ROOTS OF A COMPLEX POLYNOMIAL USING
'     NEWTON'S ITERATIVE FORMULATION
'
'     INPUTS:
'     N        ORDER OF POLYNOMIAL
'     A,B      TABLES OF DIMENSION N+1, RESPECTIVELY REAL,
'              IMAGINARY PART OF POLYNOMIAL COMPLEX COEFFICIENTS
'              C = A+I*B, STORED IN DECREASING POWERS.
'     OUTPUTS:
'     X,Y      TABLES OF DIMENSION N, RESPECTIVELY REAL,
'              IMAGINARY PART OF COMPLEX ROOTS FOUND,  Z = X+I*Y
'     IER      ERROR CODE = 0, CONVERGENCE OK
'                         = I, NO CONVERGENCE FOR I-TH ROOT
'     WORKING ZONES:
'     G,H,T,F  TABLES OF DIMENSION N
'----------------------------------------------------------------------
'     REFERENCE:
'     E.DURAND -SOLUTIONS NUMERIQUES DES EQUATIONS ALGEBRIQUES,
'               TOME 1, MASSON & CIE, 1960, PP.260-266
'----------------------------------------------------------------------
'     Labels: 5, 7, 10, 15, 30, 40, 45, 55, 60, 70, 75
'     int I,IT1,IT2,K,K1,K2,L,M
'     double EPS,EPS1,P,Q,R,TOL,X0,X1,Y0,Y1
'     double CABS,DX,DY,GD,GX,GY,R1,R2,R3,X2,Y2

      IT1 = 25: IT2 = 100: TOL = .01: EPS = 0#
      IER = 0
      IF EPS <> 0! THEN GOTO 7
      EPS = 1#
5     EPS = .5 * EPS
      EPS1 = EPS + 1#
      IF EPS1 > 1# THEN GOTO 5
7     L = N
      X0 = 1#
      Y0 = 1#
      W1(1) = A(1)
      W3(1) = A(1)
      W2(1) = B(1)
      W4(1) = B(1)
10    IF L = 1 THEN GOTO 55
      M = L + 1
      I = N - L + 1
      K1 = 0
      K2 = 0
      X1 = X0
      Y1 = Y0
15    FOR K = 2 TO M
	W1(K) = A(K) + X1 * W1(K - 1) - Y1 * W2(K - 1)
	W2(K) = B(K) + Y1 * W1(K - 1) + X1 * W2(K - 1)
      NEXT K
      FOR K = 2 TO L
	W3(K) = W1(K) + X1 * W3(K - 1) - Y1 * W4(K - 1)
	W4(K) = W2(K) + Y1 * W3(K - 1) + X1 * W4(K - 1)
      NEXT K
      GD = W3(L) * W3(L) + W4(L) * W4(L)
      GX = W1(M) * W3(L) + W2(M) * W4(L)
      GY = W2(M) * W3(L) - W1(M) * W4(L)
      DX = -GX / GD
      DY = -GY / GD
      X2 = X1 + DX
      Y2 = Y1 + DY
      R1 = (ABS(DX) + ABS(DY)) / (ABS(X2) + ABS(Y2))
      K2 = K2 + 1
      X1 = X2
      Y1 = Y2
      IF R1 < TOL THEN GOTO 30
      IF K2 <= IT2 THEN GOTO 15
      IER = I
      GOTO 60
30    K1 = K1 + 1
      IF K1 <= IT1 THEN GOTO 15
      X(I) = X2
      Y(I) = Y2
      IF ABS(X2) <= EPS THEN X2 = EPS
      R2 = ABS(Y2 / X2)
      IF R2 < .1 THEN GOTO 40
      IF R2 * EPS > 1# THEN X(I) = 0#
      X0 = X2
      Y0 = Y2
      GOTO 45
40    IF R2 <= EPS THEN Y(I) = 0#
      X0 = X2 + Y2
      Y0 = X2 - Y2
45   FOR K = 1 TO L
	A(K) = W1(K)
	B(K) = W2(K)
      NEXT K
      L = L - 1
      GOTO 10
55    R3 = W1(1) * W1(1) + W2(1) * W2(1)
      X(N) = -(W1(1) * W1(2) + W2(1) * W2(2)) / R3
      Y(N) = (W2(1) * W1(2) - W1(1) * W2(2)) / R3
      IF ABS(X(N)) <= EPS THEN X(N) = EPS
      R2 = ABS(Y(N) / X(N))
      IF R2 * EPS > 1# THEN X(N) = 0#
      IF R2 <= EPS THEN Y(N) = 0#
60    FOR L = 1 TO N
	P = X(L)
	Q = Y(L)
	R = P * P + Q * Q
	IF L = 1 THEN GOTO 70
	FOR I = L TO 2 STEP -1
	  CABS = X(I - 1) * X(I - 1) + Y(I - 1) * Y(I - 1)
	  IF R >= CABS THEN GOTO 75
	  X(I) = X(I - 1)
	  Y(I) = Y(I - 1)
	NEXT I
70      I = 1
75      X(I) = P
	Y(I) = Q
      NEXT L
    RETURN  'CNEWTON


' end of file tnewton.bas