'*****************************************************
'*  Program to demonstrate the BAIRSTOW subroutine   *
'* ------------------------------------------------- *
'* 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-10x^4+35x^3-50x^2+24x       *
'*                                                   *
'* SAMPLE RUN:                                       *
'*                                                   *
'* Input order of polynomial: 5                      *
'*                                                   *
'* Input the polynomial coefficients:                *
'*                                                   *
'*    A(0) = ? 0                                     *
'*    A(1) = ? 24                                    *
'*    A(2) = ? -50                                   *
'*    A(3) = ? 35                                    *
'*    A(4) = ? -10                                   *
'*    A(5) = ? 1                                     *
'*                                                   *
'* Convergence factor: 1e-8                          *
'* Maximum number of iterations: 20                  *
'*                                                   *
'* The roots found are:                              *
'*                                                   *
'*      X1 = 1                                       *
'*      Y1 = 0                                       *
'*                                                   *
'*      X2 = 0                                       *
'*      Y2 = 0                                       *
'*                                                   *
'* The number of iterations was: 15                  *
'*                                                   *
'*****************************************************
defint i-n
defdbl a-h,o-z
cls
dim A(10),B(10),C(10),D(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 BAIRSTOW 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

'**************************************************
'*       Bairstow complex root subroutine         *
'* ---------------------------------------------- *
'* This routine finds the complex conjugate roots *
'* of a polynomial having real coefficients.      *
'* ---------------------------------------------- *
'* Reference: Computer Methods for Science and    *
'*            Engineering by R.L. Lafara.         *
'* ---------------------------------------------- *
'* INPUTS:                                        *
'*  Polynomial coefficients      : A(0) to A(m)   *
'*  Order of polynomial (>=4)    : m              *
'*  Initial guess                : a and b        *
'*  Convergence factor           : e              *
'*  Maximum number of iterations : n              *
'* OUTPUTS:                                       *
'*  Two conjugate complex roots  : x1,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
  'Take initial estimates for aa and bb
  k=0 : B(m)=1#
  'Start iteration sequence
1100 B(m-1)=C(m-1)-aa
  for j=2 to m-1
    B(m-j)=C(m-j)-aa*B(m+1-j)-bb*B(m+2-j)
  next
  B(0)=C(0)-bb*B(2)
  D(m-1)=-1# : D(m-2)=-B(m-1)+aa
  for j=3 to m-1
    D(m-j)=-B(m+1-j)-aa*D(m+1-j)-bb*D(m+2-j)
  next
  D(0)=-bb*D(2)
  d2=-B(2)-bb*D(3)
  dd=D(1)*d2-D(0)*D(2)
  a1=-B(1)*d2+B(0)*D(2) : a1=a1/dd
  b1=-D(1)*B(0)+D(0)*B(1) : b1=b1/dd
  aa=aa+a1 : bb=bb+b1 : k=k+1
  'Test for the number of iterations
  if k>=n then goto 1200
  'Test for convergence
  if ABS(a1)+ABS(b1)>e*e then goto 1100
  'Extract roots from quadratic equation
1200 cc=aa*aa-4#*bb
  'Test to see if a complex root
  if cc>0 then goto 1300
  x1=-aa : x2=x1 : y1=SQR(-cc) : y2=-y1
  goto 1400
1300 x1=-aa+SQR(cc)
  x2=-aa-SQR(cc) : y1=0 : y2=y1
1400 x1=x1/2# : x2=x2/2# : y1=y1/2# : y2=y2/2#
return

'End of file Bairstow.bas