Giải phương trình bậc 3 có cả nghiệm phức

Từ rất lâu rồi tôi đã muốn tự giải được phương trình bậc 3 và sau 1 thời gian tôi rất vui được giới thiệu đến mọi người đoạn mã pascal của chương trình giải phương trình bậc 3 có cả nghiệm số phức. Do kiến thức hạn hẹp nên có nhiều chỗ không được tốt lắm mong bạn đọc sửa chữa cho tốt hơn.
program giaiptbac3;
uses crt;
function arccos(k:real):real; {tao ham arrcos khi |k|<=1}
var alpha:real;
begin
if k>0 then alpha:=arctan(sqrt(1/k/k-1)) else if k<0 then alpha:=pi+arctan(-sqrt(1/k/k-1)) else> {tan2=1/cos2-1}
arccos:=alpha;
end;
function can(k:real):real; {tao ham tinh can bac 3 thong qua exp va ln}
var ex:real;
begin
if k<0 then ex:=-exp(ln(abs(k))/3) else if k>0 then ex:=exp(ln(k)/3) else ex:=0;
{vi a^x=(e^ln a)^x = e^(ln a . x) a>0}
can:=ex;
end;
procedure nghiemphuc(a,b,c:real); {tao thu tuc tim nghiem phuc}
var d:real; x,y:real;
begin
d:=-b*b+4*a*c; {delta luon nho hon 0 trong truong hop nay nen lay so doi de khong can dung abs(d)}
x:=-b/2/a; y:=sqrt(d)/2/a;
writeln(` Nghiem Phuc: x2 = `,x,` + `,y,`i`);
writeln(` x3 = `,x,` - `,y,`i`);
end;
var a,b,c,d:real; s,k:real;
x1,x2,x3:real;
begin
clrscr;
write(`Nhap a,b,c,d = `);readln(a,b,c,d);
s:=b*b-3*a*c; writeln(s);
if s<>0 then
begin
k:=(9*a*b*c-2*b*b*b-27*a*a*d)/2/sqrt(abs(s*s*s)); writeln(k);
if s>0 then
begin
if abs(k)<=1 then begin
writeln(`PT co 3 nghiem thuc: `);
x1:=(2*sqrt(s)*cos(arccos(k)/3)-b)/(3*a);
x2:=(2*sqrt(s)*cos(arccos(k)/3-2*pi/3)-b)/(3*a);
x3:=(2*sqrt(s)*cos(arccos(k)/3+2*pi/3)-b)/(3*a);
writeln(`x1 = `,x1:10:12);
writeln(`x2 = `,x2:10:12);
writeln(`x3 = `,x3:10:12);
end
else
begin
x1:=sqrt(abs(s))*abs(k)/3/a/k*(can(abs(k)+sqrt(k*k-1))+can(abs(k)-sqrt(k*k-1)))-b/3/a;
writeln(` Nghiem thuc: x1 = `,x1:10:12);
a:=a; b:=x1*a+b; c:=x1*b+c; {mieu ta hooclo}
nghiemphuc(a,b,c);
end;
end
else
begin
x1:=sqrt(abs(s))/(3*a)*(can(k+sqrt(k*k+1))+can(k-sqrt(k*k+1)))-b/3/a;
writeln(` Nghiem thuc: x1 = `,x1:10:12);
a:=a; b:=x1*a+b; c:=x1*b+c; {mieu ta hooc lo}
nghiemphuc(a,b,c);
end;
end
else
begin
x1:=(-b+can(b*b*b-27*a*a*d))/3/a;
writeln(` Nghiem thuc: x = `,x1:10:12);
end