procedure tform1.third(a1,a2,a3,a4:extended; var y:otvet);
var  Q, R, t, A, B: extended;
begin
 // --> x^3 + a2*x^2 + a3*x + a4 = 0
  a2 := a2/a1;   a3 := a3/a1;   a4 := a4/a1;

  Q := (sqr(a2) - 3*a3) / 9;
  R := (a2*(2*sqr(a2) - 9*a3) + 27*a4) / 54;

    // --> 3 дійсних корені.
  if sqr(R) < sqr(Q)*Q
  then begin
    t := ArcCos(R/sqrt(sqr(Q)*Q)) / 3;
    y[1].x:= -2*sqrt(Q)*cos(t) - a2/3;            y[1].y:= 0;
    y[2].x:= -2*sqrt(Q)*cos(t+(2*pi/3)) - a2/3;   y[2].y:= 0;
    y[3].x:= -2*sqrt(Q)*cos(t-(2*pi/3)) - a2/3;   y[3].y:= 0;
  end
  else begin
    A := - Power(Abs(R)+sqrt(sqr(R)-sqr(Q)*Q), 1/3);
    if R < 0 then A := - A;

    if A <> 0
    then B := Q / A
    else B := 0;

    y[1].x:=A + B - a2/3;    y[1].y:= 0;  // 1 дійсний корінь
    if A <> B
    then begin
     // ще 1 комплексний корінь
      y[2].x:= -(A+B)/2 - a2/3;     y[2].y:= - sqrt(3)*abs(A-B)/2;
    end
    else begin
     // ще 1 дйсний корінь
     y[2].x:= - A - a2/3;    y[2].y:= 0;
    end;
     // ще 1 комплексйни корінь або дублювання дійсного
     y[3].x:=y[2].x;    y[3].y:= -y[2].y;
  end;

end;

procedure tform1.fourth(a1,a2,a3,a4,a5:extended; var f:otvet);
var
  p, q, r: extended;
  z,z3:otvet;
  a2div4,prov,prov1: TComplex;
begin
   // --> x^4 + a2*x^3 + a3*x^2 + a4*x + a5 = 0
  a2 := a2/a1;      {a1 := 1;}
  a3 := a3/a1;
  a4 := a4/a1;
  a5 := a5/a1;
               a2div4.x:=- a2 / 4;      a2div4.y:=0;
   // --> x = y - a2/4
  // --> y^4 + p*y^2 + q*y + r = 0

  p := a3 - 3*sqr(a2)/8;
  q := a4 - a2*a3/2 + a2*a2*a2/8;
  r := sqr(a2)*a3/16 - 3*sqr(sqr(a2))/256 - a2*a4/4 + a5;
   a1:=1; a2:=p/2; a3:=(sqr(p)-4*r)/16; a4:=-sqr(q)/64;

// --> z^3 + p/2*z^2 + (p^2-4*r)/16*z - q^2/64 = 0

  Third(a1,a2,a3,a4, z);
for i:=1 to 3 do
       z3[i]:=csqrt(z[i]);

 prov1.x:=z3[1].x*z3[2].x-z3[1].y*z3[2].y;
 prov1.y:=z3[1].x*z3[2].y+z3[2].x*z3[1].y;
 prov.x:=prov1.x*z3[3].x-prov1.y*z3[3].y;
 prov.y:=prov1.x*z3[3].y+prov1.y*z3[3].x;

  //   y = sqrt(z1)+sqrt(z2)+sqrt(z3),
  // а x = sqrt(z1)+sqrt(z2)+sqrt(z3) - a2/4,
  //  коли sqrt(z1)*sqrt(z2)*sqrt(z3) = - q / 8;


if ( -( q / 8  )>0) and (prov.x>0)
  then begin
      // +++
     f[1].x:=z3[1].x+z3[2].x+z3[3].x+a2div4.x;
     f[1].y:=z3[1].y+z3[2].y+z3[3].y+a2div4.y;
      // +--
     f[2].x:=z3[1].x-z3[2].x-z3[3].x+a2div4.x;
     f[2].y:=z3[1].y-z3[2].y-z3[3].y+a2div4.y;
     // -+-
     f[3].x:=-z3[1].x+z3[2].x-z3[3].x+a2div4.x;
     f[3].y:=-z3[1].y+z3[2].y-z3[3].y+a2div4.y;
     // --+
     f[4].x:=-z3[1].x-z3[2].x+z3[3].x+a2div4.x;
     f[4].y:=-z3[1].y-z3[2].y+z3[3].y+a2div4.y;
  end
  else begin
     // ++-
     f[1].x:=z3[1].x+z3[2].x-z3[3].x+a2div4.x;
     f[1].y:=z3[1].y+z3[2].y-z3[3].y+a2div4.y;
     // +-+
     f[2].x:=z3[1].x-z3[2].x+z3[3].x+a2div4.x;
     f[2].y:=z3[1].y-z3[2].y+z3[3].y+a2div4.y;
     // -++
     f[3].x:=-z3[1].x+z3[2].x+z3[3].x+a2div4.x;
     f[3].y:=-z3[1].y+z3[2].y+z3[3].y+a2div4.y;
     // ---
     f[4].x:=-z3[1].x-z3[2].x-z3[3].x+a2div4.x;
     f[4].y:=-z3[1].y-z3[2].y-z3[3].y+a2div4.y;
   end;
end;