Programmer 7500/MAX_PROGRAMMERS.iso/PASCAL/POLYNOM.ZIP/POLYNOM.PRG

Text File | 1993-04-12 | 5.0 KB | 175 lines

program PolyNom;  {Polynomial Curve Fitter pgm.}
   uses crt;

const
   RowSize   = 250;                                     {Mod. #1}
   ColSize   = 2;
   MaxDegree = 7;                                       {Mod. #2}

type
   Array2Type     = array[1..RowSize, 1..ColSize] of real;
   ArrayType      = array[1..RowSize] of real;
   CoeffArrayType = array[0..30] of real;

var
   DataArray               : Array2Type;
   XArray, YArray          : ArrayType;
   PowerArray, RHS, Coeffs : CoeffArrayType;
   Factors                 : array[1..15, 1..15] of real;
   NumDataPairs, Degree, TwoDegree, NumEqns, Code : integer;
   J, K, L, M, Last                               : integer;
   CharFlag                                       : char;
   TimeToQuit, NoMoreToDo                         : boolean;
   Work, Numer, Denom, Correlation, MeanY, X, Y   : real;

{$I GetNumI.PSL}
{$I GetNumR.PSL}
{$I Load2Arr.PSL}
{$I Power.PSL}
{$I Stats.PSL}

procedure Interpolate;

begin
   NoMoreToDo := false;
   repeat;
      writeln;
      writeln('X value to evaluate or E to End evaluations');
      GetNumR(X, CharFlag, Code);
      case Code of
         -1: if CharFlag = 'E' then
                NoMoreToDo := true;
          0: begin
                Y := 0.0;
                for K := 1 to NumEqns do
                   Y := Y + Coeffs[K] * Power(X, K - 1);
                writeln('Y =', Y)
             end
      end
   until
      NoMoreToDo
end;

procedure SolveEqns;

begin
   TwoDegree := Degree * 2;
   NumEqns   := Degree + 1;
   for J := 0 to TwoDegree do
      begin
         PowerArray[J] := 0.0;
         for K := 1 to NumDataPairs do
            PowerArray[J] := PowerArray[J] + Power(XArray[K], J)
      end;
   for J := 1 to NumEqns do
      begin
         RHS[J] := 0.0;
         for K := 1 to NumDataPairs do
            RHS[J] := RHS[J] + YArray[K] *
                               Power(XArray[K], J - 1)
      end;
   for J := 1 to NumEqns do
      for K := 1 to NumEqns do
         Factors[J,K] := PowerArray[J + K - 2];
   for K := 1 to Degree do
      begin
         M := K + 1;
         L := K;
         repeat
            if abs(Factors[M,K]) > abs(Factors[L,K]) then
               L := M;
            inc(M)
         until
            M > NumEqns;
         for J := K to NumEqns do
            begin
               Work         := Factors[K,J];
               Factors[K,J] := Factors[L,J];
               Factors[L,J] := Work;
            end;
         Work := RHS[K];  RHS[K] := RHS[L];  RHS[L] := Work;
         M := K + 1;
         repeat
            Work         := Factors[M,K] / Factors[K,K];
            Factors[M,K] := 0.0;
            for J := K + 1 to NumEqns do
               Factors[M,J] := Factors[M,J] -
                               Work * Factors[K,J];
            RHS[M] := RHS[M] - Work * RHS[K];
            inc(M)
         until
            M > NumEqns
      end;
   Coeffs[NumEqns] := RHS[NumEqns] / Factors[NumEqns,NumEqns];
   for M := Degree downto 1 do
      begin
         Work := 0.0;
         for J := M + 1 to NumEqns do
            begin
               Work      := Work + Factors[M,J] * Coeffs[J];
               Coeffs[M] := (RHS[M] - Work) / Factors[M,M]
            end
      end;
   writeln;
   writeln('X Power', 'Coefficient':22);
   for J := 1 to NumEqns do
      writeln(J - 1:3, Coeffs[J]:28);
   Numer := 0.0;
   Denom := 0.0;
   for J := 1 to NumDataPairs do
      begin
         Work := 0.0;
         for K := 1 to NumEqns do
            Work := Work + Coeffs[K] * Power(XArray[J], K - 1);
         Numer := Numer + sqr(YArray[J] - Work);
         Denom := Denom + sqr(YArray[J] - MeanY)
      end;
   if Denom = 0.0 then
      Correlation := 1.0
   else
      Correlation := sqrt(1.0 - Numer / Denom);
   writeln;
   writeln('Correlation =', Correlation)
end;

BEGIN
   clrscr;
   writeln('PolyFit - Polynomial curve fitter');
   writeln;
   writeln('The input data must now be read from a disk file.');
   NumDataPairs := 0;
   Load2Arr(DataArray, NumDataPairs, Code);
   writeln;
   if Code <> 0 then
      begin
         writeln('File was not read successfully.');
         exit
      end;
   for J := 1 to NumDataPairs do
      begin
         XArray[J] := DataArray[J,1];
         YArray[J] := DataArray[J,2]
      end;
   Stats(YArray, NumDataPairs, MeanY, Work, Work, Work, Work);
   Last := NumDataPairs - 1;
   if Last > MaxDegree then
      Last := MaxDegree;
   TimeToQuit := false;
   repeat
      writeln;
      writeln('Degree to fit (0-', Last,
                      ') or Q to Quit the program');
      GetNumI(Degree, CharFlag, Code);
      case Code of
         -1: if CharFlag = 'Q' then
                TimeToQuit := true;
          0: if (Degree >= 0) and (Degree <= Last) then
                begin
                   SolveEqns;
                   Interpolate
                end
      end
   until
      TimeToQuit
END.