Category : Science and Education
Archive   : EEPUB05.ZIP
Filename : BAIRSTOW.PAS

 
Output of file : BAIRSTOW.PAS contained in archive : EEPUB05.ZIP



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.


$