Построить интерполяционный многочлен Ньютона

program NewtonPolynomialForward;
 
const
  K = 6;
type
  TList = array [0..K - 1] of real;
const
  TableX: TList = (-2.5, -2.0, -1.0, 0.5, 1.0, 2.0);
 
  function sh(x: real): real;
  var
    t: real;
  begin
    t  := exp(x);
    sh := (1 / t + t) / 2;
  end;
 
  procedure FillTable(var AY: TList);
  var
    i: integer;
  begin
    for i := 0 to K - 1 do
      AY[i] := sh(TableX[i]);
  end;
 
  procedure NewtonCoefficientCalc(Amount: integer; const X, Y: TList; var Coeff: TList);
  var
    i, j: integer;
  begin
    Coeff := Y;
    for i := 1 to Amount - 1 do
      for j := Amount - 1 downto i do
        Coeff[j] := (Coeff[j] - Coeff[j - 1]) / (X[j] - X[j - i]);
  end;
 
  function NewtonPolynomial(x: real; TableLen: integer;
  const TableX, Coefficients: TList): real;
  var
    F: real;
    i: integer;
  begin
    F := Coefficients[TableLen - 1];
    for i := TableLen - 2 downto 0 do
      F := F * (X - TableX[i]) + Coefficients[i];
    NewtonPolynomial := F;
  end;
 
  procedure Show(XStart, XFinish, DeltaX: real; TableLen: integer;
  const TableX, Coefficients: TList);
  var
    x, YInter, YSharp: real;
  begin
    writeln('x': 6, 'sh(x)': 8, 'P(x)': 8, '| sh(x) - P(x) |': 20);
    x := XStart;
    repeat
      YSharp := sh(x);
      YInter := NewtonPolynomial(x, K, TableX, Coefficients);
      writeln(x: 8: 1, YSharp: 8: 4, YInter: 8: 4, abs(YSharp - YInter): 8: 4);
      x := x + DeltaX;
    until x > XFinish;
  end;
 
var
  TableY: TList;
  Coefficients: TList;
begin
  FillTable(TableY);
  NewtonCoefficientCalc(K, TableX, TableY, Coefficients);
  Show(-3, 3, 0.1, K, TableX, Coefficients);
end.