Find the Roots of a Polynomial.

Updated 2 Tallied Votes 3K Views

This snippet shows how to find the complex roots of a polynomial using python (tested with python 2.6). You need the scipy or numpy module.

TrustyTony commented: Nice collection of libraries +13
``````from numpy import poly1d

P = poly1d([3.4, -5.22, 2, 1, -7.1], variable = 'x')

print("Here is a polynomial\n")
print(P)
print("\nIts complex roots as a numpy array:")
print(P.r)
print("\nIts complex roots as a python list:")
print(list(P.r))

"""my output -->
Here is a polynomial

4        3     2
3.4 x - 5.22 x + 2 x + 1 x - 7.1

Its complex roots as a numpy array:
[ 1.57667814+0.j          0.43601402+1.12245146j  0.43601402-1.12245146j
-0.91341205+0.j        ]

Its complex roots as a python list:
[(1.5766781395927345+0j), (0.43601401578206267+1.1224514554841667j), (0.43601401578206267-1.1224514554841667j), (-0.91341205350980392+0j)]
"""``````
Gribouillis 1,391

If you want roots to arbitrary precision, a very interesting library is rpncalc.ratfun (you will need another module for arbitrary precision called clnum). Here is an example with the same polynomial on a linux machine:

``````from rpncalc.ratfun import Polynomial

P = Polynomial(-7.1, 1, 2, -5.22, 3.4)
print("Here is a polynomial\n")
print P
print("\nIts complex roots as a list of multiprecision floating points numbers:\n")
print P.roots(eps=1e-30)

"""  my output ---->
Here is a polynomial

3.4*x^4-5.22*x^3+2*x^2+x-7.1

Its complex roots as a list of multiprecision floating points numbers:

[cmpf('-0.913412053509803965200954224858149742316-5.23858817606991616489822814712984590827e-312j',prec=36), cmpf('0.436014015782063442817853261862879518184-1.12245145548416568676706088778633149012j',prec=36), cmpf('0.436014015782063442817853261862879518184+1.12245145548416568676706088778633149012j',prec=36), mpf('1.57667813959273587005688087837686649466',36)]
"""``````
Gribouillis 1,391

If you have the sympy module, this is another way to find the roots. If the polynomial has rational coefficients and small degree, sympy can also compute exact roots. However I had some memory leak problems with sympy.

``````import sympy

print("Floating point roots:n")
print(sympy.roots([3.4, -5.22, 2, 1, -7.1]))

from sympy import Integer as Int

print("Exact roots as sympy expressions:")
print(sympy.roots([Int(34)/10, -Int(522)/100, 2, 1, -Int(71)/10]))

"""my output --->
Floating point roots:

{-0.913412053509803: 1, 0.436014015782064 + 1.12245145548416*I: 1, 0.436014015782064 - 1.12245145548416*I: 1, 1.57667813959273: 1}
Exact roots as sympy expressions:
{261/680 - (68363/346800 - 13501/(10404*(-38639957/424483200 + 12534291435179913031527229440000**(1/2)/6673555077120000)**(1/3)) + 2*(-38639957/424483200 + 12534291435179913031527229440000**(1/2)/6673555077120000)**(1/3))**(1/2)/2 - (68363/173400 + 11528419/(19652000*(68363/346800 - 13501/(10404*(-38639957/424483200 + 12534291435179913031527229440000**(1/2)/6673555077120000)**(1/3)) + 2*(-38639957/424483200 + 12534291435179913031527229440000**(1/2)/6673555077120000)**(1/3))**(1/2)) + 13501/(10404*(-38639957/424483200 + 12534291435179913031527229440000**(1/2)/6673555077120000)**(1/3)) - 2*(-38639957/424483200 + 12534291435179913031527229440000**(1/2)/6673555077120000)**(1/3))**(1/2)/2: 1, 261/680 + (68363/346800 - 13501/(10404*(-38639957/424483200 + 12534291435179913031527229440000**(1/2)/6673555077120000)**(1/3)) + 2*(-38639957/424483200 + 12534291435179913031527229440000**(1/2)/6673555077120000)**(1/3))**(1/2)/2 + (68363/173400 - 11528419/(19652000*(68363/346800 - 13501/(10404*(-38639957/424483200 + 12534291435179913031527229440000**(1/2)/6673555077120000)**(1/3)) + 2*(-38639957/424483200 + 12534291435179913031527229440000**(1/2)/6673555077120000)**(1/3))**(1/2)) + 13501/(10404*(-38639957/424483200 + 12534291435179913031527229440000**(1/2)/6673555077120000)**(1/3)) - 2*(-38639957/424483200 + 12534291435179913031527229440000**(1/2)/6673555077120000)**(1/3))**(1/2)/2: 1, 261/680 + (68363/173400 + 11528419/(19652000*(68363/346800 - 13501/(10404*(-38639957/424483200 + 12534291435179913031527229440000**(1/2)/6673555077120000)**(1/3)) + 2*(-38639957/424483200 + 12534291435179913031527229440000**(1/2)/6673555077120000)**(1/3))**(1/2)) + 13501/(10404*(-38639957/424483200 + 12534291435179913031527229440000**(1/2)/6673555077120000)**(1/3)) - 2*(-38639957/424483200 + 12534291435179913031527229440000**(1/2)/6673555077120000)**(1/3))**(1/2)/2 - (68363/346800 - 13501/(10404*(-38639957/424483200 + 12534291435179913031527229440000**(1/2)/6673555077120000)**(1/3)) + 2*(-38639957/424483200 + 12534291435179913031527229440000**(1/2)/6673555077120000)**(1/3))**(1/2)/2: 1, 261/680 + (68363/346800 - 13501/(10404*(-38639957/424483200 + 12534291435179913031527229440000**(1/2)/6673555077120000)**(1/3)) + 2*(-38639957/424483200 + 12534291435179913031527229440000**(1/2)/6673555077120000)**(1/3))**(1/2)/2 - (68363/173400 - 11528419/(19652000*(68363/346800 - 13501/(10404*(-38639957/424483200 + 12534291435179913031527229440000**(1/2)/6673555077120000)**(1/3)) + 2*(-38639957/424483200 + 12534291435179913031527229440000**(1/2)/6673555077120000)**(1/3))**(1/2)) + 13501/(10404*(-38639957/424483200 + 12534291435179913031527229440000**(1/2)/6673555077120000)**(1/3)) - 2*(-38639957/424483200 + 12534291435179913031527229440000**(1/2)/6673555077120000)**(1/3))**(1/2)/2: 1}
"""
``````
Gribouillis 1,391

