program bairstow; { Polynomial root Finder - Bairstow Method } { Written by Jeff Crawford September 20, 1984 Transported to TRW 1/17/86 References : [1] "Applied Numerical Methods for Digital Computation", M. L. James G. M. Smith, J. C. Wolford, Harper & Row, 1977. [2] "Applied Numerical Analysis ", Curtis F. Gerald, Patrick Wheatley Addison-Wesley, Third Edition, 1984. } type matrix = array[0..20] of real; var n : integer; { Order of Polynomial - is } { Gradually Reduced During } { the Procedure } n_save : integer; { Same as "n" but retains the } { Value of the Original Order } i : integer; { Counter Used in For Loops } a : matrix; { Coefficient Matrix of Poly } bn : matrix; { Coefficients of Reduced Poly} cn : matrix; { Partial Derivatives Used } Re_root,Im_root : matrix; { Real & Imaginary Roots } Iteration : matrix; { Number of Iterations Required} dummy : real; { Variable to Receive Input } u,v : real; { Variables Used in Iteration } { to Find Roots of Polynomial } rad : real; { Number Under Radical } count : integer; { Iteration Counter } Re_remainder : matrix; { Remainder When Root is Subst} Im_remainder : matrix; { Remainder When Root is Subst} { } { } label 10,20,30; { <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> } { <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> } procedure bn_cn(var a,bn,cn: matrix; n: integer); { Calculates Bn's & Cn's as } { Described in [1] above, Eqs. 2-69 & 2-74 } { } { } var delta_u : real; { Fractional Part Added to Each u Each Iteration} delta_v : real; { Fractional Part Added to Each v Each Iteration} i : integer; {Loop Counter } Denom : real; { Denominator Used in Masons Rule to Determine Delta_u and Delta_v } stat : real; { Aborts Loop if > 100 Iteration } begin if a[n-2] <> 0.0 then begin u:= a[n-1]/a[n-2]; { Initialize Starting Values } v:= a[n]/a[n-2]; { Initialize Starting Values } end; if a[n-2] = 0.0 then { If the Above Values are Zero, Starting Values of 0.5 } begin { Are Used Instead } u:= 0.5; v:= 0.5; end; delta_u:= 1.0; count := 1; stat:= 1.0; while stat > 1.e-9 do begin bn[0]:= 1.0; bn[1]:= a[1] - u; for i:= 2 to n do bn[i]:= a[i]-bn[i-1]*u - bn[i-2]*v; { Eq. 2-69 of [1]} cn[0]:= 1.0; cn[1]:= bn[1] - u; for i:= 2 to n-1 do cn[i]:= bn[i] - cn[i-1]*u - cn[i-2]*v;{Eq. 2-74 of [1]} if n <= 3 then Denom:= sqr(cn[n-2])-(cn[n-1]-bn[n-1]){Eqs.(f),(g)pg.165[1]} else Denom:= sqr(cn[n-2])-(cn[n-1]-bn[n-1])*cn[n-3]; delta_u:= (bn[n-1]*cn[n-2] - bn[n]*cn[n-3])/Denom; delta_v:= (cn[n-2]*bn[n] - (cn[n-1]-bn[n-1])*bn[n-1])/Denom; u:= u + delta_u; { Eq. 2-80 of [1] } v:= v + delta_v; stat:= abs(delta_u) + abs(delta_v); if count >= 100 then stat:= 1.e-9; { Aborts Loop for Over 100 Iterations} count:= count + 1; end; Iteration[n]:= count; Iteration[n-1]:= count; rad:= sqr(u) - 4*v; if rad >= 0.0 then begin rad:= sqrt(rad); Re_root[n-1]:= (-u + rad)/2.0; Re_root[n]:= (-u - rad)/2.0; Re_remainder[n-1]:= Re_root[n-1]*bn[n-1]; Re_remainder[n]:= Re_remainder[n-1]; end; if rad < 0.0 then begin rad:= -1*rad; rad:= sqrt(rad); Re_root[n-1]:= -u/2.0; Re_root[n]:= Re_root[n-1]; Im_root[n-1]:= rad/2.0; Im_root[n]:= -rad/2.0; Re_remainder[n]:= (Re_root[n] + u)*bn[n-1] + bn[n]; Im_remainder[n]:= Im_root[n]*bn[n-1]; Re_remainder[n-1]:= Re_remainder[n]; Im_remainder[n-1]:= Im_remainder[n]; end; end; { <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> } { <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> } begin clrscr; gotoXY(5,10); 30:write('Enter the Order of the Polynomial '); readln(dummy); for i:= 0 to 15 do begin Re_root[i]:= 0.0; Im_root[i]:= 0.0; Re_remainder[i]:= 0.0; Im_remainder[i]:= 0.0; end; n:= Round(dummy); n_save:= n; clrscr; gotoXY(5,10); writeln('Enter Polynomial Coefficients Beginning with Highest Power of X'); clrscr; gotoXY(1,1); for i:= 0 to n do begin gotoXY(25,2*i+1); write('Enter Coefficient x^',n-i,' = '); readln(a[i]); end; if a[0] <> 1.0 then begin Dummy:= a[0]; for i:= 0 to n do begin a[i]:= a[i]/Dummy; end; end; { Beginning of Root Evaluation - First for X to (1) Power, then X to (2) } { and then for X Raised to 3 or More Power } 20: if n = 1 then { Case of Single Real Root } begin Re_root[n]:= -a[1]/a[0]; Im_root[n]:= 0.0; Iteration[n]:= 1; n:= n - 1; if n = 0 then goto 10; end; if n = 2 then { Calculates Root for Order (n) = 2 } begin Iteration[n]:= 1; Iteration[n-1]:= 1; rad:= sqr(a[1]) - 4*a[0]*a[2]; if rad >= 0.0 then { Positive Roots } begin rad:= sqrt(rad); Re_root[n-1]:= (-a[1]+rad)/(2*a[0]); Re_root[n]:= (-a[1] - rad)/(2*a[0]); end; if rad < 0.0 then { Quantity Under Radical is (-) so Complex Roots } begin rad:= -1*rad; rad:= sqrt(rad); Re_root[n-1]:= -a[1]/(2*a[0]); Re_root[n]:= Re_root[n-1]; Im_root[n-1]:= rad/(2*a[0]); Im_root[n]:= -rad/(2*a[0]); end; n:= n - 2; if n = 0 then goto 10; end; if n >= 3 then { Case for Order > 3 } begin bn_cn(a,bn,cn,n); n:= n - 2; if n = 0 then goto 10; end; for i:= 0 to n do { Synthetic Division Performed & Begin With } begin { New Polynomial for Root Evaluation } a[i]:= bn[i]; end; goto 20; 10:clrscr; writeln; writeln(' Roots of Polynomial by Bairstows Method '); writeln; i:= n_save; { Calculation Completed - Output Results } write(' Real Imaginary Iterations Real '); writeln(' Imag '); write(' Remainder '); writeln(' Remainder'); writeln; while i > 0 do begin write('[',n_save-i+1:3,'] ',Re_root[i]:11:7,' ',Im_root[i]:11:7); write(' ',Iteration[i]:3:0,' ', Re_Remainder[i]:10,' '); writeln(Im_remainder[i]:10); i:= i - 1; end; writeln; goto 30; end. $