PROCEDURE zbrac(VAR x1,x2: real; VAR succes: boolean);
(* Programs using routine ZBRAC must externally define a
function fx(x:real):real which is to be bracketed. *)
LABEL 99;
CONST
   factor=1.6;
   ntry=50;
VAR
   j: integer;
   f2,f1: real;
BEGIN
   IF (x1 = x2) THEN BEGIN
      writeln('pause in routine ZBRAC');
      writeln('you have to guess an initial range'); readln
   END;
   f1 := fx(x1);
   f2 := fx(x2);
   succes := true;
   FOR j := 1 to ntry DO BEGIN
      IF (f1*f2 < 0.0) THEN GOTO 99;
      IF (abs(f1) < abs(f2)) THEN BEGIN
         x1 := x1+factor*(x1-x2);
         f1 := fx(x1)
      END ELSE BEGIN
         x2 := x2+factor*(x2-x1);
         f2 := fx(x2)
      END
   END;
   succes := false;
99:   END;
PROCEDURE zbrak(x1,x2: real; n: integer;
      VAR xb1,xb2: glnbmax; VAR nb: integer);
(* Programs using routine ZBRAK must externally define a function
fx(x:real):real which is to be analyzed for roots. Also, the
calling program must define type
TYPE
   glnbmax = ARRAY [1..nbmax] OF real;
where nbmax is the maximum number to be required for 'nb'. *)
LABEL 99;
VAR
   nbb,i: integer;
   x,fp,fc,dx: real;
BEGIN
   nbb := nb;
   nb := 0;
   x := x1;
   dx := (x2-x1)/n;
   fp := fx(x);
   FOR i := 1 to n DO BEGIN
      x := x+dx;
      fc := fx(x);
      IF ((fc*fp) < 0.0)  THEN BEGIN
         nb := nb+1;
         xb1[nb] := x-dx;
         xb2[nb] := x
      END;
      fp := fc;
      IF (nbb = nb) THEN GOTO 99;
   END;
99:   END;
FUNCTION zbrent(x1,x2,tol: real): real;
(* Programs using routine ZBRENT must externally define a function
fx(x:real):real a root of which is to be found. *)
LABEL 99;
CONST
   itmax=100;
   eps=3.0e-8;
VAR
   a,b,c,d,e: real;
   min1,min2,min: real;
   fa,fb,fc,p,q,r: real;
   s,tol1,xm: real;
   iter: integer;
BEGIN
   a := x1;
   b := x2;
   fa := fx(a);
   fb := fx(b);
   IF (fb*fa > 0.0) THEN BEGIN
      writeln('pause in routine ZBRENT');
      writeln('root must be bracketed'); readln
   END;
   fc := fb;
   FOR iter := 1 to itmax DO BEGIN
      IF (fb*fc > 0.0) THEN BEGIN
         c := a;
         fc := fa;
         d := b-a;
         e := d
      END;
      IF (abs(fc) < abs(fb)) THEN BEGIN
         a := b;
         b := c;
         c := a;
         fa := fb;
         fb := fc;
         fc := fa
      END;
      tol1 := 2.0*eps*abs(b)+0.5*tol;
      xm := 0.5*(c-b);
      IF ((abs(xm) <= tol1) OR (fb = 0.0)) THEN BEGIN
         zbrent := b; GOTO 99 END;
      IF ((abs(e) >= tol1) AND (abs(fa) > abs(fb))) THEN BEGIN
         s := fb/fa;
         IF (a = c)  THEN BEGIN
            p := 2.0*xm*s;
            q := 1.0-s
         END ELSE BEGIN
            q := fa/fc;
            r := fb/fc;
            p := s*(2.0*xm*q*(q-r)-(b-a)*(r-1.0));
            q := (q-1.0)*(r-1.0)*(s-1.0)
         END;
         IF (p > 0.0) THEN  q := -q;
         p := abs(p);
         min1 := 3.0*xm*q-abs(tol1*q);
         min2 := abs(e*q);
         IF (min1 < min2) THEN min := min1 ELSE min := min2;
         IF (2.0*p < min) THEN BEGIN
            e := d;
            d := p/q
         END ELSE BEGIN
            d := xm;
            e := d
         END
      END ELSE BEGIN
         d := xm;
         e := d
      END;
      a := b;
      fa := fb;
      IF (abs(d) > tol1) THEN BEGIN
         b := b+d
      END ELSE BEGIN
         IF (xm > 0) THEN BEGIN
            b := b+abs(tol1)
         END ELSE BEGIN
            b := b-abs(tol1)
         END
      END;
      fb := fx(b)
   END;
   writeln('pause in routine ZBRENT');
   writeln('maximum number of iterations exceeded'); readln;
   zbrent := b;
99:   END;
PROCEDURE zroots(a: glcarray; m: integer; VAR roots: glcarray;
       polish: boolean);
(* Programs using routine ZROOTS must define the types
TYPE
   glcarray = ARRAY [1..2*m+2] OF real;
   gl2array = ARRAY [1..2] OF real;
in the main routine. *)
LABEL 10;
CONST
   eps=2.0e-6;
VAR
   jj,j,i: integer;
   dum: real;
   b,c,x: gl2array;
   ad: glcarray;
BEGIN
   FOR j := 1 to 2*(m+1) DO BEGIN
      ad[j] := a[j]
   END;
   FOR j := m DOWNTO 1 DO BEGIN
      x[1] := 0.0;
      x[2] := 0.0;
      laguer(ad,j,x,eps,false);
      IF (abs(x[2]) <= (2.0*eps*abs(x[1]))) THEN BEGIN
         x[2] := 0.0
      END;
      roots[2*j-1] := x[1];
      roots[2*j] := x[2];
      b[1] := ad[2*j+1];
      b[2] := ad[2*j+2];
      FOR jj := j DOWNTO 1 DO BEGIN
         c[1] := ad[2*jj-1];
         c[2] := ad[2*jj];
         ad[2*jj-1] := b[1];
         ad[2*jj] := b[2];
         dum := b[1];
         b[1] := b[1]*x[1]-b[2]*x[2]+c[1];
         b[2] := dum*x[2]+b[2]*x[1]+c[2]
      END
   END;
   IF (polish) THEN BEGIN
      FOR j := 1 to m DO BEGIN
         x[1] := roots[2*j-1];
         x[2] := roots[2*j];
         laguer(a,m,x,eps,true);
         roots[2*j-1] := x[1];
         roots[2*j] := x[2]
      END
   END;
   FOR j := 2 to m DO BEGIN
      x[1] := roots[2*j-1];
      x[2] := roots[2*j];
      FOR i := j-1 DOWNTO 1 DO BEGIN
         IF (roots[2*i-1] <= x[1]) THEN GOTO 10;
         roots[2*i+1] := roots[2*i-1];
         roots[2*i+2] := roots[2*i]
      END;
      i := 0;
10:      roots[2*i+1] := x[1];
      roots[2*i+2] := x[2]
   END
END;
