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!