program bug;

{$mode objfpc}{$H+}

uses
  Classes, SysUtils, math;
  { you can add units after this }

type tpointarray=array of tpoint;
procedure smallestCircel(var x,y,rsqr: float; const p1i,p2i,p3i: longint; const points: tpointarray);
var //p1i,p2i,p3i: longint;
    //hull: TLongintArray;
    f1,f2,besti,temp:longint;
    nx,ny,nrsqr:float;
    p0,p1,p2: ^tpoint;
begin
  //(x-x0)^2+(y-y0)^2=r^2=x^2-2xx0+x0^2+y^2-2yy0+y0^2
  //(x-x1)^2+(y-y1)^2=r^2=x^2-2xx1+x1^2+y^2-2yy1+y1^2
  //(x-x2)^2+(y-y2)^2=r^2=x^2-2xx2+x2^2+y^2-2yy2+y2^2
  //=> r^2 - r^2=x^2-2xx0+x0^2+y^2-2yy0+y0^2 - x^2+2xx1-x1^2-y^2+2yy1-y1^2
  //-2xx0+x0^2-2yy0+y0^2 +2xx1-x1^2 + 2yy1-y1^2 = 0
  //2x(x1-x0) + 2y(y1-y0) = x1^2 + y1^2 - x0^2 - y0^2
  //2x(x2-x1) + 2y(y2-y1) = x2^2 + y2^2 - x1^2 - y1^2
  //p1i:=hull[i];
  p0:=@points[p1i];
  //p2i:=hull[j];
  p1:=@points[p2i];
  //p3i:=hull[k];
  p2:=@points[p3i];
  nrsqr:=-1;
  if p0^.y=p1^.y then begin
    if (p1^.y=p2^.y) or (p0^.x=p1^.x) then begin
      //Fall 0: points on a line
      ny:=p0^.y;
      if (p0^.x<=p1^.x) and (p1^.x<=p2^.x) then begin
        nx:=(p0^.x + p2^.x) / 2;
        nrsqr:=sqr(p0^.x - p2^.x) /4;
      end else if (p1^.x<=p0^.x) and (p0^.x<=p2^.x) then begin
        nx:=(p1^.x + p2^.x) / 2;
        nrsqr:=sqr(p1^.x - p2^.x) /4;
      end else begin//if (p0^.x<p2^.x) and (p2^.x<p1^.x) then begin
        nx:=(p0^.x + p1^.x) / 2;
        nrsqr:=sqr(p0^.x - p1^.x) /4;
      end;
    end else begin
      //=> Fall 1: y0=y1
      // 2x(x1-x0) = x1^2 + x0^2 => x = x1^2 + x0^2 / (2(x1-x0))
      nx:=(sqr(p1^.x) - sqr(p0^.x)) / (2*(p1^.x-p0^.x));
      // 2y(y2-y1) = x2^2 + y2^2 - x1^2 - y1^2 - 2x(x2-x1)
      ny := (sqr(p2^.x)+sqr(p2^.y)-sqr(p1^.x)-sqr(p1^.y) - 2*nx*(p2^.x-p1^.x))/(2*(p2^.y-p1^.y));
    end;
  end else begin
    //=> Fall 3: y0<>y1
    //2x(x1-x0)*(x2-x1) + 2y(y1-y0)*(x2-x1) = (x1^2 + y1^2 - x0^2 - y0^2) * (x2-x1)
    //2x(x1-x0)*(x2-x1) + 2y(y2-y1)*(x1-x0) = (x2^2 + y2^2 - x1^2 - y1^2) * (x1-x0)
    //=> 2y(y1-y0)*(x2-x1)-2y(y2-y1)*(x1-x0) = (x1^2 + y1^2 - x0^2 - y0^2) * (x2-x1)-(x2^2 + y2^2 - x1^2 - y1^2) * (x1-x0)
    //=> y= ((x1^2 + y1^2 - x0^2 - y0^2) * (x2-x1)-(x2^2 + y2^2 - x1^2 - y1^2) * (x1-x0)) /
    // (2(y1-y0)*(x2-x1)-2(y2-y1)*(x1-x0))
    temp:=(p1^.y-p0^.y)*(p2^.x-p1^.x)-(p2^.y-p1^.y)*(p1^.x-p0^.x);
    if temp=0 then begin
      //=>(p1^.y-p0^.y)/(p1^.x-p0^.x) = (p2^.y-p1^.y)/(p2^.x-p1^.x)
      //=> p0->p1 parallel to p1->p2
      //=> all points lie on a single line
      //where is p2 on the line from p0 to p1?
      //Eu: (u-p0^.x) * (p1^.y-p0^.y)/(p1^.x-p0^.x) + p0^.y = p2^.y
      //=> u = (p2^.y - p0^.y) * (p1^.x-p0^.x) / (p1^.y-p0^.y) + -p0^.x
      nx:=(p2^.y - p0^.y) * (p1^.x-p0^.x) / (p1^.y-p0^.y) -p0^.x; //exists, y checked above
      if nx > 1 then begin
        nx:=(p0^.x + p2^.x) / 2;
        ny:=(p0^.y + p2^.y) / 2;
        nrsqr:=sqr(p0^.x - p2^.x) /4;
      end else if nx<1 then begin
        nx:=(p1^.x + p2^.x) / 2;
        ny:=(p1^.y + p2^.y) / 2;
        nrsqr:=sqr(p1^.x - p2^.x) /4;
      end else begin
        nx:=(p0^.x + p1^.x) / 2;
        ny:=(p0^.y + p1^.y) / 2;
        nrsqr:=sqr(p0^.x - p1^.x) /4;
      end;
    end else begin
      ny:=((sqr(p1^.x) + sqr(p1^.y) - sqr(p0^.x) - sqr(p0^.y)) * (p2^.x-p1^.x) - (sqr(p2^.x) + sqr(p2^.y) - sqr(p1^.x) - sqr(p1^.y)) * (p1^.x-p0^.x)) / (2*temp);
      //2x(x1-x0) = x1^2 + y1^2 - x0^2 - y0^2 - 2y(y1-y0)
      //2x(x2-x1) = x2^2 + y2^2 - x1^2 - y1^2 - 2y(y2-y1)
      if p0^.x <> p1^.x then
        nx:= (sqr(p1^.x) + sqr(p1^.y) - sqr(p0^.x) - sqr(p0^.y) - 2*ny*(p1^.y-p0^.y)) / (2*(p1^.x-p0^.x))
      else
        nx:= (sqr(p2^.x) + sqr(p2^.y) - sqr(p1^.x) - sqr(p1^.y) - 2*ny*(p2^.y-p1^.y)) / (2*(p2^.x-p1^.x))
    end;
  end;
  if nrsqr < 0 then
    nrsqr := sqr(nx-p0^.x) + sqr(ny-p0^.y);

  x:=nx;
  y:=ny;
  rsqr:=nrsqr;
end;

var x,y,rsqr: float;
    points: tpointarray;
    s: string;
begin
  setlength(points,4);
  points[0].x:=55; points[0].y:=14;
  points[1].x:=63; points[1].y:=61;
  points[2].x:=10; points[2].y:=64;
  points[3].x:=20; points[3].y:=44;
  smallestCircel(x,y,rsqr,2,1,0,points);
  writestr(s,x:4:2,' ',y:4:2,' ',rsqr:4:2);
  if (s<>'35.31 41.53 1145.57') then
    halt(1);
end.