home *** CD-ROM | disk | FTP | other *** search
/ World of Shareware - Software Farm 2 / wosw_2.zip / wosw_2 / PASCAL / NRPAS13.ZIP / ZROOTS.PAS < prev   
Pascal/Delphi Source File  |  1991-04-29  |  2KB  |  63 lines
  1. PROCEDURE zroots(a: glcarray; m: integer; VAR roots: glcarray;
  2.        polish: boolean);
  3. (* Programs using routine ZROOTS must define the types
  4. TYPE
  5.    glcarray = ARRAY [1..2*m+2] OF real;
  6.    gl2array = ARRAY [1..2] OF real;
  7. in the main routine. *)
  8. LABEL 10;
  9. CONST
  10.    eps=2.0e-6;
  11. VAR
  12.    jj,j,i: integer;
  13.    dum: real;
  14.    b,c,x: gl2array;
  15.    ad: glcarray;
  16. BEGIN
  17.    FOR j := 1 TO 2*(m+1) DO BEGIN
  18.       ad[j] := a[j]
  19.    END;
  20.    FOR j := m DOWNTO 1 DO BEGIN
  21.       x[1] := 0.0;
  22.       x[2] := 0.0;
  23.       laguer(ad,j,x,eps,false);
  24.       IF (abs(x[2]) <= (2.0*eps*abs(x[1]))) THEN BEGIN
  25.          x[2] := 0.0
  26.       END;
  27.       roots[2*j-1] := x[1];
  28.       roots[2*j] := x[2];
  29.       b[1] := ad[2*j+1];
  30.       b[2] := ad[2*j+2];
  31.       FOR jj := j DOWNTO 1 DO BEGIN
  32.          c[1] := ad[2*jj-1];
  33.          c[2] := ad[2*jj];
  34.          ad[2*jj-1] := b[1];
  35.          ad[2*jj] := b[2];
  36.          dum := b[1];
  37.          b[1] := b[1]*x[1]-b[2]*x[2]+c[1];
  38.          b[2] := dum*x[2]+b[2]*x[1]+c[2]
  39.       END
  40.    END;
  41.    IF (polish) THEN BEGIN
  42.       FOR j := 1 TO m DO BEGIN
  43.          x[1] := roots[2*j-1];
  44.          x[2] := roots[2*j];
  45.          laguer(a,m,x,eps,true);
  46.          roots[2*j-1] := x[1];
  47.          roots[2*j] := x[2]
  48.       END
  49.    END;
  50.    FOR j := 2 TO m DO BEGIN
  51.       x[1] := roots[2*j-1];
  52.       x[2] := roots[2*j];
  53.       FOR i := j-1 DOWNTO 1 DO BEGIN
  54.          IF (roots[2*i-1] <= x[1]) THEN GOTO 10;
  55.          roots[2*i+1] := roots[2*i-1];
  56.          roots[2*i+2] := roots[2*i]
  57.       END;
  58.       i := 0;
  59. 10:      roots[2*i+1] := x[1];
  60.       roots[2*i+2] := x[2]
  61.    END
  62. END;
  63.