'*********************************************
'*        Solving algebraic equations         *
'*           by Bairstow's Method             *
'* ------------------------------------------ *
'* Ref.: "Mathematiques par l'informatique    *
'*        individuelle - Programmes en basic  *
'*        by H. Lehning and D. Jakubowicz,    *
'*        MASSON Paris, 1982" [BIBLI 06].     *
'*                                            *
'* SAMPLE RUN:                                *
'* (Find real and complex roots of equation:  *
'*  x^6 -127x^5 +215x^4 +28x^3 -39x^2 +20x    *
'*  -15 = 0)                                  *
'*                                            *
'* Input degree of equation: 6                *
'*   A(6) = 1                                 *
'*   A(5) = -127                              *
'*   A(4) = 215                               *
'*   A(3) = 28                                *
'*   A(2) = -39                               *
'*   A(1) = 20                                *
'*   A(0) = -15                               *
'*                                            *
'*             ROOTS                ERROR     *
'*   0.03989619     0.44667179  1.874778E-011 *
'*   0.03989619    -0.44667179  1.874778E-011 *
'*   0.52383379     0.00000000  1.297698E-006 *
'*  -0.64575255     0.00000000  3.497111E-006 *
'* 125.28210889     0.00000000  3.889075E-009 *
'*   1.76001412     0.00000000  1.438088E-006 *
'*                                            *
'* NOTE: if a better precision is desired,    *
'*       you can set e to a smaller value or  *
'*       use Newton's method for a real root. *
'**********************************************
'See explanation file Bairstow.txt
'---------------------------------
defdbl a-h,o-z
defint i-n
dim A(20),B(20),C(20),P(20)
F$="####.########"
'Enter polynomial A(i) and copy in P(i)
cls
print
input " Input degree of equation: ", N
M=N
for I=0 to N
  print "   A(";N-i;") = "; : input "", A(i)
  P(I)=A(I)
next I
print
print "             ROOTS                 ERROR"
'Init
p1=0 : q=0 : k=100 : e=0.001
'Factoring loop
100 if N<=2 then goto 500
J=0
200 if J>K then goto 1000
J=J+1
'calculate B(i) and C(i)
B(0)=0 : B(1)=0 : C(0)=0 : C(1)=0
for I=2 to N+2
  B(I)=A(I-2)-p1*B(I-1)-q*B(I-2)
  C(I)=-B(I-1)-p1*C(I-1)-q*C(I-2)
next I
'calculate dp=a1 and dq=b1
x=B(N+1) : y=B(N+2) : z=C(N)
t=C(N+1) : u=C(N+2)
d=t*t-z*(u+x)
if d=0 then goto 1000
a1=(z*y-x*t)/d
b1=(-x*(q*z+p1*t)-y*t)/d
'New p1 and q
p1=p1+a1 : q=q+b1
f=(ABS(a1)+ABS(b1))/(ABS(p1)+ABS(q))
if F>E then goto 200
'A factor has been found
gosub 700
'Change polynomial
N=N-2
for I=0 to N
  A(I)=B(I+2)
next I
goto 100
500 'Last factor
if N=2 then goto 530
x=-A(1)/A(0) : Y=0
gosub 800 'print results
goto 550
530 p1=A(1)/A(0) : q=A(2)/A(0)
gosub 700
550 END

700 'Solve x^2+p1x+q=0
  d=p1*p1-4#*q
  if d<0 then goto 750
  d=SQR(d) : y=0
  x=(-p1+d)/2# : gosub 800
  x=(-p1-d)/2# : gosub 800
return
750 d=SQR(-d)/2# : x=-p1/2#
  y=d : gosub 800
  y=-d : gosub 800
return

800 'print rsults
  print using F$; x; : print "  ";
  print using F$; y;
'calculate error estimation
  u=P(0) : v=0
  for I=1 to M
    r=u*x-v*y : v=u*y+v*x : u=r+P(I)
  next I
  a1=SQR(u*u+v*v)
  u=0 : v=0
  for I=1 to M
    r=u*x-v*y : v=u*y+v*x : u=r+(M-I+1)*P(I-1)
  next I
  if A1=0 then goto 850
  a1=a1/SQR(u*u+v*v)
850 print "  "; a1
return

1000 print " Process not convergent !"
END

'end of file bairsto1.bas