'''This program will display the real and imaginary roots for a polynomial of the form f(x)=A0 + A1*x + A2*x^2 +... ''spagetti code in http://www3.sympatico.ca/ltoms/html/bairstow_s_method.html modified by: ''Adolfo Leon Sepulveda ''Aug-01-2004 ''Version 1.0 ''Test case: ''real and imaginary roots of polynomials ''using bairstow''s method ''form is y = A0 + A1*x + A2*x^2 + A3*x^3 +.... ''degree of polynomial: 4 ''coefficent of A(0)= 12 ''coefficent of A(1)= -19 ''coefficent of A(2)= 12 ''coefficent of A(3)= -6 ''coefficent of A(4)= 1 ''Roots: ''0.5 + 1.6583124 i ''0.5 - 1.6583124 i ''4.0 ''1.0 rem roots of polynomials - bairstow''s method cls print print "real and imaginary roots of polynomials" print "using bairstow''s method" print "form is y = A0 + A1*x + A2*x^2 + A3*x^3 +...." print dim a(22) Const err=.0001 Const cero=1e-19 Const iter=100 InputData a,n print print "Roots:" Bairstow a,n,err,iter End sub InputData( ByRef a(), ByRef n) local i,xa input "degree of polynomial: ";n print for i=0 to n print "coefficent of A(";i;")= "; input " ";xa rem stored in array a in reverse order a(n-i+1)=xa next i end Sub Bairstow(a,n,err,iter) dim b(22) dim e(22) local p, q, p1, q1, r0, r1, v0, v1, s0, s1, d0, d1, d2, s, t local m1, n1, k, i, h, j local const cero = 1e-19 rem branch for special treatment of 1st and 2nd degree equations if n<=2 then OrderLessThan2 a, n Exit Sub endif a(n+2)=0 n1=2*int((n+1)/2) m1 = 1 While true If m1 >= n1 / 2 Then Exit loop p=1 q=1 For k=1 to iter while true rem store coefficients in array b for i=1 to n1+1 b(i)=a(i) next i for j=n1-2 to n1-4 step -2 for i=1 to j+1 b(i+1) = b(i+1) - p*b(i) b(i+2) = b(i+2) - q*b(i) next i next j r0 = b(n1+1) r1 = b(n1) s0 = b(n1-1) s1 = b(n1-2) v0 = -q*s1 v1 = s0 - s1*p d0 = v1*s0 - v0*s1 if abs(d0) >= cero then Exit loop p=p+5 q=q+5 wend d1 = s0*r1 - s1*r0 d2 = r0*v1 - v0*r1 p1=d1/d0 q1=d2/d0 p=p+p1 q=q+q1 if !(Abs(r0)>=err or abs(r1)>=err) then e(m1)=1 Exit For endif if !(abs(p1)>=err or abs(q1)>=err) then e(m1)=2 Exit For endif if p=0 then if !(q=0) Then if !(abs(q1/q)>=err) Then e(m1)=3 Exit For endif endif Else if !(abs(p1/p)>=err) Then if !(q=0) Then if !(abs(q1/q)>=err) Then e(m1)=3 Exit For endif endif endif endif Next k If k > iter Then e(m1) = 4 ''Is Order 2 endif While true s = -p/2 t = s^2 - q if !(t<0) then t=Sqr(t) PrintReal s,t Else t=Sqr(-1*t) PrintImag s,t endif If e(m1) = 4 then Exit Sub for j=1 to n1-1 a(j+1) = a(j+1) - p*a(j) a(j+2) = a(j+2) - q*a(j) next j n1 = n1 - 2 if !(n1>1) Then Exit Sub if n1>3 then Exit loop m1 = m1 + 1 e(m1)=1 p=a(2)/a(1) q=a(3)/a(1) Wend Wend ''m1 end Sub OrderLessThan2(a(),n) If n = 2 Then Order2 a Else print -a(2)/a(1) Endif End sub order2(a()) local s,t a(3) = a(2)*a(2) - 4*a(1)*a(3) s = -a(2) / (2*a(1)) t=Sqr(abs(a(3))) / (2*a(1)) if sgn(a(3))<0 then PrintImag s,t Else PrintReal s,t endif End Sub PrintImag(s,t) Print If Abs(s) > cero print s;" + ";t;" i" print s;" - ";t;" i" Else Print 0;" + ";t;" i" print 0;" - ";t;" i" endif End Sub PrintReal(s,t) print If Abs(s+t) > cero print s+t endif If Abs(s-t) > cero print s-t endif End '