DELPHI - How to make Lagrange interpolation

{ ---------------------------------------------------------------------------
  Calc Lagrange polynomial.
  -------------------------------------------------------------------------- }
function TForm1.CalcLagrange( _iX : double ) : double;
var
  i, j : integer;
  iNumerator, iDenominator : double;
begin
  result := 0;

  for i := 1 to pX.Count - 1 do
    begin
      iNumerator := 1;
      iDenominator := 1;

      for j := 1 to pX.Count - 1 do
        if j <> i then
          begin
            iNumerator := iNumerator * ( _iX - pX.GetValue( j ) );
            iDenominator := iDenominator * ( pX.GetValue( i ) - pX.GetValue( j ) );
          end;

      result := result + ( pY.GetValue( i ) * ( iNumerator / iDenominator ) );    
    end;
end;

{ ---------------------------------------------------------------------------
  Lagrange interpolation.
  -------------------------------------------------------------------------- }
procedure TForm1.Lagrange;
var
  j, i, iIndex : integer;
  iNumSegmentsCount : integer;
  iT, i1_T, iB0, iB1, iB2, iB3 : double;
  iNewX, iNewY : double;
  pSerieLine : TLineSeries;
  iValue : extended;
begin

  pSerieLine := TLineSeries( Chart1.Series[1] );
  pSerieLine.Clear;  

  { draw count of cubic curves }
  for iIndex := 0 to px.Count - 2 do
    begin

       { every curve generate from I parts }
       for i := 0 to 10 - 1 do
         begin
           iT := i / (10-1);

           iNewX := pX.GetValue( iIndex ) + 
                    ( (pX.GetValue( iIndex + 1 ) - pX.GetValue( iIndex )) * iT );

           iNewY := CalcLagrange( iNewX );          

           pSerieLine.AddXY( iNewX, iNewY );
         end;

  end;
end;