Sympy solve of cubic equation appears to be wrong and includes imaginary units. What is wrong?

Question

from sympy import *
from sympy.abc import x

f = x**3 - 3*x + 1
res = solve(f)
print(res)

However, it gives an incorrect answer with imaginary units. Checking via subs doesn't return 0, either.

UPD:

res:

[-3/((-1/2 - sqrt(3)*I/2)*(27/2 + 27*sqrt(3)*I/2)**(1/3)) - (-1/2 - sqrt(3)*I/2)*(27/2 + 27*sqrt(3)*I/2)**(1/3)/3,
 -(-1/2 + sqrt(3)*I/2)*(27/2 + 27*sqrt(3)*I/2)**(1/3)/3 - 3/((-1/2 + sqrt(3)*I/2)*(27/2 + 27*sqrt(3)*I/2)**(1/3)),
 -(27/2 + 27*sqrt(3)*I/2)**(1/3)/3 - 3/(27/2 + 27*sqrt(3)*I/2)**(1/3)]

subs check:

>>> print(f.subs(x, res[0]))
1 + (-1/2 - sqrt(3)*I/2)*(27/2 + 27*sqrt(3)*I/2)**(1/3) + (-3/((-1/2 - sqrt(3)*I/2)*(27/2 + 27*sqrt(3)*I/2)**(1/3)) - (-1/2 - sqrt(3)*I/2)*(27/2 + 27*sqrt(3)*I/2)**(1/3)/3)**3 + 9/((-1/2 - sqrt(3)*I/2)*(27/2 + 27*sqrt(3)*I/2)**(1/3))

The plot:

>>> plot(f, xlim=(-2, 2), ylim=(-2, 2))

Answer 1

Let's try to convert the symbolic solution computed by sympy to floating point numbers:

from sympy import *
var("x")
f = x**3 - 3*x + 1
sol = solve(f)
print([s.n() for s in sol])
# [0.347296355333861 - 0.e-23*I, 1.53208888623796 + 0.e-20*I, -1.87938524157182 + 0.e-23*I]

Here you can see very small imaginary parts. These are rounding errors. We can remove them with this:

print([s.n(chop=True) for s in sol])
# [0.347296355333861, 1.53208888623796, -1.87938524157182]

Now, I'm going to substitute the first symbolic solution back into the equation:

print(f.subs(x, sol[0]))
# 1 + (-1/2 - sqrt(3)*I/2)*(27/2 + 27*sqrt(3)*I/2)**(1/3) + (-3/((-1/2 - sqrt(3)*I/2)*(27/2 + 27*sqrt(3)*I/2)**(1/3)) - (-1/2 - sqrt(3)*I/2)*(27/2 + 27*sqrt(3)*I/2)**(1/3)/3)**3 + 9/((-1/2 - sqrt(3)*I/2)*(27/2 + 27*sqrt(3)*I/2)**(1/3))

It doesn't appear to be zero, but we can ask sympy to simplify it, like this:

print(f.subs(x, sol[0]).simplify())
# 0

As you can see, the symbolic solutions computed by sympy are correct.