Радиус и центр окружности.

uses crt;
const nmax=28;
      e=0.01;
type
Point=record {точка}
      x,y:real;{координаты}
      end;
Okruz=record {окружность}
      o:Point;{координаты центра}
      r:real;{радиус}
    end;
MassP=array[1..nmax] of Point; {множество-массив точек}
procedure Tochki(var t:MassP;var n:integer);{создание множества точек}
var i:integer;
begin
repeat
write('Количество точек 3 do ',nmax,' n=');
readln(n);
until n in [3..nmax];
for i:=1 to n do
 begin
  t[i].x:=10*random;
  t[i].y:=10*random;
 end;
write('X:');
for i:=1 to n do
write(t[i].x:5:2);
writeln;
write('Y:');
for i:=1 to n do
write(t[i].y:5:2);
writeln;
end;
{вычисление параметров окружности, проходящей через три точки}
procedure Okr(t1,t2,t3:Point;var ok:Okruz);
var a,b,x,y:real;
    k0,k1,k2,m0,m1,m2:real;
begin
k0:=sqr(t1.x)-sqr(t2.x)+sqr(t1.y)-sqr(t2.y);
k1:=2*(t1.y-t2.y);
k2:=2*(t1.x-t2.x);
m0:=sqr(t1.x)-sqr(t3.x)+sqr(t1.y)-sqr(t3.y);
m1:=2*(t1.y-t3.y);
m2:=2*(t1.x-t3.x);
a:=k2*m0-k0*m2;
b:=k2*m1-k1*m2;
if b=0 then exit;
y:=a/b;
ok.o.y:=y;
if abs(m2) > e then x:=(m0-y*m1)/m2
else
  if abs(k2) > e then x:=(k0-y*k1)/k2
  else exit;
ok.o.x:=x;
ok.r:=sqrt(sqr(t1.x-x)+sqr(t1.y-y));
end;
{определение принадлежности точки окружности}
function Prin(a:Point;ok:Okruz):boolean;
begin
if (sqr(a.x-ok.o.x)+sqr(a.y-ok.o.y))-sqr(ok.r)<-e then Prin:=true
else Prin:=false
end;
 
var t:MassP;
    ok:Okruz;
    n,i,j,k,l,kol,mx:integer;
    mxr,mxx,mxy:real;
begin
clrscr;
randomize;
Tochki(t,n);
mx:=0;
mxr:=0;
mxx:=0;
mxy:=0;
for i:=1 to n-2 do
for j:=i+1 to n-1 do
for k:=j+1 to n do
 begin
  Okr(t[i],t[j],t[k],ok);
  kol:=0;
  for l:=1 to n do
  if Prin(t[l],ok)then kol:=kol+1;
  if kol>mx then
   begin
    mx:=kol;
    mxr:=ok.r;
    mxx:=ok.o.x;
    mxy:=ok.o.y;
   end;
 end;
if mx=0 then write('Окружности, содержащей точки внутри нет!')
else
 begin
  writeln('Окружность: координаты центра x0=',mxx:0:2,' Y0=',mxy:0:2,' R=',mxr:0:2);
  writeln('Внутри точек=',mx);
 end;
readln;
end.

program Triangle;
 
TYPE
  Points=array[1..2] of real;
  mat3=array[1..3,1..3] of real;
  mat4=array[1..4,1..4] of real;
var
    p:array[1..30] of points;
    i,j,k,l,{счетчики}
    n,{количество точек}
    vnutri,vne, {количество точек внутри и вне окружности}
    vnutriTemp, vneTemp {количество точек внутри и вне окружности, временные переменные}
   :integer;
    matr3:mat3;{квадратная матрица 3-го порядка}
    matr4:mat4; {квадратная матрица 4-го порядка}
    x0,y0,{центр окружности}
    a,b,c,{стороны треугольника}
    r,square,{радиус и площадь}
    det,{определитель}
    pol,{полупериметр}
    pk{подкоренное ввыражение, формуле герона}
    :real;
 
 
    {функция нахождения длины между точками}
    function dlina(p1,p2:points):real;
      begin
        dlina:=sqrt(sqr(p1[1]-p2[1])+sqr(p1[2]-p2[2]));
      end;
 
    {функция нахождения определителя 3-го порядка}
    function determine3(A:mat3):real;
      begin
      determine3:=a[1,1]*a[2,2]*a[3,3]+a[2,1]*a[3,2]*a[1,3]+
        a[3,1]*a[1,2]*a[2,3]-a[3,1]*a[2,2]*a[1,3]-a[2,1]*a[1,2]*a[3,3]-
        a[1,1]*a[3,2]*a[2,3];
      end;
 
 
   procedure GetMatr(a:mat4; var b:mat3; m:integer);
    { Вычеркивание из матрицы m-го столбца и первой строки, процедура нужнаая
    при нахождении определителя 4-го порядка }
    var
      i,j:integer;
    begin
      for i:=2 to 4 do
        for j:= 1 to 4 do
          begin
            if j=m
              then
                  continue
              else
                if j<m
                  then
                    b[i-1,j]:=a[i,j]
                  else
                    b[i-1,j-1]:=a[i,j];
          end;
    end;
 
    {определитель 4-го порядка}
    function determine4(A:mat4):real;
      var
        i,st:integer;
        tA:mat3;
        S:real;
      begin
        st:=1;
        s:=0;
        for i:= 1 to 4 do
          begin
            GetMatr(a,tA,i);
            s:=s+st*determine3(tA)*a[1,i];
            st:=st*(-1);
          end;
        determine4:=s;
      end;
 
 
 
