home *** CD-ROM | disk | FTP | other *** search
/ HAM Radio 1 / HamRadio.cdr / tech / eepub05 / bairstow.pas < prev    next >
Pascal/Delphi Source File  |  1986-01-16  |  8KB  |  214 lines
  1.  
  2. program bairstow;              {  Polynomial root Finder - Bairstow Method  }
  3.  
  4. {    Written by Jeff Crawford                  September 20, 1984
  5.  
  6.                        Transported to TRW 1/17/86
  7.  
  8.      References :
  9.  
  10.      [1]  "Applied Numerical Methods for Digital Computation", M. L. James
  11.           G. M. Smith, J. C. Wolford, Harper & Row, 1977.
  12.      [2]  "Applied Numerical Analysis ", Curtis F. Gerald, Patrick Wheatley
  13.           Addison-Wesley, Third Edition, 1984.                              }
  14.  
  15.  
  16. type    matrix =  array[0..20] of real;
  17.  
  18.  
  19. var          n                    : integer;   {  Order of Polynomial -  is }
  20.                                                {  Gradually Reduced During  }
  21.                                                {  the Procedure             }
  22.              n_save               :  integer;  {  Same as "n" but retains the }
  23.                                                {  Value of the Original Order }
  24.              i                    :  integer;  {  Counter Used in For Loops }
  25.              a                    :  matrix;   {  Coefficient Matrix of Poly  }
  26.              bn                   :  matrix;   {  Coefficients of Reduced Poly}
  27.              cn                   :  matrix;   {  Partial Derivatives Used    }
  28.              Re_root,Im_root      :  matrix;   {  Real & Imaginary Roots      }
  29.              Iteration            :  matrix;   {  Number of Iterations Required}
  30.              dummy                :  real;     {  Variable to Receive Input   }
  31.              u,v                  :  real;     {  Variables Used in Iteration }
  32.                                                {  to Find Roots of Polynomial }
  33.              rad                  :  real;     {  Number Under Radical        }
  34.              count                :  integer;  {  Iteration Counter           }
  35.              Re_remainder         :  matrix;   {  Remainder When Root is Subst}
  36.              Im_remainder         :  matrix;   {  Remainder When Root is Subst}
  37. {                                                                             }
  38. {                                                                             }
  39. label 10,20,30;
  40.  
  41. {  <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> }
  42. {  <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> }
  43.  
  44. procedure bn_cn(var a,bn,cn: matrix; n: integer); { Calculates Bn's & Cn's as }
  45.                    { Described in [1] above, Eqs. 2-69 & 2-74                 }
  46. {                                                                             }
  47. {                                                                             }
  48. var      delta_u   : real; { Fractional Part Added to Each u Each Iteration}
  49.          delta_v   : real; { Fractional Part Added to Each v Each Iteration}
  50.          i         : integer; {Loop Counter }
  51.          Denom     : real; { Denominator Used in Masons Rule to Determine
  52.                              Delta_u and Delta_v }
  53.          stat      : real; { Aborts Loop if > 100 Iteration }
  54. begin
  55. if a[n-2] <> 0.0 then
  56.     begin
  57.     u:= a[n-1]/a[n-2];                {  Initialize Starting Values    }
  58.     v:= a[n]/a[n-2];                  {  Initialize Starting Values    }
  59.     end;
  60. if a[n-2] = 0.0 then   { If the Above Values are Zero, Starting Values of 0.5 }
  61.     begin              {    Are Used Instead                                  }
  62.     u:= 0.5;
  63.     v:= 0.5;
  64.     end;
  65. delta_u:= 1.0;
  66. count := 1;  stat:= 1.0;
  67. while stat > 1.e-9 do
  68.     begin
  69.     bn[0]:= 1.0;
  70.     bn[1]:= a[1] - u;
  71.     for i:= 2 to n do  bn[i]:= a[i]-bn[i-1]*u - bn[i-2]*v; { Eq. 2-69 of [1]}
  72.     cn[0]:= 1.0;
  73.     cn[1]:= bn[1] - u;
  74.     for i:= 2 to n-1 do cn[i]:= bn[i] - cn[i-1]*u - cn[i-2]*v;{Eq. 2-74 of [1]}
  75.     if n <= 3 then Denom:= sqr(cn[n-2])-(cn[n-1]-bn[n-1]){Eqs.(f),(g)pg.165[1]}
  76.        else Denom:= sqr(cn[n-2])-(cn[n-1]-bn[n-1])*cn[n-3];
  77.     delta_u:= (bn[n-1]*cn[n-2] - bn[n]*cn[n-3])/Denom;
  78.     delta_v:= (cn[n-2]*bn[n] - (cn[n-1]-bn[n-1])*bn[n-1])/Denom;
  79.     u:= u + delta_u;                           { Eq. 2-80 of [1] }
  80.     v:= v + delta_v;
  81.     stat:= abs(delta_u) + abs(delta_v);
  82.     if count >= 100 then stat:= 1.e-9;   { Aborts Loop for Over 100 Iterations}
  83.     count:= count + 1;
  84.     end;
  85.     Iteration[n]:= count;
  86.     Iteration[n-1]:= count;
  87.     rad:= sqr(u) - 4*v;
  88.     if rad >= 0.0 then
  89.          begin
  90.          rad:= sqrt(rad);
  91.          Re_root[n-1]:= (-u + rad)/2.0;
  92.          Re_root[n]:= (-u - rad)/2.0;
  93.          Re_remainder[n-1]:= Re_root[n-1]*bn[n-1];
  94.          Re_remainder[n]:= Re_remainder[n-1];
  95.          end;
  96.     if rad < 0.0 then
  97.          begin
  98.          rad:= -1*rad;
  99.          rad:= sqrt(rad);
  100.          Re_root[n-1]:= -u/2.0;
  101.          Re_root[n]:= Re_root[n-1];
  102.          Im_root[n-1]:= rad/2.0;
  103.          Im_root[n]:= -rad/2.0;
  104.          Re_remainder[n]:= (Re_root[n] + u)*bn[n-1] + bn[n];
  105.          Im_remainder[n]:= Im_root[n]*bn[n-1];
  106.          Re_remainder[n-1]:= Re_remainder[n];
  107.          Im_remainder[n-1]:= Im_remainder[n];
  108.          end;
  109. end;
  110.  
  111. { <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> }
  112. { <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> }
  113.  
  114. begin
  115. clrscr;
  116. gotoXY(5,10);
  117. 30:write('Enter the Order of the Polynomial  ');
  118. readln(dummy);
  119. for i:= 0 to 15 do
  120.     begin
  121.     Re_root[i]:= 0.0;
  122.     Im_root[i]:= 0.0;
  123.     Re_remainder[i]:= 0.0;
  124.     Im_remainder[i]:= 0.0;
  125.     end;
  126. n:= Round(dummy);
  127. n_save:= n;
  128. clrscr;
  129. gotoXY(5,10);
  130. writeln('Enter Polynomial Coefficients Beginning with Highest Power of X');
  131. clrscr;
  132. gotoXY(1,1);
  133. for i:= 0 to n do
  134.     begin
  135.     gotoXY(25,2*i+1);
  136.     write('Enter Coefficient  x^',n-i,'  =  ');
  137.     readln(a[i]);
  138.     end;
  139. if a[0] <> 1.0 then
  140.     begin
  141.     Dummy:= a[0];
  142.     for i:= 0 to n do
  143.          begin
  144.          a[i]:= a[i]/Dummy;
  145.          end;
  146.     end;
  147.  
  148. {   Beginning of Root Evaluation  -  First for X to (1) Power, then X to (2)  }
  149. {        and then for X  Raised to 3 or More Power                            }
  150.  
  151. 20: if n = 1 then   {  Case of Single Real Root }
  152.     begin
  153.     Re_root[n]:= -a[1]/a[0];
  154.     Im_root[n]:= 0.0;
  155.     Iteration[n]:= 1;
  156.     n:= n - 1;
  157.     if n = 0 then goto 10;
  158.     end;
  159. if n = 2 then                  { Calculates Root for Order (n) = 2 }
  160.     begin
  161.     Iteration[n]:= 1;
  162.     Iteration[n-1]:= 1;
  163.     rad:= sqr(a[1]) - 4*a[0]*a[2];
  164.     if rad >= 0.0 then                 {  Positive Roots  }
  165.          begin
  166.          rad:= sqrt(rad);
  167.          Re_root[n-1]:= (-a[1]+rad)/(2*a[0]);
  168.          Re_root[n]:= (-a[1] - rad)/(2*a[0]);
  169.          end;
  170.     if rad < 0.0 then      { Quantity Under Radical is (-) so Complex Roots }
  171.          begin
  172.          rad:= -1*rad;
  173.          rad:= sqrt(rad);
  174.          Re_root[n-1]:= -a[1]/(2*a[0]);
  175.          Re_root[n]:= Re_root[n-1];
  176.          Im_root[n-1]:= rad/(2*a[0]);
  177.          Im_root[n]:= -rad/(2*a[0]);
  178.          end;
  179.     n:= n - 2;
  180.     if n = 0 then goto 10;
  181.     end;
  182. if n >= 3 then        { Case for Order > 3  }
  183.     begin
  184.     bn_cn(a,bn,cn,n);
  185.     n:= n - 2;
  186.     if n = 0 then goto 10;
  187.     end;
  188.     for i:= 0 to n do            {  Synthetic Division Performed & Begin With }
  189.          begin                   {  New Polynomial for Root Evaluation        }
  190.          a[i]:= bn[i];
  191.          end;
  192.     goto 20;
  193. 10:clrscr;
  194. writeln;
  195. writeln('                    Roots of Polynomial by Bairstows Method ');
  196. writeln;
  197. i:= n_save;            { Calculation Completed - Output Results }
  198. write('            Real      Imaginary   Iterations       Real      ');
  199. writeln('      Imag ');
  200. write('                                                 Remainder ');
  201. writeln('     Remainder');
  202. writeln;
  203. while i > 0 do
  204.     begin
  205.     write('[',n_save-i+1:3,']    ',Re_root[i]:11:7,'  ',Im_root[i]:11:7);
  206.     write('    ',Iteration[i]:3:0,'         ', Re_Remainder[i]:10,'     ');
  207.     writeln(Im_remainder[i]:10);
  208.     i:= i - 1;
  209.     end;
  210.     writeln;
  211.     goto 30;
  212. end.
  213.  
  214. $