home *** CD-ROM | disk | FTP | other *** search
/ Programmer 7500 / MAX_PROGRAMMERS.iso / PASCAL / ROOTS.ZIP / ROOTS.PAS
Encoding:
Pascal/Delphi Source File  |  1989-06-05  |  2.5 KB  |  117 lines
  1. program roots(input,output);
  2.  
  3. type complex = record
  4.                 r:real;
  5.                 i:real;
  6.                end;
  7.  
  8.      poly = array[0..50] of complex;
  9.  
  10. var P,dP:poly;
  11.     degree,i,i1,i2:integer;
  12.     xpower:array[1..30] of complex;
  13.     f,df,g1,g2,guess:complex;
  14.     conv:boolean;
  15.     ans:char;
  16.     err:real;
  17.  
  18. procedure ad(var outvar:complex; v1,v2:complex);
  19. var work:complex;
  20. begin
  21.   work.r:=v1.r+v2.r;
  22.   work.i:=v1.i+v2.i;
  23.   outvar:=work;
  24. end;
  25.  
  26. procedure mult(var outvar:complex; v1,v2:complex);
  27. var work:complex;
  28. begin
  29.   work.r:=v1.r*v2.r-(v1.i*v2.i);
  30.   work.i:=v1.i*v2.r+(v1.r*v2.i);
  31.   outvar:=work;
  32. end;
  33.  
  34. function mag(v1:complex):real;
  35. begin
  36.   mag:=sqrt(sqr(v1.r)+sqr(v1.i));
  37. end;
  38.  
  39. procedure recip(var v2:complex);
  40. var work:complex;
  41.     d:real;
  42.  
  43. begin
  44.   d:=sqr(v2.r)+sqr(v2.i);
  45.   work.r:=v2.r/d;
  46.   work.i:=-v2.i/d;
  47.   v2:=work;
  48. end;
  49.  
  50. procedure GetComplex(var c:complex);
  51. begin
  52.   writeln;
  53.   write('Re:');
  54.   readln(c.r);
  55.   write('Im:');
  56.   readln(c.i);
  57. end;
  58.  
  59. procedure eval(var res:complex; c:complex;p1:poly;d:integer);
  60. var count:integer;
  61.     t1,power:complex;
  62. begin
  63.   res.r:=0;res.i:=0;
  64.   power.r:=1;power.i:=0;
  65.   for count:=0 to d do begin
  66.     mult(t1,p1[count],power);
  67.     ad(res,res,t1);
  68.     mult(power,power,c);
  69.   end;
  70. end;
  71.  
  72.  
  73. begin   {Main Program}
  74.   err:=0.000001;
  75.   clrscr;
  76.   writeln;
  77.   writeln('  This program finds complex roots of real polynomials.');
  78.   writeln('  By Roger R. Darr - May 27, 1990');
  79.   writeln;
  80.   write('Enter the degree of the polynomial (1-50):');readln(degree);
  81.   writeln;
  82.   for i:=0 to 50 do begin
  83.     p[i].r:=0;p[i].i:=0;
  84.     dp[i].r:=0;dp[i].i:=0;
  85.   end;
  86.   for i:=0 to degree do begin
  87.     write('A',i,'= ');readln(p[i].r);
  88.     if i>0 then dp[i-1].r:=p[i].r*i;
  89.   end;
  90.   repeat
  91.     writeln('Enter an initial guess.  Note: to find complex-conjugate');
  92.     writeln('roots, you must include an imaginary part in your guess.');
  93.     getcomplex(guess);
  94.     g2:=guess;
  95.     i1:=0;
  96.     repeat
  97.       g1:=g2;
  98.       eval(f,g1,p,degree);
  99.       eval(df,g1,dp,degree-1);
  100.       recip(df);
  101.       mult(g2,f,df);
  102.       g2.r:=-g2.r;
  103.       g2.i:=-g2.i;
  104.       ad(g2,g2,g1);
  105.       i1:=i1+1;
  106.       conv:= (abs(g2.i-g1.i)<err) and (abs(g2.r-g1.r)<err);
  107.     until conv or (i1>100);
  108.     if not conv then writeln('Does not converge.')
  109.     else begin
  110.       writeln('Root =  ',g2.r:10:5,'  +   j',g2.i:10:5);
  111.     end;
  112.     writeln;
  113.     write('Again ? ');readln(ans);
  114.   until (ans='n') or (ans='N');
  115. end.
  116.  
  117.