Формула численного интегрирования интерполяционного типа

var a,zn,h,b,ab,x0,x1,x2,fx0:real;
    f:array[0..100] of real;
    x:array[0..100] of real;
    Pn,Pn2:array[1..4] of real;
    i,j,k,M,l,r,n,right,w:integer;
 
begin
   {  readln(a,b,n);          }
     a:=0;
     b:=1;
     h:=0.1;
     ab:=a;
     x0:=0.35;
     fx0:=3*x0*x0*x0-x0*x0+5*x0-4;
 
     repeat
           x[n]:=ab;
           f[n]:=3*ab*ab*ab-ab*ab+5*ab-4;
           ab:=ab+h;
           n:=n+1;
     until ab>b;
     dec(n);
 
     i:=0;
     repeat
           i:=i+1;
           r:=i;
     until x[i]>x0;
     right:=r;
 
     l:=r;
     for k:=1 to 4 do
     begin
          if k mod 2 > 0 then
            l:=l-1
          else
            r:=r+1;
 
          for j:=l to r do
          begin
               zn:=1;
               i:=l;
               while i<=r do
               begin
                    if i<>j then
                      zn:=zn*(-x[i]+x[j]);
                    inc(i);
               end;
               Pn[k]:=Pn[k]+f[j]/zn;
          end;
     end;
 
     for k:=2 to 4 do
     begin
          ab:=ab*k;
          Pn[k]:=Pn[k]*ab;
     end;
 
     writeln('f^(1)=',Pn[1]:3:5,'    Nev= ',abs(Pn[1]-(9*x0*x0-2*x0+5)):3:5);
     writeln('f^(2)=',Pn[2]:3:5,'    Nev= ',abs(Pn[2]-(18*x0-2)):3:5);
     writeln('f^(3)=',Pn[3]:3:5,'    Nev= ',abs(Pn[3]-18):3:5);
     writeln('f^(4)=',Pn[4]:3:5,'    Nev= ',abs(Pn[4]-0):3:5);
 
     r:=right;
     l:=r-1;
     for k:=1 to 4 do
     begin
          for j:=l to r+w do
          begin
               zn:=1;
               i:=l;
               if j<=r then
                 x2:=x[j]
               else
                 x2:=x0;
 
               while i<=r+w do
               begin
                    if i<=r then
                      x1:=x[i]
                    else
                      x1:=x0;
 
                    if i<>j then
                      zn:=zn*(x2-x1);
                    inc(i);
               end;
 
               if j<=r then
                 Pn2[k]:=Pn2[k]+f[j]/zn
               else
                 Pn2[k]:=Pn2[k]+fx0/zn;
          end;
 
          if k mod 2 > 0 then
            w:=1
          else
          begin
               l:=l-1;
               r:=r+1;
               w:=0;
          end;
     end;
 
     ab:=1;
     for k:=2 to 4 do
     begin
          ab:=ab*k;
          Pn2[k]:=Pn2[k]*ab;
     end;
 
     writeln(#10,'Tochka povishennoy tochnosti:');
     writeln('f^(1)=',Pn2[1]:3:5,'    Nev= ',abs(Pn2[1]-(9*x0*x0-2*x0+5)):3:5);
     writeln('f^(2)=',Pn2[2]:3:5,'    Nev= ',abs(Pn2[2]-(18*x0-2)):3:5);
     writeln('f^(3)=',Pn2[3]:3:5,'    Nev= ',abs(Pn2[3]-18):3:5);
     writeln('f^(4)=',Pn2[4]:3:5,'    Nev= ',abs(Pn2[4]-0):3:5);
 
     readln;
end.