(*$B- Aufgabe 4.20 *)
PROGRAM laguerre;
CONST nn=20; eps = 2e-12;
TYPE koeff = ARRAY[0..nn] OF real;
VAR ar,ai,br,bi,c : koeff;  n,i,grad : integer;
    yr,yi,xr,xi,f1,f2,fs1,fs2,fss1,fss2,ffs1,ffs2,fssf1,fssf2,h,yb : real;
    tf: text; stri: string[20];
PROCEDURE csqrt(a,b:real; VAR x,y:real);
(* von Aufgabe 2.21 *)
BEGIN
  IF a>0 THEN y:=b/sqrt(2*(a+sqrt(sqr(a)+sqr(b))))
         ELSE y:= sqrt((-a+sqrt(sqr(a)+sqr(b)))/2);
  IF y=0 THEN
    IF a<0 THEN BEGIN x:=0; y:=sqrt(-a) END
           ELSE x:=sqrt(a)
  ELSE
    BEGIN
       IF (b<0) AND (y>0) THEN y:=-y;
       x:= b/2/y;
    END;
END;
PROCEDURE cmul(a,b,c,d:real; VAR x,y:real);
BEGIN  x:= (a*c-b*d); y:= (b*c+a*d) END;
PROCEDURE cdiv(a,b,c,d:real; VAR x,y:real);
BEGIN
  x:= sqr(c)+sqr(d);
  y:= (b*c-a*d)/x; x:= (a*c+b*d)/x;
END;
PROCEDURE deflation(yr,yi:real);
BEGIN
  br[n-1]:=ar[n]; bi[n-1]:=ai[n];
  FOR i:= n-1 DOWNTO 1 DO
  BEGIN
    br[i-1]:= ar[i]+br[i]*yr-bi[i]*yi;
    bi[i-1]:= ai[i]+br[i]*yi+bi[i]*yr;
  END;
  f1:= ar[0]+br[0]*yr-bi[0]*yi;
  f2:= ai[0]+br[0]*yi+bi[0]*yr;
  writeln(tf,'Nst Nr. ',grad+1-n,'  ',yr,yi,f1,f2);
  br[n]:=yr; bi[n]:=yi; ar:=br; ai:=bi
END;
PROCEDURE laguerreiteration;
VAR fertig: boolean;
BEGIN
   FOR i:=0 TO n DO c[i]:=sqrt(sqr(ar[i]) + sqr(ai[i]));
   REPEAT
     fss1:=0; fss2:=0; fs1:=0; fs2:=0;f1:=0; f2:=0;
     FOR i:=n DOWNTO 0 DO
     BEGIN
       h:=2*fs1+fss1*xr-fss2*xi;
       fss2:=2*fs2+fss2*xr+fss1*xi; fss1:=h;
       h:=f1+fs1*xr-fs2*xi;
       fs2:=f2+fs2*xr+fs1*xi; fs1:=h;
       h:=ar[i]+f1*xr-f2*xi;
       f2 :=ai[i]+f2*xr+f1*xi; f1:=h;
     END;
     writeln(tf,'x =',xr:12:7,xi:12:7,'  P(x) =',f1,f2);
     yb :=sqrt(sqr(xr)+sqr(xi)); h:=0;
     for i:=n downto 0 do h := h*yb+c[i];
     fertig := sqrt(sqr(f1)+sqr(f2))<10*eps*h;
     if not fertig then
     BEGIN
       cdiv(f1,f2,fs1,fs2,ffs1,ffs2);
       cdiv(fss1,fss2,fs1,fs2,fssf1,fssf2);
       cmul(ffs1,ffs2,fssf1,fssf2,yr,yi);
       yr:=sqr(n-1)-n*(n-1)*yr; yi:=-n*(n-1)*yi;
       csqrt(yr,yi,yr,yi);
       cdiv(n*ffs1,n*ffs2,1+yr,yi,yr,yi);
       yr:=xr-yr; yi:=xi-yi; xr:=yr; xi:=yi;
     END
   UNTIL fertig
END;
BEGIN
  writeln('Wohin mit dem Output ?');
  readln(stri); assign(tf,stri); rewrite(tf);
  writeln('Grad a[0],...,a[n]'); read(grad);
  FOR i:=0 TO grad DO read(ar[i],ai[i]);
  yr:=1; yi:=-1;
  FOR n:= grad DOWNTO 1 DO
  BEGIN
    xr:=yr; xi:=-yi;
    laguerreiteration; deflation(xr,xi)
  END;
  writeln(tf,'Nullstellen');
  FOR i:=1 TO grad DO writeln(tf,ar[i],' + I*',ai[i]);
  close(tf);
END.