'***************************************************** '* 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 