

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.


$ 