PROGRAM Dempol;
{Demonstrace vypoctu polynomialni interpolace pomoci Nevillova algoritmu
 Procedure Polint z Numerical Recipes.
 Zada se rovnomerne n bodu jmezi xmin, xmax
 y dano funkci exp(x/alf) a pro ruzna x se demonstruje postup vypoctu
 hodnoty interpolacniho polynomu }
CONST
   np = 30;
TYPE
   RealArrayNP = ARRAY [1..np] OF real;

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]);
   Writeln('  1.krok  ');
   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];
      Writeln(' c[',i:2,'] = ',c^[i]:11:-4,' d[',i:2,'] = ',d^[i]:11:-4);
   END;
   y := ya[ns];
   Writeln(' ns= ',ns:2,' hodnota Y v 1.kroku je ',y);
   ns := ns-1;
   FOR m := 1 TO n-1 DO BEGIN
      Readln;
      Writeln(' Krok ',m+1:2);
      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;
         Writeln(' c[',i:2,'] = ',c^[i]:11:-4,' d[',i:2,'] = ',d^[i]:11:-4);
      END;
      IF 2*ns < n-m THEN BEGIN
         dy := c^[ns+1];
         Writeln(' dy = c[',ns+1:2,'] = ',dy); END
      ELSE BEGIN
         dy := d^[ns];
         Writeln(' dy = d[',ns:2,'] = ',dy);
         ns := ns-1;
      END;
      y := y+dy;
      Writeln(' Hodnota y v ',m+1:2,'-tem kroku je ',y);
   END;
   dispose(d);
   dispose(c)
END;

Function Hod(xh,alf:real):Real;
BEGIN
     Hod:=exp(xh/alf);
END;

VAR
   XHOD,YHOD:RealArrayNP;
   x,y,xmin,xmax,dx,alf,dy:Real;
   i,j,n:Integer;

BEGIN
   Writeln('Zadej stupen polynomu,xmin,xmax,konstantu alf ve fci');
   Readln(n,xmin,xmax,alf);
   writeln('zadane hodnoty');
   dx := (xmax-xmin)/n;
   FOR i:=0 TO n DO BEGIN
      j := i+1;
      IF i < n THEN xhod[j]:=xmin + i*dx
               ELSE xhod[j]:=xmax;
      yhod[j] := Hod(xhod[j],alf);
      writeln(j,xhod[j],yhod[j]);
      END;
   Repeat
      Write(' Zadej x (x <= -1.e+30 konci) >>');
      Readln(x);
      IF x > -1.e+30 THEN BEGIN
          polint(xhod,yhod,n+1,x, y,dy);
          Writeln(' Pro x = ',x:11:-4,' se y = ',y:11:-4,
                  ' pri odhadovane chybe= ',abs(dy):11:-4);
          writeln(' hodnota funkce ',Hod(x,alf):11:-4)
          END;
   Until x <= -1.e+30;
END.