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