Интерполяционный полином Лагранжа (Pascal)

[adada5555]

Добрый день. Задача: по заданным значениям х [7;8;9] и y[7;8;7] вывести коэффициенты полинома лагранжа. На бумаге получается y=-x^2+16x-56, то есть коэффициенты -1, 16, -56.
Есть код программы, есть блок схема. Но не могу разобрать действие с 21 по 39 строку. Объясните, пожалуйста, как происходит расчет в этих строках. Хотя бы как одно значение считается.

Var n,k,m,I,j : integer;
    ct,u:real;
    x,y,r,pt:array[1..100] of real;
begin
repeat
write('Число узлов интерполирования  n=');
readln(n);
until n in [2..100];
writeln('Введите массив узлов интерполирования:');
for i:=1 to n do
 begin
  write('x[',i,']=');
  readln(x[i]);
 end;
writeln('Введите массив значений функции в узлах:');
for i:=1 to n do
 begin
  write('y[',i,']=');
  readln(y[i]);
 end;
for i:=1 to n do R[i]:=0;
for i:=1 to n do
 begin
  ct:=1;pt[1]:=1;k:=1;
  for j:=1 to n do
  if i<>j then
   begin
    ct:=ct*(x[i]-x[j]);
    pt[k+1]:=1;
    if k>1 then
    for m:=k downto 2 do
    pt[m]:=pt[m-1]-pt[m]*x[j];
    pt[1]:=-pt[1]*x[j];
    k:=k+1;
   end;
  u:=y[i]/ct;
  for m:=1 to n do
  r[m]:=r[m]+pt[m]*u;
 end;
writeln('Массив коэффициентов:');
for i:=1 to n do
write('R[', i, ']=', r[i]:6:4,' ');
readln
end.