CODEWITHSUNDEEP
pythonsympy
Mathita314
Mathita314

Reputation: 1

Sympy real roots

Anyone know how the real_roots function works in sympy? Because on the official page I don't find how sympy evaluates the solution. In particular, I am only interested in the solution of a polynomial function

Upvotes: 0

Views: 1522

Answers (1)

Oscar Benjamin
Oscar Benjamin

Reputation: 14470

The real_roots function supports finding the real roots of an univariate polynomial with rational coefficients e.g.:

In [12]: real_roots(x**2 - 1)
Out[12]: [-1, 1]

In [13]: real_roots(x**2 + 1)
Out[13]: []

Internally this is done in several stages:

  1. The polynomial can be completely factorised over the rationals into irreducibles.
  2. For each irreducible expressions for the roots are computed as either rational, radical or CRootOf.
  3. Root isolation algorithms are used to identify which of those roots is real.

Note that step 1 immediately identifies all rational roots since they are reduced to linear factors. In step 2 the cubic and quartic formulas are not used but the quadratic formula is and decomposition and other techniques are possibly tried (not sure). You can use the roots function for the cubic and quartic formulas but they usually produce excessively complicated output. Anything not handled becomes CRootOf.

In step 3 the implementation can work for CRootOf meaning that it can also work for quintics and higher degree polynomials (it is not limited by Abel-Ruffini):

In [6]: [r] = real_roots(x**5 - x - 1)

In [7]: r
Out[7]: 
       ⎛ 5           ⎞
CRootOf⎝x  - x - 1, 0⎠

In [8]: r.evalf()
Out[8]: 1.16730397826142

The root isolation code is here: https://github.com/sympy/sympy/blob/01a5776b9feaf01d06c94165ccebe4ac4b617a7b/sympy/polys/rootisolation.py

That code implements some fairly complex algorithms that are apparently quite fast but the implementation in sympy is very slow for higher degree polynomials. I am not sure if those algorithms are optimal for something like univariate polynomials or whether the implementation is just not well optimised. I think that an implementation based on the interval-Newton method could be much faster:

https://github.com/sympy/sympy/issues/19164

Upvotes: 2

Related Questions