f[x_] := x^3 - 2 x + 1;
Solve[f[x] == 0, x]
{{x -> 1}, {x -> 1/2 (-1 - Sqrt[5])}, {x -> 1/2 (-1 + Sqrt[5])}}

But the 2nd element of the list is not a solution in fact:

N[f[1/2 (-1 + Sqrt[5])]]

Reduce has the same issue, so how to force both of those functions to only return exact solutions as -5.55112*10^-17 is only an approximation but is not a solution.

Can Solve and Reduce be trusted, and what is the real meaning of that solution 1/2 (-1 + Sqrt[5])?

Is there a way to ask Mathematica to provide a detailed explanation of how it got to that approximation of a solution (kind of the reverse tracing of computations that occurred)?

  • 1
    $\begingroup$ Please read the guidelines to the bugs tag. $\endgroup$ – Michael E2 Apr 15 '19 at 2:57
  • 4
    $\begingroup$ Have you heard of round-off error? Floating-point numbers are discrete and cannot exactly represent most numbers (e.g. Sqrt[5] and 1/7 and so forth). Try N[f[1/2 (-1 + Sqrt[5])], 100] and Simplify[f[1/2 (-1 + Sqrt[5])]] $\endgroup$ – Michael E2 Apr 15 '19 at 2:59
  • 1
    $\begingroup$ It will be impossible to have a bug in Solve for a third order polynomial, since there is an exact formula for the roots of 3rd order and it has been known for around 600 years or so. Found by Girolamo Cardano somewhere in the 1500's. $\endgroup$ – Nasser Apr 15 '19 at 5:01
  • $\begingroup$ I get the technical reason now, so about the comment on question being off-topic, I tried to delete the question all together but the system refuses because a formal answer has been provided. If you can delete it as an admin of the site, feel free to do so. I get the feeling it is not a site for people that are slowly returning to Mathematics after 25 years and new to Mathematica. All the best to you. $\endgroup$ – user3851217 Apr 15 '19 at 22:32

No you have not found any bugs. The solutions Solve returned are exact. The problem you ran into lies with using N in a naive way and not with Solve or Reduce.

The following code shows the solutions are exact.

f[x_] := x^3 - 2 x + 1
sol = x /. Solve[f[x] == 0, x] // Flatten;
Simplify[f /@ sol]

{0, 0, 0}


Not the answer you're looking for? Browse other questions tagged or ask your own question.