{$A+,B-,D+,E+,F-,G-,I+,L+,N-,O-,P-,Q-,R-,S+,T-,V+,X+}
{$M 16384,0,655360}
{Program ukazky Rombergovy integrace pri pouziti lichobeznikove
metody (vlastni meze s dobrym vypoctem funkce na okraji) - f(x)=cos(x)
Vysledek - int = vypoctena hodnota int., dss = odhad chyby,
jhod = posl. cyklus pridani bodu (jhod = 5, uzito 17 bodu)
aroz = chyba vypoctu bez polynom extrapolace (dela se poprve pro j=5)
teoret. hodnota intregalu}
CONST
   np = 100;
TYPE
   RealArrayNP = ARRAY [1..np] OF real;
VAR
   TrapzdIt: Integer;
   as : RealArrayNP;
   eps,xmin,xmax,sum,adss,aroz:real;
   jhod :integer;
FUNCTION func(x: real): real;
BEGIN
   func := cos(x);
END;
PROCEDURE polint(VAR xa,ya: RealArrayNP;
                         n: integer;
                         x: real;
                  VAR y,dy: real);
{
 Pocita hodnotu y interpolacniho polynomu  n-1  stupne v bodu x,
 dy je odhad chyby y
 zadane hodnoty funkce jsou xa[n],ya[n]
}
VAR
   ns,m,i: integer;
   w,hp,ho,dift,dif,den: real;
   c,d: ^RealArrayNP;
BEGIN
   new(c);
   new(d);
   ns := 1;
   dif := abs(x-xa[1]);
   FOR i := 1 TO n DO BEGIN
      dift := abs(x-xa[i]);
      IF dift < dif THEN BEGIN
         ns := i;
         dif := dift
      END;
      c^[i] := ya[i];
      d^[i] := ya[i]
   END;
   y := ya[ns];
   ns := ns-1;
   FOR m := 1 TO n-1 DO BEGIN
      FOR i := 1 TO n-m DO BEGIN
         ho := xa[i]-x;
         hp := xa[i+m]-x;
         w := c^[i+1]-d^[i];
         den := ho-hp;
         IF den = 0.0 THEN BEGIN
            writeln ('pause in routine POLINT');
            readln
         END;
         den := w/den;
         d^[i] := hp*den;
         c^[i] := ho*den
      END;
      IF 2*ns < n-m THEN
         dy := c^[ns+1]
      ELSE BEGIN
         dy := d^[ns];
         ns := ns-1
      END;
      y := y+dy
   END;
   dispose(d);
   dispose(c)
END;

PROCEDURE trapzd(a,b: real;
               VAR s: real;
                   n: integer);
{Tato procedura pocita integral funkce FUNC od A do B, pro N=1 je
intgral dan S = 0.5 * (B-A) *(FUNC(A) + FUNC(B)), pri nasledujicich
volanich s N = 2, 3, .. dochazi k postupnemu zpresnovani pridanim
2**(N-2) vnitrnich bodu. Je vnitrni procedurou pro integrovani pomoci
QTRAP, QSIMP, QROMB }
VAR
   j: integer;
   x,tnm,sum,del: real;
BEGIN
   IF n = 1 THEN BEGIN
      s := 0.5*(b-a)*(func(a)+func(b));
      TrapzdIt := 1
   END
   ELSE BEGIN
      tnm := TrapzdIt;
      del := (b-a)/tnm;
      x := a+0.5*del;
      sum := 0.0;
      FOR j := 1 TO TrapzdIt DO BEGIN
         sum := sum+func(x);
         x := x+del
      END;
      s := 0.5*(s+(b-a)*sum/tnm);
      TrapzdIt := 2*TrapzdIt
   END
END;
PROCEDURE qromb(a,b: real;
             VAR ss: real);
{Pocita integral S funkce FUNC od A do B pomoci Rombrgovy metody radu 2*K
s presnosti zadanou parametrem EPS, maximalni pocet volani TRAPZD je zadan
parametrem JMAX, tj. maximalni pocet bodu je 2**(JMAX-1)+1 }
LABEL 99;
CONST
   jmax = 16;
   jmaxp = 17;
   k = 5;
TYPE
   RealArrayJMAXP = ARRAY [1..jmaxp] OF real;
VAR
   i,j: integer;
   dss: real;
   h,s: ^RealArrayJMAXP;
   c,d: ^RealArrayNP;
BEGIN
   new(h);
   new(s);
   new(c);
   new(d);
   h^[1] := 1.0;
   FOR j := 1 TO jmax DO BEGIN
      trapzd(a,b,s^[j],j);
      as[j] := s^[j];
      IF j >= k THEN BEGIN
         FOR i := 1 TO k DO BEGIN
            c^[i] := h^[j-k+i];
            d^[i] := s^[j-k+i]
         END;
         polint(c^,d^,k,0.0,ss,dss);
         adss := dss;
         jhod := j;
         IF (j>k) and (abs(dss) < eps*abs(ss)) THEN GOTO 99
      END;
      s^[j+1] := s^[j];
      h^[j+1] := 0.25*h^[j]
   END;
   writeln('pause in QROMB - too many steps');
   readln;
99:
   dispose(d);
   dispose(c);
   dispose(s);
   dispose(h)
END;
  Var teorint:real;
BEGIN
   writeln('zadej xmin,xmax,eps');
   readln(xmin,xmax,eps);
   qromb(xmin,xmax,sum);
   aroz := as[jhod] - as[jhod-1];
   writeln('int =',sum,' dss= ',adss,' jhod= ',jhod,' aroz= ',aroz);
   teorint := sin(xmax) - sin(xmin);
   writeln('teoreticka hodnota integralu = ',teorint);
   readln;
END.