PROGRAM Chebvyp2;
{Provadi Chebyshevovu interpolaci Rungeho funkce 1/(1 + 25 x^2)
od Xmin do Xmax, pri vypoctu koeficientu se zada rad Nrad
zadana presnost Eps -> Mrad - pocet pouzitych clenu
pocita se funkce, derivave a integral
tiskne se srovnani funkce, Cheb.aprox funkce, derivace, Ch.aprox.der
a Ch.aprox integralu}
const
  np = 60;
TYPE
   Double = Real;
   RealArrayNP = ARRAY [1..np] OF real;
FUNCTION func(y: real): real;
BEGIN
        func :=  1./(1.+25*sqr(y))
END;
PROCEDURE chebft(a,b: real;
               VAR c: RealArrayNP;
                   n: integer);
{Chebyshevuv fit funkce FUNC v intervalu [A,B] stupne maximalne N,
 procedura pocita koeficieny C[1..N] takove, ze
  FUNC ~= (suma od 1 do N z C[K]*T[K-1](y)) - C[1]/2,
      kde y = x-0.5*[b+a]/0.5[b-a]
 Tato procedura ma byt  uzivana s neprilis velkymi N (30 - 50),
 pole C by melo byt zkraceno na mensi hodnote M tak ze C[M+1]
 a dalsi jsou zanedbatelne.}
CONST
   pi = 3.141592653589793;
VAR
   k,j: integer;
   y,fac,bpa,bma: real;
   sum: double;
   f: ^RealArrayNP;
LABEL 1;
BEGIN
   new(f);
   bma := 0.5*(b-a);
   bpa := 0.5*(b+a);
   FOR k := 1 TO n DO BEGIN
      y := cos(pi*(k-0.5)/n);
      f^[k] := func(y*bma+bpa)
   END;
   fac := 2.0/n;
   FOR j := 1 TO n DO BEGIN
      sum := 0.0;
      FOR k := 1 TO n DO
         sum := sum+f^[k]*cos(pi*(j-1)*(k-0.5)/n);
      c[j] := fac*sum
   END;
   dispose(f)
END;
FUNCTION chebev(a,b: real;
              VAR c: RealArrayNP;
                  m: integer;
                  x: real): real;
{Pocita hodnotu funkce v bodu x pomoci Chebyshevova polynomu
 v bode y urceneho z bodu x v intevalu A, B.
 Koeficienty Chebyshevova polynomu byly spocitany pomoci CHEBFT.}
VAR
   d,dd,sv,y,y2: real;
   j: integer;
BEGIN
   IF (x-a)*(x-b) > 0.0 THEN BEGIN
      writeln('pause in CHEBEV - x not in range.');
      readln
   END;
   d := 0.0;
   dd := 0.0;
   y := (2.0*x-a-b)/(b-a);
   y2 := 2.0*y;
   FOR j := m DOWNTO 2 DO BEGIN
      sv := d;
      d := y2*d-dd+c[j];
      dd := sv
   END;
   chebev := y*d-dd+0.5*c[1]
END;
PROCEDURE chder(a,b: real;
         VAR c,cder: RealArrayNP;
                  n: integer);
{Pocita Chebyshevovy konficienty CDER[1,N] derivace funkce
 s pomoci Chebyshevovych polynomu  s koeficienty C[1..N] spoctenymi v CHEBFT. }
VAR
   j: integer;
   con: real;
BEGIN
   cder[n] := 0.0;
   cder[n-1] := 2*(n-1)*c[n];
   IF n >= 3 THEN
      FOR j := n-2 DOWNTO 1 DO
         cder[j] := cder[j+2]+2*j*c[j+1];
   con := 2.0/(b-a);
   FOR j := 1 TO n DO
      cder[j] := cder[j]*con
END;
PROCEDURE chint(a,b: real;
         VAR c,cint: RealArrayNP;
                  n: integer);
{Pocita Chebyshevovy konficienty CINT[1,N] integralu funkce od A do x
 s pomoci Chebyshevovych polynomu  s koeficienty C[1..N] spoctenymi v CHEBFT. }
VAR
   j: integer;
   sum,fac,con: real;
BEGIN
   con := 0.25*(b-a);
   sum := 0.0;
   fac := 1.0;
   FOR j := 2 TO n-1 DO BEGIN
      cint[j] := con*(c[j-1]-c[j+1])/(j-1);
      sum := sum+fac*cint[j];
      fac := -fac
   END;
   cint[n] := con*c[n-1]/(n-1);
   sum := sum+fac*cint[n];
   cint[1] := 2.0*sum
END;

var
   C,cint,cder : RealArrayNP;
   xmax,xmin,delt,eps,sumeps,dx1,yf,ych,y1f,y1ch,yif,yich,xh : Real;
   herr,merr,rerr,mrerr:real;
   nrad,mrad,ndel,i : integer;
   f1 : text;
LABEL 1;
BEGIN
   ASSIGN(f1,'chebv2a.dat');
   REWRITE(f1);
   Writeln('ZADEJ xmin,xmax,nrad');
   Readln(xmin,xmax,nrad);
   Chebft(xmin,xmax,C,nrad);
   Chder(xmin,xmax,C,CDER,nrad);
   Chint(xmin,xmax,C,Cint,nrad);
   mrad := nrad;
   {WHILE true DO BEGIN}
      Writeln('Zadej pozadovanou presnost Cheb. approx');
      Readln(eps);
      IF (eps <= 0) THEN GOTO 1;
      i := nrad+1;
      sumeps := 0;
      REPEAT
        i := i - 1;
        sumeps := sumeps + abs(c[i])
      UNTIL sumeps >= eps;
      mrad := i;
      WRITELN('eps= ',eps,' nrad= ',nrad,' mrad= ',mrad,' sumeps= ',sumeps);
      {END;}
 1:   Writeln(' zadej ndel');
      Readln(ndel);
      dx1 := (xmax-xmin)/(ndel - 1);
      writeln('    x         yt       ych     y1t   y1ch   yich');
      merr:=0;
      mrerr:=0;
   FOR i := 1 to ndel DO BEGIN
       xh := xmin+(i-1)*dx1;
       yf := 1./(1.+25*sqr(xh));
       ych:= chebev(xmin,xmax,c,mrad,xh);
       herr:=abs(ych-yf);
       rerr:=herr/yf;
       if herr>merr then merr:=herr;
       if rerr>mrerr then mrerr:=rerr;
       y1f:= -50*xh/sqr(1.+25*sqr(xh));
       y1ch:= chebev(xmin,xmax,cder,mrad,xh);
       yich:= chebev(xmin,xmax,cint,mrad,xh);
       {writeln(xh,yf,ych,y1f,y1ch,yich);}
       writeln(f1,xh:10:6,' ',yf:13:-4,' ',ych:13:-4,' ',
                  y1f:13:-4,' ',y1ch:13:-4,' ',yich:13:-4)
       END;
   CLOSE(f1);
   Writeln('Max. absolut. chyba  ',merr,'  Max.relat.chyba  ',mrerr);
   Writeln(' V  ',ndel,' bodech');
   Readln;
END.