{ BEGINENVIRON}
PROGRAM POLYNFIT;

USES Gamma;

CONST
   ndatap = 50;
   map =   20;
TYPE
   RealArrayNDATA = ARRAY [1..ndatap] OF real;
   RealArrayMA = ARRAY [1..map] OF real;
   IntegerArrayMFIT = ARRAY [1..map] OF integer;
   RealArrayMAbyMA = ARRAY [1..map,1..map] OF real;
   RealArrayNPbyMP = ARRAY [1..map,1..1] OF real;
   IntegerArrayNP = ARRAY [1..map] OF integer;

PROCEDURE gaussj(VAR a: RealArrayMAbyMA;
                     n: integer;
                 VAR b: RealArrayNPbyMP;
                     m: integer);
VAR
   big,dum,pivinv: real;
   i,icol,irow,j,k,l,ll: integer;
   indxc,indxr,ipiv: ^IntegerArrayNP;
BEGIN
   new(indxc);
   new(indxr);
   new(ipiv);
   FOR j := 1 TO n DO ipiv^[j] := 0;
   FOR i := 1 TO n DO BEGIN
      big := 0.0;
      FOR j := 1 TO n DO
         IF ipiv^[j] <> 1 THEN
            FOR k := 1 TO n DO
               IF ipiv^[k] = 0 THEN
                  IF abs(a[j,k]) >= big THEN BEGIN
                     big := abs(a[j,k]);
                     irow := j;
                     icol := k
                  END
               ELSE IF ipiv^[k] > 1 THEN BEGIN
                  writeln('pause 1 in GAUSSJ - singular matrix');
                  readln
               END;
      ipiv^[icol] := ipiv^[icol]+1;
      IF irow <> icol THEN BEGIN
         FOR l := 1 TO n DO BEGIN
            dum := a[irow,l];
            a[irow,l] := a[icol,l];
            a[icol,l] := dum
         END;
         FOR l := 1 TO m DO BEGIN
            dum := b[irow,l];
            b[irow,l] := b[icol,l];
            b[icol,l] := dum
         END
      END;
      indxr^[i] := irow;
      indxc^[i] := icol;
      IF a[icol,icol] = 0.0 THEN BEGIN
         writeln('pause 2 in GAUSSJ - singular matrix');
         readln
      END;
      pivinv := 1.0/a[icol,icol];
      a[icol,icol] := 1.0;
      FOR l := 1 TO n DO
         a[icol,l] := a[icol,l]*pivinv;
      FOR l := 1 TO m DO
         b[icol,l] := b[icol,l]*pivinv;
      FOR ll := 1 TO n DO
         IF ll <> icol THEN BEGIN
            dum := a[ll,icol];
            a[ll,icol] := 0.0;
            FOR l := 1 TO n DO
               a[ll,l] := a[ll,l]-a[icol,l]*dum;
            FOR l := 1 TO m DO
               b[ll,l] := b[ll,l]-b[icol,l]*dum
         END
   END;
   FOR l := n DOWNTO 1 DO
      IF indxr^[l] <> indxc^[l] THEN
         FOR k := 1 TO n DO BEGIN
            dum := a[k,indxr^[l]];
            a[k,indxr^[l]] := a[k,indxc^[l]];
            a[k,indxc^[l]] := dum
         END;
   dispose(ipiv);
   dispose(indxr);
   dispose(indxc)
END;
{END GAUSSJ}
PROCEDURE covsrt(VAR covar: RealArrayMAbyMA;
                        ma: integer;
                 VAR lista: IntegerArrayMFIT;
                      mfit: integer);
VAR
   j,i: integer;
   swap: real;
BEGIN
   FOR j := 1 TO ma-1 DO
      FOR i := j+1 TO ma DO covar[i,j] := 0.0;
   FOR i := 1 TO mfit-1 DO BEGIN
      FOR j := i+1 TO mfit DO
         IF lista[j] > lista[i] THEN
            covar[lista[j],lista[i]] := covar[i,j]
         ELSE
            covar[lista[i],lista[j]] := covar[i,j]
   END;
   swap := covar[1,1];
   FOR j := 1 TO ma DO BEGIN
      covar[1,j] := covar[j,j];
      covar[j,j] := 0.0
   END;
   covar[lista[1],lista[1]] := swap;
   FOR j := 2 TO mfit DO
      covar[lista[j],lista[j]] := covar[1,j];
   FOR j := 2 TO ma DO
      FOR i := 1 TO j-1 DO
         covar[i,j] := covar[j,i]