begin
  { TODO -oUser -cConsole Main : Insert code here }
  {количество точек}
  writeln('Введите количество точек');
  readln(N);
 
  {ввод координат}
  for i:= 1 to N do
    begin
      writeln('Введите координаты  ',i,'-й вершины (x,y)');
      readln(p[i,1],p[i,2]);
    end;
 
    {инициализация переменных, точнее просто заполнение нулями}
    vnutri:=0;
    vne:=N;
    vnutriTemp:=0;
    vneTemp:=0;
    x0:=0;
    r:=0;
    y0:=0;
    square:=0;
    {перебор всех точек}
    for i:= 1 to N do
      for j:= 1 to N do
        for k:= 1 to N do
          {условие чтоб точки не совпадали}
          if (i<>j) and (i<>k) and (j<>k) then
            begin
             {нахождение площади и радиуса}
             a:=dlina(p[i],p[j]);
             b:=dlina(p[k],p[j]);
             c:=dlina(p[k],p[i]);
             pol:=(a+b+c)/2;
             pk:=pol*(pol-a)*(pol-b)*(pol-c);
             if pk<=0 then continue else square:=sqrt(pk);
             r:=(a*b*c)/(4*square);{радиус окружности}
              {цикл для перебора всех точек кроме точек треугольника}
              for l:= 1 to N do
                {четвертая точка чтоб не совпадала с точками треугольника}
                if  (l<>i) and (l<>j) and (l<>k)  then
                 begin
                 {заполнение матрицы}
                  matr4[1,1]:=sqr(p[l,1])+sqr(p[l,2]); matr4[1,2]:=p[l,1]; matr4[1,3]:=p[l,2]; matr4[1,4]:=1;
                  matr4[2,1]:=sqr(p[i,1])+sqr(p[i,2]); matr4[2,2]:=p[i,1]; matr4[2,3]:=p[i,2]; matr4[2,4]:=1;
                  matr4[3,1]:=sqr(p[j,1])+sqr(p[j,2]); matr4[3,2]:=p[j,1]; matr4[3,3]:=p[j,2]; matr4[3,4]:=1;
                  matr4[4,1]:=sqr(p[k,1])+sqr(p[k,2]); matr4[4,2]:=p[k,1]; matr4[4,3]:=p[k,2]; matr4[4,4]:=1;
                  det:=determine4(matr4);
                  {находим определитель и сравниваем знак}
                  if det>0 then inc(vneTemp) else inc(vnutriTemp);
                 end;
             {если при вершинах треугольника i,j,k разность точек минимальна
             то находим центр окрцжности }
             if abs(vne-vnutri)>abs(vneTemp-vnutriTemp)
                  then
                    begin
                      {нахождение центра окружности}
                      matr3[1,1]:=sqr(p[i,1])+sqr(p[i,2]); matr3[1,2]:=p[i,2]; matr3[1,3]:=1;
                      matr3[2,1]:=sqr(p[j,1])+sqr(p[j,2]); matr3[2,2]:=p[j,2]; matr3[2,3]:=1;
                      matr3[3,1]:=sqr(p[k,1])+sqr(p[k,2]); matr3[3,2]:=p[k,2]; matr3[3,3]:=1;
                      x0:=determine3(matr3)/(4*square);
                      matr3[1,2]:=p[i,1]; matr3[2,2]:=p[j,1]; matr3[3,2]:=p[k,1];
                      y0:=(-1)*(determine3(matr3))/(4*square);
                      vnutri:=vnutriTemp;
                      vne:=vneTemp;
                    end;
                vnutriTemp:=0;
                vneTemp:=0;
                end;
      {вывод данных}
      writeln('x0=',x0:5:3);
      writeln('y0=',y0:5:3);
      writeln('r=',r:5:3);
     readln;
end.