program Lagrange_Prog;
{$I-}        { Disable I/O error trapping }
{$R+}        { Enable range checking }
uses
  Dos, Crt, Common;
const
TNArraySize = 50;            { Size of vectors }
TNNearlyZero = 1E-015;
type
  TNvector = array[0..TNArraySize] of Real;
var
  NumPoints : integer;              { Number of data points }
  XData, YData : TNvector;          { Data points (X,Y) }
  Poly : TNvector;                  { The constructed polynomial }
  Error : byte;                     { Flags if something went wrong }

procedure GetTwoVectorsFromKeyboard(var NumPoints : integer;
                                    var XData     : TNvector;
                                    var YData     : TNvector);
var
  Term : integer;
begin
  NumPoints := 0;
  Writeln;
  repeat
    Write('Number of points (0-', TNArraySize, ')? ');
    Readln(NumPoints);
    IOCheck;
  until ((NumPoints >= 0) and (NumPoints <= TNArraySize) and not IOerr);
  Writeln;
  Write('Enter the X ');
  Writeln('and Y values, separated by a space (not a comma):');
  for Term := 1 to NumPoints do
  repeat
    Write('X[', Term, '], Y[', Term, ']:');
    Read(XData[Term], YData[Term]);
    Writeln;
    { Read in the XData and YData }
    IOCheck;
  until not IOerr;
end; { procedure GetTwoVectorsFromKeyboard }

procedure Lagrange(NumPoints : integer;
               var XData      : TNvector;
               var YData      : TNvector;
               var Poly       : TNvector;
               var Error      : byte);

var
  Degree : integer;        { Degree of resulting polynomial }
                           { one less than NumPoints        }

procedure SynDiv(Degree : integer;
             var Poly   : TNvector;
                 X      : Real;
             var Y      : Real);
var
  Index : integer;

begin
  Y := Poly[Degree];
  for Index := Degree-1 downto 0 do
    Y := Y * X + Poly[Index];
end; { procedure SynDiv }

procedure MakePolynomial(Degree : integer;
                     var XData  : TNvector;
                     var YData  : TNvector;
                     var Poly   : TNvector;
                     var Error  : byte);

var
  T, D : TNvector;       { Iterative variables for computing }
                         { the Lagrange Polynomial           }
  Num, Denom, Coef : Real;
  Term, Term2 : integer;

begin
  for Term := 1 to Degree + 1 do
  begin
    T[Term] := 0;
    D[Term] := 0
  end;
  T[0] := YData[1];
  D[0] := -XData[1];
  D[1] := 1;
  for Term := 1 to Degree do
  begin
    { Evaluate D at XData[Term] }
    SynDiv(Term, D, XData[Term + 1], Denom);
    if ABS(Denom) < TNNearlyZero then
      Error := 1    { Points not unique }
    else
      begin
        { Evaluate T at XData[Term] }
        SynDiv(Term - 1, T, XData[Term + 1], Num);
        Coef := (YData[Term + 1] - Num) / Denom;
        for Term2 := Term downto 0 do
        begin
          T[Term2] := Coef * D[Term2] + T[Term2];
          D[Term2 + 1] := D[Term2] - XData[Term + 1] * D[Term2 + 1];
        end;
        D[0] := -XData[Term + 1] * D[0];
      end;
  end;
  Poly := T;
end; { procedure MakePoly }

begin { procedure Lagrange }
  Degree := NumPoints - 1;
  Error := 0;
  if NumPoints < 1 then
    Error := 2
  else
    begin
      MakePolynomial(Degree, XData, YData, Poly, Error);
    end;
end; { procedure Lagrange }

procedure Results(var XData : TNvector;
                  var YData : TNvector;
                  var Poly  : TNvector;
                  Error     : byte);
var
  Term : integer;

begin
  Writeln('The Data : ');
  for Term := 1 to NumPoints do
    Writeln(XData[Term] : 12 : 7, '             ', YData[Term]);
  if Error >= 1 then
    DisplayError;

  case Error of
    0 : begin
          Writeln(OutFile, 'The polynomial : ');
          for Term := NumPoints - 1 downto 0 do
            Writeln('Poly[', Term : 2, ']=', Poly[Term]);
          Writeln('    X                  Interpolated Y value');
        end;
    1 : Writeln('The data points must be unique.');

    2 : Writeln('There must be at least one data point.');

  end;
  readln;
end; { procedure Results }
begin
GetTwoVectorsFromKeyboard(NumPoints,XData,YData);
Lagrange(NumPoints,XData,YData,Poly,Error);
Results(XData,YData,Poly,Error);
readln;
end.