'*****************************************************
'*     Program to demonstrate the Lin's method       *
'* ------------------------------------------------- *
'* Reference: BASIC Scientific Subroutines, Vol. II  *
'* By F.R. Ruckdeschel, BYTE/McGRAWW-HILL, 1981 [1]. *
'* ------------------------------------------------- *
'* Example: Find two complex roots of polynomial:    *
'*            f(x) = x^5+5x^4+10x^3+10x^2+5x+1       *
'*                                                   *
'* SAMPLE RUN:                                       *
'*                                                   *
'* Input order of polynomial: 5                      *
'*                                                   *
'* Input the polynomial coefficients:                *
'*                                                   *
'*    A(0) = ? 1                                     *
'*    A(1) = ? 5                                     *
'*    A(2) = ? 10                                    *
'*    A(3) = ? 10                                    *
'*    A(4) = ? 5                                     *
'*    A(5) = ? 1                                     *
'*                                                   *
'* Convergence factor: 0                             *
'* Maximum number of iterations: 1000                *
'*                                                   *
'* The roots found are:                              *
'*                                                   *
'*      X1 = -.93763285                              *
'*      Y1 = 2.5730332E-02                           *
'*                                                   *
'*      X2 = -.93763285                              *
'*      Y2 = -2.2730332E-02                          *
'*                                                   *
'* The number of iterations was: 1000                *
'*                                                   *
'*****************************************************
defint i-n
defdbl a-h,o-z
cls
dim A(10),B(10),C(10)
print
input " Order of polynomial: ", m
print
print " Input the polynomial coefficients:"
print
for i=0 to m
  print "     A(";i;") = "; : input A(i)
next
print
input " Convergence factor: ", e
input " Maximum number of iterations: ",n
aa=3.1415926535# : bb=SQR(2#)

gosub 1000   'Call LIN subroutine

print
print " The roots found are:"
print
print "      X1 = "; x1
print "      Y1 = "; y1
print
print "      X2 = "; x2
print "      Y2 = "; y2
print
print " The number of iterations was: "; k
print
end

'*************************************************
'*     Polynomial complex roots subroutine       *
'* --------------------------------------------- *
'*     This routine uses the Lin's method        *
'* --------------------------------------------- *
'* Reference: A Practical Guide to Computer      *
'*            Method for Engineers by Shoup      *
'* --------------------------------------------- *
'* INPUTS:                                       *
'*  Polynomial coefficients      : A(0) to A(m)  *
'*  Order of polynomial          : m             *
'*  Initial guess                : a and b       *
'*  Convergence factor           : e             *
'*  Maximum number of iterations : n             *
'* OUTPUTS:                                      *
'*  Two complex roots            : x,y1  x2,y2   *
'*  Number of iterations         : k             *
'*************************************************
1000 'Normalize the A(i) series
  for i=0 to m
    C(i)=A(i)/A(m)
  next
  'Start iteration
  'Set initial guess for the quadratic coefficients
  B(0)=0 : B(1)=0
1100 B(m-1)=C(m-1)-aa
  B(m-2)=C(m-2)-aa*B(m-1)-bb
  for j=3 to m
    B(m-j)=C(m-j)-aa*B(m+1-j)-bb*B(m+2-j)
  next
  'Guard against divide by zero
  if B(2)<>0 then goto 1200
  aa=aa+1e-7 : bb=bb-1e-7
1200 a1=(C(1)-bb*B(3))/B(2)
  b1=C(0)/B(2)
  k=k+1
  'Test for the number of iterations
  if k>=n then goto 1300
  'Test for convergence
  if ABS(aa-a1)+ABS(bb-b1) < e*e then goto 1300
  aa=a1 : bb=b1
  'Next iteration
  goto 1100
1300 aa=a1 : bb=b1
  cc=aa*aa-4#*bb
  'Is there an imaginary part?
  if cc>0 then goto 1400
  y1=SQR(-cc) : y2=-y1 : x1=-aa : x2=x1
  goto 1500
1400 y1=0 : y2=y1 : x1=-aa+SQR(cc) : x2=-aa-SQR(cc)
1500 x1=x1/2#
  x2=x2/2# : y1=y1/2# : y2=y2/2#
return

'End of file Lin.bas