END;
{END COVSRT}
PROCEDURE funcs(x: real;
        VAR afunc: RealArrayMA;
               ma: integer);
VAR
  if1 : INTEGER;
BEGIN
  afunc[1] := 1;
  FOR if1 := 2 TO ma DO
      afunc[if1] := x * afunc[if1 - 1];
END;
{END FUNCS }
PROCEDURE lfit(VAR x,y,sig: RealArrayNDATA;
                     ndata: integer;
                     VAR a: RealArrayMA;
                        ma: integer;
                 VAR lista: IntegerArrayMFIT;
                      mfit: integer;
                 VAR covar: RealArrayMAbyMA;
                 VAR chisq: real);
VAR
   k,kk,j,ihit,i: integer;
   ym,wt,sum,sig2i: real;
   beta: ^RealArrayNPbyMP;
   afunc: ^RealArrayMA;
BEGIN
   new(beta);
   new(afunc);
   kk := mfit+1;
   FOR j := 1 TO ma DO BEGIN
      ihit := 0;
      FOR k := 1 TO mfit DO
         IF lista[k] = j THEN ihit := ihit+1;
      IF ihit = 0 THEN BEGIN
         lista[kk] := j;
         kk := kk+1
      END
      ELSE IF ihit > 1 THEN BEGIN
         writeln('pause in routine LFIT');
         writeln('improper permutation in LISTA');
         readln
      END
   END;
   IF kk <> ma+1 THEN BEGIN
      writeln('pause in routine LFIT');
      writeln('improper permutation in LISTA');
      readln
   END;
   FOR j := 1 TO mfit DO BEGIN
      FOR k := 1 TO mfit DO covar[j,k] := 0.0;
      beta^[j,1] := 0.0
   END;
   FOR i := 1 TO ndata DO BEGIN
      funcs(x[i],afunc^,ma);
      ym := y[i];
      IF mfit < ma THEN
         FOR j := mfit+1 TO ma DO
            ym := ym-a[lista[j]]*afunc^[lista[j]];
      sig2i := 1.0/sqr(sig[i]);
      FOR j := 1 TO mfit DO BEGIN
         wt := afunc^[lista[j]]*sig2i;
        FOR k := 1 TO j DO
            covar[j,k] := covar[j,k]+wt*afunc^[lista[k]];
         beta^[j,1] := beta^[j,1]+ym*wt
      END
   END;
   IF mfit > 1 THEN BEGIN
      FOR j := 2 TO mfit DO
         FOR k := 1 TO j-1 DO covar[k,j] := covar[j,k]
   END;
   gaussj(covar,mfit,beta^,1);
   FOR j := 1 TO mfit DO
      a[lista[j]] := beta^[j,1];
   chisq := 0.0;
   FOR i := 1 TO ndata DO BEGIN
      funcs(x[i],afunc^,ma);
      sum := 0.0;
      FOR j := 1 TO ma DO
         sum := sum+a[j]*afunc^[j];
      chisq := chisq+sqr((y[i]-sum)/sig[i])
   END;
   covsrt(covar,ma,lista,mfit);
   dispose(afunc);
   dispose(beta)
END;
{END LFIT}

{ZACATEK  ********************* MAIN *************************}
CONST
  nbodobr = 100;

VAR
   Vstfile, Vystfile : text;
   tex1,tex2:string[12];
   x, y, sig  : RealArrayNDATA;
   mwt, mfit, i1, i2 , ndata, npev, ma, mam1, ifit: integer;
   chi2,q,xmin,xmax,dx,r : real;
   xbod, ybod : real;
   covar: RealArrayMAbyMA;
   a: RealArrayMA;
   lista: IntegerArrayMFIT;