You can also use the Gnu Scientific Library (http://www.gnu.org/software/gsl/). There is a python interface to gsl called pygsl http://pygsl.sourceforge.net/ . To install it you need an ansi C compiler and the development files of gsl and python-numeric. On a linux system, this is not a problem. Here is how to compute the roots using pygsl:

``````from pygsl.poly import poly_complex

L = [-7.1, 1, 2, -5.22, 3.4]
roots = poly_complex(len(L)).solve(L)

print("Using the Gnu Scientific Library\n")
print("Complex roots as a numpy array:")
print(roots)
print("\nComplex roots as a python list:")

""" my output ---->
Using the Gnu Scientific Library

Complex roots as a numpy array:
[-0.91341205+0.j          0.43601402+1.12245146j  0.43601402-1.12245146j
1.57667814+0.j        ]

Complex roots as a python list:
[(-0.91341205350980292+0j), (0.43601401578206328+1.122451455484166j), (0.43601401578206328-1.122451455484166j), (1.5766781395927358+0j)]
"""``````
Gribouillis 1,391

The PARI computer algebra system also has a function to compute the roots. To use it with python you need the python interface to pari http://code.google.com/p/pari-python/ . To install it on my Mandriva system, I had to install the packages pari and pari-devel first. It should work in Windows too with cygwin. Here are the roots of the same polynomial computed with pari

``````from pari import polroots
print polroots("-7.1 + x + 2*x^2 -5.22*x^3 + 3.4*x^4", 0)

""" My output --->
[-0.9134120535098039691852671888 + 0.E-38*I, 1.576678139592735906083229898 + 0.E-38*I, 0.4360140157820634433157245276 + 1.122451455484165684444607827*I, 0.4360140157820634433157245276 - 1.122451455484165684444607827*I]~
"""``````

Pari seems able to compute multiprecision roots too.

Gribouillis 1,391

Yet another module to find the roots: mpmath.polyroots() finds the roots with arbitrary precision and can also give an error estimate, using the Durand-Kerner method:

``````>>> from mpmath import *
>>> mp.dps = 60
>>> for r in polyroots([1, 0, -10, 0, 1]): # x^4 -10 x^2 + 1
...     print r
...
-3.14626436994197234232913506571557044551247712918732870123249
-0.317837245195782244725757617296174288373133378433432554879127
0.317837245195782244725757617296174288373133378433432554879127
3.14626436994197234232913506571557044551247712918732870123249
>>>
>>> sqrt(3) + sqrt(2)
3.14626436994197234232913506571557044551247712918732870123249
>>> sqrt(3) - sqrt(2)
0.317837245195782244725757617296174288373133378433432554879127
``````

It seems that `polyroot()`is also available in sympy as `sympy.mpmath.polyroot()`.

Gribouillis 1,391

It's been 2 years since I last updated this collection, but today, I found a new way to compute the roots of a polynomial from python: one can use the python interface to scilab, called sciscipy. Here is the code

``````>>> from scilab import scilab as sci
>>> sci.roots([3.4, -5.22, 2, 1, -7.1])
array([-0.91341205+0.j        ,  1.57667814+0.j        ,
0.43601402+1.12245146j,  0.43601402-1.12245146j])
``````

In ubuntu, install sciscipy with

``````sudo aptitude install python-sciscipy
``````

I did not test sciscipy in windows.

Be a part of the DaniWeb community

We're a friendly, industry-focused community of developers, IT pros, digital marketers, and technology enthusiasts meeting, networking, learning, and sharing knowledge.