{Program ukazky Rombergovy integrace pri pouziti obdelnikove metody
(nevlastni meze) - f(x)= x * cos(aksi*x)/(1 + x^6)
od a = 1 do b= \infty (1.e30)
Vysledek - int = vypoctena hodnota int., dss = odhad chyby,
jm = posl. cyklus pridani bodu (jhod = 5, uzito 81 bodu)
vypocet funkce - 3^(jm-1) bodu
}

CONST
   np = 20;
TYPE
   RealArrayNP = ARRAY [1..np] OF real;
VAR
   MidinfIt,jm: integer;
   aksi,eps,dss :real;
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: integer;
   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 = 11;
   jmaxp = 12;
   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);
   h^[1] := 1.0;
   FOR j := 1 TO jmax DO BEGIN
      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;
      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, eps (konec pro eps <= 0 >>');
      READLN(aksi,eps);
      IF eps <= 0. THEN GOTO 1;
      qromo(a,b,sum);
      writeln(' jm= ',jm,' dss= ', dss,' int=  ',sum);
      END;
 1:
END.