BEGIN
   writeln('Zadej nazev vstupni file');
   readln(tex1);
   ASSIGN (Vstfile,tex1);
   writeln('Zadej nazev vystupni file');
   readln(tex2);
   ASSIGN (Vystfile,tex2);
   RESET(Vstfile);
   REWRITE(Vystfile);
   READLN(Vstfile,ndata,mwt);
   i1 := 0;
   WHILE (i1 < ndata) AND (NOT EOF(Vstfile)) DO
   BEGIN
      i1 := i1 + 1;
      IF mwt = 0
         THEN BEGIN READLN (Vstfile, x[i1],y[i1]); sig[i1]:=1. END
         ELSE READLN (Vstfile, x[i1],y[i1],sig[i1]);
      IF mwt = 0
         THEN WRITELN (i1,x[i1],y[i1])
         ELSE WRITELN (i1,x[i1],y[i1],sig[i1]);
   END;
   IF i1=0 THEN Writeln('nic jsem neprecetl')
           ELSE IF i1 < ndata
                THEN BEGIN
                     Writeln('nacteno',i1,' udaju < ndat =',ndata);
                     ndata := i1
                     END
                 ELSE Writeln('nacteno',ndata,' udaju');
                 ;
   xmin := x[1];
   xmax := y[1];
   FOR i1 := 2 to ndata DO
       BEGIN
       IF x[i1] < xmin THEN xmin:=x[i1];
       if x[i1] > xmax THEN xmax:=x[i1];
       END;
   WRITELN ('xmin,xmax ',xmin,xmax);

   WRITE('zadej stupen polynomu >> ');
   READLN (mam1);
   ma := mam1 + 1;
   WRITE('zadej pocet pevnych paramatru >> ');
   READLN (npev);
   mfit := ma - npev;
   IF npev > 0 THEN
      BEGIN
         WRITE(' zadej ',npev:2,' mocnin s pevn.koef. >> ');
         FOR i1 := ma DOWNTO mfit+1 DO
                BEGIN
                   READ (lista[i1]);
                   lista[i1] := lista[i1] + 1
                END;
         READLN;
         WRITE(' zadej ',npev:2,' pevnych.koef. >> ');
         FOR i1 := ma DOWNTO mfit+1 DO
                    READ (a[lista[i1]]);
         READLN;
         ifit := 1;
         FOR i1 := 1 TO MA DO BEGIN
           i2 := mfit + 1;
           WHILE (i1 <> lista[i2]) AND (i2 <= ma) DO
                  i2 := i2 + 1;
           IF i2 > ma THEN BEGIN
              lista[ifit] := i1;
              ifit := ifit + 1
              END;
           END;
         ifit := ifit - 1;
         IF ifit > mfit THEN BEGIN
             WRITELN ('zadano malo pevnych mocnin, menime mfit z',mfit:2,
                      'na mfit=',ifit:2 );
             mfit := ifit
             END;
         END
      ELSE
         FOR i1:= 1 TO mfit DO
              lista[i1] := i1;
       ;
   lfit(x,y,sig,ndata,a,ma,lista,mfit,covar,chi2);

   for i1:=1 TO ma DO
             WRITE('a[',i1:2,']=',a[i1]);
   WRITELN;

   for i1:=1 TO ma DO
             WRITE('s[',i1:2,']=',SQRT(covar[i1,i1]));
   WRITELN;

   FOR i1:=1 TO ma-1 DO
      IF covar[i1,i1] <> 0 THEN BEGIN
         FOR i2 := i1+1 TO ma DO
                     IF covar[i2,i2] <> 0 THEN BEGIN
                        r:= covar[i1,i2]/sqrt(covar[i1,i1]*covar[i2,i2]);
                        WRITE('r[',i1:2,',',i2:2,']=',r)
                     END; ;
          WRITELN;
          END; ;
   q := gammq(0.5*(ndata - mfit),0.5*chi2);

   WRITELN('chi2=',chi2:9:3,'   q=',q);

   dx := (xmax - xmin)/(nbodobr - 1);
   xbod := xmin - dx;
   FOR i1 := 1 TO nbodobr DO
       BEGIN
         xbod := xbod + dx;
         ybod := a[ma];
         FOR i2:= ma-1 DOWNTO 1 DO
             ybod:= ybod*xbod + a[i2];
         IF i1 <= ndata
            THEN WRITELN(Vystfile,xbod,' ',ybod,' ',x[i1],' ',y[i1])
            ELSE WRITELN(Vystfile,xbod,' ',ybod)
       END;
   CLOSE(Vystfile);
   Readln
END.