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