{Program ukazky Rombergovy integrace pri pouziti obdelnikove metody
(nevlastni meze) - f(x)= x * cos(aksi*x)/(1 + x^6)
Zadavam pocet cyklu  !!
Vysledek - int = vypoctena hodnota int., dss = odhad chyby,
jm = cyklus pridani bodu (od jm=5 se provadi extrapolace na h=0
}

CONST
   np = 50;
TYPE
   RealArrayNP = ARRAY [1..np] OF real;
VAR
   MidinfIt:Longint;
   jm, jmax1: integer;
   aksi,eps,dss :real;
   psum, pdss, zsum: RealArrayNP;
FUNCTION func(x: real): real;
BEGIN
      func:= x/(1 + x*x*x*x*x*x) * cos(aksi*x)
END;
PROCEDURE midinf(aa,bb: real;
                 VAR s: real;
                     n: integer);
{Procedura totozna s MIDPNT az na to, ze body integrace jsou voleny
rovnomerne v 1/X a nikoli v X. To dovoluje, aby hranice BB byla
nejvetsi pripustne kladne realne cislo a hranice AA muze byt zaporne
cislo s nejvetsi absolutni hodnotou. V jednom volani musi mit ale
AA a BB stejne znamenko.}
VAR
   j: Longint;
   x,tnm,sum,del,ddel,b,a: real;

FUNCTION funk(x: real): real;
BEGIN
   funk := func(1.0/x)/sqr(x)
END;

BEGIN
   b := 1.0/aa;
   a := 1.0/bb;
   IF n = 1 THEN BEGIN
      s := (b-a)*funk(0.5*(a+b));
      MidinfIt := 1
   END
   ELSE BEGIN
      tnm := MidinfIt;
      del := (b-a)/(3.0*tnm);
      ddel := del+del;
      x := a+0.5*del;
      sum := 0.0;
      FOR j := 1 TO MidinfIt DO BEGIN
         sum := sum+funk(x);
         x := x+ddel;
         sum := sum+funk(x);
         x := x+del
      END;
      s := (s+(b-a)*sum/tnm)/3.0;
      MidinfIt := 3*MidinfIt
   END
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 qromo(a,b: real;
             VAR ss: real);
{Pocita integral funkce FUNC Rombergovou integraci v otevrenem
intervalu A,B s presnosti EPS, maximalni pocet volani integracni
procedury je JMAX pri zhruba 3 nasobnemu zvetseni poctu pouzitych
bodu pri kazdem volani, K je rad Rombergovy metody. Podle charakteru
integralu je mozno v danem miste programu (oznacenem poznamkou)
zvolit integracni metodu MIDPNT, MIDINF, MIDSQL, MIDSQU}
LABEL 99;
CONST
   jmax = 20;
   jmaxp = 21;
   k = 5;
   km = 4;
TYPE
   RealArrayJMAXP = ARRAY [1..jmaxp] OF real;
VAR
   i,j: integer;
   h,s: ^RealArrayJMAXP;
   c,d: ^RealArrayNP;
BEGIN
   new(h);
   new(s);
   new(c);
   new(d);
   IF jmax1 > jmax THEN BEGIN
          WRITELN (' zadano prilis velke mnozstvi kroku ');
          jmax1 := jmax;
          WRITELN (' zmeneno na maximalni hodnotu = jmax = ',jmax)
          end;
   h^[1] := 1.0;
   FOR j := 1 TO jmax1 DO BEGIN
      writeln (' j v qromo =', j);
      midinf(a,b,s^[j],j);   {ZDE JE MOZNO ZVOLIT INTEGRACNI PROCEDURU}
      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);
        { IF abs(dss) < eps*abs(ss) THEN GOTO 99}
      END;
      IF j < k THEN BEGIN
             psum[j] := s^[j];
             zsum[j] := s^[j];
             IF j = 1 THEN BEGIN pdss[1] := 1.e30 END
                      else BEGIN pdss[j] := abs(s^[j] - psum[j-1]) end
             END
        ELSE BEGIN
             psum[j] := ss;
             zsum[j] := s^[j];
             pdss[j] := abs(dss)
             END;
      s^[j+1] := s^[j];
      h^[j+1] := h^[j]/9.0
   END;
   {writeln('pause in QROMO - too many steps');
   readln;}
99: jm:=j;
   dispose(d);
   dispose(c);
   dispose(s);
   dispose(h)
END;
{zacatek hlavniho programu}
CONST
  a = 1.0;
  b = 1.e30;
LABEL
   1;
VAR
   sum:real;
BEGIN
   WHILE true DO BEGIN
      WRITE (' ZADEJ aksi, jmax1 (konec pro jmax1 <= 0 >>');
      READLN(aksi,jmax1);
      IF jmax1 <= 0. THEN GOTO 1;
      qromo(a,b,sum);
      for jm := 1 TO jmax1 DO
      writeln('jm= ',jm,'dss= ', pdss[jm],'int=  ',psum[jm],
              ' intbp= ',zsum[jm]);
      END;
 1:
END.