4

I'm trying to use Haskell for evaluating a few functions. One of them is:

k x = 27 * x ^ 4 + 1792 * x ^ 3 + 2592 * x ^ 2 - 324 * x - 243

Using wolfram alpha I've found it's root:

27x^4 + 1792x^3 + 2592x^2 - 324x - 243 = 0
=> x ≃ 0.331478863194230

But haskell gives the wrong answer when evaluated with the root:

k 0.331478863194230 == -6.252776074688882e-13

I'm sure it's not Wolfram's issue, as Desmos gives the same result. I'm not that familiar with Haskell and I wonder what might be the cause of this quirky behavior.

Tomasz Surowiec
  • 49
  • 1
  • 5
    Looks very close to zero to me. Perhaps a rounding or floating point error? – DarthFennec Oct 04 '19 at 16:55
  • 7
    Quickly, before someone recommends that you familiarise yourself with *What Every Computer Scientist Should Know About Floating Point Arithmetic*, do some research on how floating point numbers work and why you might think that `-6.252776074688882e-13` is a good approximation to `0`. – High Performance Mark Oct 04 '19 at 16:55
  • 4
    Also note that Wolfram tells you an approximate solution for `x ≃ 0.331478863194230` – lehins Oct 04 '19 at 16:57
  • 4
    And this is a good place to start your research ... [Is floating point math broken?](https://stackoverflow.com/questions/588004/is-floating-point-math-broken) – High Performance Mark Oct 04 '19 at 17:00
  • `x ≃ 0.331478863194230` Since this is an approximation of a root of the polynomial, you shouldn't be surprised that `f(x)` is also *approximately* zero. Note the `e-13` on the end. This means that the answer is on the order of `10^-13`, so very small. – Code-Apprentice Oct 04 '19 at 17:54
  • 2
    What if __inside Wolfram__ you ask back for the value of the polynomial function at its _approximate_ root ? Does Wolfram tell you it is exactly zero ?? – jpmarinier Oct 04 '19 at 17:57
  • Using [Horner's Method](https://en.wikipedia.org/wiki/Horner's_method) (using `k x = (((27*x + 1792)*x + 2592)*x - 324)*x - 243)`) will get you the slightly more accurate answer `-5.684341886080801e-13`. – Daniel Wagner Oct 04 '19 at 18:42
  • Using more digits the Horner formula gives 0.0 for the floating point evaluation at x=0.33147886319423027. The original monomial formula gives there the value -2.84e-14 while a little to the right x=0.33147886319423033 the floating point evaluation gives 5.68e-14. At about x=0.331478863194230283 the value jumps from the lower to the higher. – Lutz Lehmann Oct 08 '19 at 10:42

1 Answers1

6

By the Rational Root Theorem, it's easy to see this polynomial has no rational roots. So Wolfram would never be able to show you an exact answer in decimal form. It of course has to round to a near approximation, and even indicates it's doing so with the ≃ symbol. It's not surprising that plugging that approximate root into your function doesn't output 0.

In fact, if you avoid floating point numbers you can see with arbitrary precision what the output would be when passing that input:

 > k x = 27 * x ^ 4 + 1792 * x ^ 3 + 2592 * x ^ 2 - 324 * x - 243
 > k 0.331478863194230 :: Rational
=> (-54003681976002235123002052941041334918065893) % 100000000000000000000000000000000000000000000000000000000

which is

-5.4003681976002235123002052941041334918065893e-13
Amit Kumar Gupta
  • 17,184
  • 7
  • 46
  • 64
  • You could go even further than `Rational` with `CReal` here, I think, if you wanted to. – Rein Henrichs Oct 05 '19 at 16:32
  • @Daniel Wagner could you share some complete code that I could put into https://repl.it/languages/haskell or https://www.tryhaskell.org/, I have not been able to get your edits to my answer to work, and not sure how to respond to Rein Henrichs suggestion. – Amit Kumar Gupta Oct 05 '19 at 17:44
  • @AmitKumarGupta The only additional line of code needed is from the question itself: `k x = 27 * x ^ 4 + 1792 * x ^ 3 + 2592 * x ^ 2 - 324 * x - 243`. (No modifications are needed for this to support both `Double` and `Rational`.) I think `CReal` is not a super interesting next step from `Rational`, personally; I suspect it was suggested because it is another type for arbitrary-precision calculations, but `Rational` already gives an exact answer, so why bother with the extra complications/slowness of `CReal`? – Daniel Wagner Oct 05 '19 at 21:22