Finding the roots of a polynomial

Question

Given a polynomial of degree n, how can I reliably find all values of x where y = 0.

I'm currently using the Math.Polynomial library, which doesn't seem to have this function built in. However, I could be missing something here, it seems like this would be a widely used function.

Thanks

`separateRoots` *should* get you most of the way there, but it doesn't seem to work correctly. E.g. `separateRoots BE [1,0,-1]` fails to separate. — luqui, Apr 11 '18 at 02:49
Note that once *n* is greater than 5, there is no algebraic exact way to find the roots: https://en.wikipedia.org/wiki/Abel%E2%80%93Ruffini_theorem — Willem Van Onsem, Apr 11 '18 at 08:21
@WillemVanOnsem Interesting. I was sure this theorem was originally due to Galois. — Stéphane Laurent, Apr 11 '18 at 09:35
Well the theorem has been proved in 1824, at that moment, Galois was 13 years old. Although Galois was very young when he died (as a result of a love affair and in a duel if I recall correctly), but as far as I know, he did not wrote that theorem when he was 13. — Willem Van Onsem, Apr 11 '18 at 09:41
@WillemVanOnsem No general way, there however exceptions, x^5 = 1 has only roots which can be algebraically found for example. — Elmex80s, Apr 11 '18 at 10:45
@luqui thats what I was thinking, but wasn't sure if that was expected behavior or not. I may file a bug report. — , Apr 11 '18 at 15:16
@WillemVanOnsem I read somewhere (don't remember where) that Abel proved the case degree=5 only (this contradicts Wikipedia), while Galois proved the case higher than 5. — Stéphane Laurent, Jul 24 '18 at 09:43

alias · Answer 1 · 2018-04-11T06:08:46.773

If you don't mind using an external SMT solver, you can use the SBV package:

Prelude Data.SBV> allSat $ \x -> 2*x^3-19*x^2+15*x+72 .== (0::SReal)
Solution #1:
  s0 = 3.0 :: Real
Solution #2:
  s0 = 8.0 :: Real
Solution #3:
  s0 = -1.5 :: Real
Found 3 different solutions.

You are guaranteed that the roots will be precise, i.e., no rounding error will happen. The Real type corresponds to infinitely precise algebraic reals. If the result cannot be printed in a finite form, it'll be written as the root of a simple equation and will be expanded to the current print precision:

Prelude Data.SBV> allSat $ \x -> x^2-2 .== (0::SReal)
Solution #1:
  s0 = root(1, x^2 = 2) = -1.414213562373095... :: Real
Solution #2:
  s0 = root(2, x^2 = 2) = 1.414213562373095... :: Real
Found 2 different solutions.

The number of digits can be arbitrarily long and is configurable:

Prelude Data.SBV> allSatWith z3{printRealPrec=20} $ \x -> x^2-2 .== (0::SReal)
Solution #1:
  s0 = root(1, x^2 = 2) = -1.4142135623730950488... :: Real
Solution #2:
  s0 = root(2, x^2 = 2) = 1.4142135623730950488... :: Real
Found 2 different solutions.

Note that SMT solvers will only get you real-roots; no complex solutions will be found:

Prelude Data.SBV> allSat $ \x -> x^2+1 .== (0::SReal)
No solutions found.

If only some roots are real, those will be found:

Prelude Data.SBV> allSat $ \x -> x^3+2*x-1 .== (0::SReal)
Solution #1:
  s0 = root(1, x^3+2x = 1) = 0.4533976515164037... :: Real
This is the only solution.

PS. To make this all work, do not forget to first install Z3 on your machine. You can get a recent copy from https://github.com/Z3Prover/z3/releases

score 3 · Accepted Answer · answered Apr 11 '18 at 09:26

You can use the Polynomial.Roots module of the dsp package. It deals with complex values (complex coefficients and complex roots).

For instance, to solve 1+2x+2x²=0:

Prelude> import Polynomial.Roots 
Prelude Polynomial.Roots> roots 1e-16 1000 [1,2,2]
[(-0.5) :+ (-0.5),(-0.5) :+ 0.5]

As a sanity check, R gives the same result:

> polyroot(c(1,2,2))
[1] -0.5+0.5i -0.5-0.5i

Finding the roots of a polynomial

2 Answers2