{ Центр вписанной в тетраэдр сферы }
type  vector = array[1..3] of real;
{-----------------------------------------------------------------}
procedure vpro(a,b: vector; var c: vector);
begin
  c[1]:=a[2]*b[3] - a[3]*b[2];
  c[2]:=a[3]*b[1] - a[1]*b[3];
  c[3]:=a[1]*b[2] - a[2]*b[1];
  end;
{-----------------------------------------------------------------}
function spro(a,b: vector): real;
begin    spro:=a[1]*b[1] + a[2]*b[2] + a[3]*b[3];    end;
{-----------------------------------------------------------------}
procedure norm(var a: vector);
{ Нормирует вектор  a }
var  n: real;
begin           n:=sqrt(spro(a,a));
  a[1]:=a[1]/n;    a[2]:=a[2]/n;    a[3]:=a[3]/n;    end;
{-----------------------------------------------------------------}
procedure vadd(a,b: vector; t: real; var c: vector);
{ c = a + t*b }
begin
  c[1]:=a[1]+t*b[1];    c[2]:=a[2]+t*b[2];    c[3]:=a[3]+t*b[3];    end;
{-----------------------------------------------------------------}
procedure nbiss(a,b, c,d: vector; var n: vector);
{ Нормаль биссектрисы двугранного угла, проходящего через ребро  AB }
var  n1,n2, ab,bc,bd: vector;
begin
  vadd(b,a,-1, ab);    vadd(c,b,-1, bc);    vadd(d,b,-1, bd);
  vpro(ab,bc, n1);    norm(n1);
  vpro(ab,bd, n2);    norm(n2);
  vadd(n1,n2,1, n);     end;
{-----------------------------------------------------------------}
function dist(a,b,c, x: vector): real;
{ Расстояние от точки  x  до плоскости, содержащей точки a,b,c }
var  n, ab,ac,ax: vector;
begin
  vadd(b,a,-1, ab);    vadd(c,a,-1, ac);    vadd(x,a,-1, ax);
  vpro(ab,ac, n);      norm(n);
  dist:=spro(ax,n);    end;
{-----------------------------------------------------------------}
{-----------------------------------------------------------------}
function det(a1,a2,a3: vector): real;
var  x: vector;
begin    vpro(a2,a3, x);    det:=spro(a1,x);    end;
{-----------------------------------------------------------------}
procedure kramer(a1,a2,a3, f: vector; var x: vector);
{ Метод Крамера:  a1,a2,a3 - строки системы,  f - правые части }
var  d: real;    b1,b2,b3: vector;
begin
  d:=det(a1,a2,a3);
  b1:=a1;    b2:=a2;    b3:=a3;
  b1[1]:=f[1];    b2[1]:=f[2];    b3[1]:=f[3];
  x[1]:=det(b1,b2,b3)/d;
  b1:=a1;    b2:=a2;    b3:=a3;
  b1[2]:=f[1];    b2[2]:=f[2];    b3[2]:=f[3];
  x[2]:=det(b1,b2,b3)/d;
  b1:=a1;    b2:=a2;    b3:=a3;
  b1[3]:=f[1];    b2[3]:=f[2];    b3[3]:=f[3];
  x[3]:=det(b1,b2,b3)/d;
  end;
{-----------------------------------------------------------------}
{-----------------------------------------------------------------}
procedure vrand(var x: vector);
begin    x[1]:=random;    x[2]:=random;    x[3]:=random;    end;
{-----------------------------------------------------------------}
{=================================================================}
var  a,b,c,d, x, a1,a2,a3, f: vector;
{=================================================================}
begin
  vrand(a);    vrand(b);    vrand(c);    vrand(d);
  nbiss(a,b,c,d, a1);    nbiss(b,c,d,a, a2);    nbiss(c,d,a,b, a3);
  f[1]:=spro(a1,a);      f[2]:=spro(a2,b);      f[3]:=spro(a3,c);
  kramer(a1,a2,a3, f, x);
  writeln('=====================================');
  writeln(dist(a,b,c, x):21:15, dist(a,b,c, d):21:15);
  writeln(dist(a,b,d, x):21:15, dist(a,b,d, c):21:15);
  writeln(dist(a,c,d, x):21:15, dist(a,c,d, b):21:15);
  writeln(dist(b,c,d, x):21:15, dist(b,c,d, a):21:15);
  readln;
  end.