{Hled nˇ ©e¨enˇ rovnice atan((x-C)/Del) + D = 0 Brentovou metodou,
porovn v  se po‡et krok– s p–lenˇm interval–
 }
VAR
   Chod,Dhod,Delth:Real;

FUNCTION fx(x: real): real;
BEGIN
     fx := ArcTan((x-Chod)/Delth) + Dhod;
END;

PROCEDURE zbrac(VAR x1,x2: real;
               VAR succes: boolean);
{Hleda oblast korene funkce FX. Vychazi z pocatecniho odhadu
 okraju X1, X2. Postupne rozsiruje oblast dokud nema funkce v obou
 krajnich bodech opacna znamenka nebo dokud neni interval neprijatelne
 velky - pak je SUCCES = false.}
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;

FUNCTION rtbis(x1,x2,xacc: real): real;
{Nalezne koren funkce FX ohraniceny intervalem X1, X2 metodou puleni intervalu.
 Koren urceny funkci RTBIS ma zadanou presnost +- xacc }
LABEL 99;
CONST
   jmax = 40;
VAR
   dx,f,fmid,xmid,rtb: real;
   j: integer;
BEGIN
   fmid := fx(x2);
   f := fx(x1);
   IF f*fmid >= 0.0 THEN BEGIN
      writeln('pause in RTBIS');
      writeln('Root must be bracketed for bisection.');
      readln
   END;
   IF f < 0.0 THEN BEGIN
      rtb := x1;
      dx := x2-x1
   END
   ELSE BEGIN
      rtb := x2;
      dx := x1-x2
   END;
   FOR j := 1 TO jmax DO BEGIN
      dx := dx*0.5;
      xmid := rtb+dx;
      fmid := fx(xmid);
      IF fmid <= 0.0 THEN rtb := xmid;
      IF (abs(dx) < xacc) OR (fmid = 0.0) THEN GOTO 99
   END;
   writeln('pause in RTBIS - too many bisections');
   readln;
99: WRITELN(' pocet iteraci v bisekci= ',j);
   rtbis := rtb

END;

FUNCTION zbrent(x1,x2,tol: real): real;
{Nalezne koren funkce FX ohraniceny intervalem X1, X2 Brentovou metodou .
 Koren urceny funkci ZBRENT ma zadanou absolutni presnost tol, relativni
 presnost je urcena konstantou eps - strojova presnost}
LABEL 99;
CONST
   itmax = 100;   {Maximalni pocet iteraci}
   {eps = 3.0e-8;  presnost cisel}
   eps = 3.0e-12;  {zmensena az prilis kvuli demo - vhodne 1.e-10}
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;                 {prejmenovani a,b,c, a stanoveni}
         fc := fa;               {delky intervalu d = b -a}
         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;     {pri stanoveni presnosti}
      xm := 0.5*(c-b);                    {vyuzivam zadanou i strojovou}
      IF (abs(xm) <= tol1) OR (fb = 0.0) THEN BEGIN
         zbrent := b;                      {dosazeni presnosti}
         GOTO 99
      END;
      IF (abs(e) >= tol1) AND (abs(fa) > abs(fb)) THEN BEGIN
         s := fb/fa;                  {pokus o inverzni kvadratickou interp}
         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;              {prijeti kvadr.interp.}
            d := p/q;
            writeln(iter,'. it - provedena kvadraticka interpolace');
         END
         ELSE BEGIN
            d := xm;
            e := d;
            writeln(iter,'. it - provedeno puleni intervalu');
         END;
      END
      ELSE BEGIN
         d := xm;
         e := d;
            writeln(iter,'. it - provedeno puleni intervalu');
      END;
      a := b;
      fa := fb;
      IF abs(d) > tol1 THEN b := b+d
      ELSE BEGIN
         IF xm >= 0 THEN b := b+abs(tol1)
         ELSE b := b-abs(tol1)
      END;
      fb := fx(b);
      writeln('xm,d ',xm,'  ',d,' xiter ',b,'  fiter ',fb);
      writeln('  a  b  c ',a,b,c);
      writeln(' fa fb fc ',fx(a),fx(b),fx(c));
      IF ITER MOD 3 = 0 THEN BEGIN
            Writeln('Pro pokracovani stiskni ENTER');Readln;END
   END;
   writeln('pause in routine ZBRENT');
   writeln('maximum number of iterations exceeded');
   readln;
   zbrent := b;
99: Writeln('Pocet iteraci v procedure Zbrent= ',iter)
END;

VAR
   i:integer;
   tol,x1,x2,xk,f1,f2,k,xz1,xz2:Real;
   nalez:Boolean;
LABEL 1,2;
BEGIN
    WHILE true DO BEGIN
      Writeln('Zadej parametry funkce Chod,Delth,Dhod (|Chod|>0.9E30 konci)');
      Readln(Chod,Delth,Dhod);
      IF abs(Chod) > 0.9e30  THEN GOTO 1;
      WHILE true DO BEGIN
             Writeln('Zadej meze x1,x2,presnost (|x1|>0.9E30 konci)');
             Readln(x1,x2,tol);
             IF abs(x1) > 0.9e30  THEN GOTO 2;
             f1:=fx(x1); f2:=fx(x2);
             Writeln('f1=',f1:12,'   f2',f2:12);
             nalez:=false;
             Repeat
               zbrac(x1,x2,nalez)
             Until nalez;
             Writeln('x1=',x1:12,'   x2',x2:12);
             f1:=fx(x1); f2:=fx(x2);
             Writeln('f1=',f1:12,'   f2',f2:12);
             xz1:=x1;xz2:=x2;
             xk:=zbrent(x1,x2,tol);
             Writeln('Vysledek');
             f1:=fx(xk);
             Writeln('xk=',xk:12,'   f1',f1:12);
             xk:=rtbis(xz1,xz2,tol);
             Writeln('Vysledek bisekce');
             f1:=fx(xk);
             Writeln('xk=',xk:12,'   f1',f1:12);
             Readln;
      END;
2:    Writeln('Prechzim na nove zadani konstant funkce');
      END;
1:    Writeln('Koncim...     Stiskni ENTER');Readln
END.