polynomials – Real solutions of third and fourth degree equations

A few hours ago I “discovered” that if a third or fourth degree equation has distinct real solutions, it’s possible to calculate them avoiding complex numbers.

In particular, we have:

poly = (x - 1)(x - 2)(x - 3);
{c, b, a} = CoefficientList(poly, x)((1 ;; 3));
d = a^2/3 - b;
e = 2 a^3/27 - a b/3 + c;
f = ArcCos(-3 Sqrt(3) e/(2 d Sqrt(d)));
N({-a/3 + 2 Sqrt(3 d) Cos((f - 4 Pi)/3)/3,
   -a/3 + 2 Sqrt(3 d) Cos((f - 2 Pi)/3)/3,
   -a/3 + 2 Sqrt(3 d) Cos((f - 0 Pi)/3)/3})

{1., 2., 3.}

poly = (x - 1)(x - 2)(x - 3)(x - 4);
{d, c, b, a} = CoefficientList(poly, x)((1 ;; 4));
e = 3 a^2/4 - 2 b;
f = 2 c - a b + a^3/4;
g = b^2 + 12 d - 3 a c;
h = 27 a^2 d - 9 a b c + 2 b^3 - 72 b d + 27 c^2;
i = ArcCos(h/(2 g Sqrt(g)));
j = (e + 2 Sqrt(g) Cos(i/3))/3;
N({-a/4 - (Sqrt(j) + Sqrt(e - j + f/Sqrt(j)))/2,
   -a/4 - (Sqrt(j) - Sqrt(e - j + f/Sqrt(j)))/2,
   -a/4 + (Sqrt(j) - Sqrt(e - j - f/Sqrt(j)))/2,
   -a/4 + (Sqrt(j) + Sqrt(e - j - f/Sqrt(j)))/2})

{1., 2., 3., 4.}

After a few moments of joy, however, I noticed that, for example, if I write (x - 10^-15) instead of (x - 1), I get respectively 6.66134*10^-16 and 8.88178*10^-16 instead of 10^-15 as the solution.

I intuitively believe that it’s a numerical problem and that the main cause is ArcCos, but I’m not too sure and also I’m not able to establish if something can be done to solve the issue, or if I have to give up and I must necessarily rely on it to the good Newton method.

Thank you!