# any Python equivalent of Math::Polynomial::Solve?

Discussion in 'Python' started by Just, Feb 26, 2005.

1. ### JustGuest

While googling for a non-linear equation solver, I found
Math:olynomial::Solve in CPAN. It seems a great little module, except
it's not Python... I'm especially looking for its poly_root()
functionality (which solves arbitrary polynomials). Does anyone know of
a Python module/package that implements that?

Just
Just, Feb 26, 2005

2. ### Nick CoghlanGuest

Just wrote:
> While googling for a non-linear equation solver, I found
> Math:olynomial::Solve in CPAN. It seems a great little module, except
> it's not Python... I'm especially looking for its poly_root()
> functionality (which solves arbitrary polynomials). Does anyone know of
> a Python module/package that implements that?
>
> Just

Does SciPy do what you want? Specifically Scientific.Functions.FindRoot [1] &
Scientific.Functions.Polynomial [2]

Regards,
Nick.

[1]
http://starship.python.net/~hinsen/ScientificPython/ScientificPythonManual/Scientific_9.html
[2]
http://starship.python.net/~hinsen/ScientificPython/ScientificPythonManual/Scientific_13.html

--
Nick Coghlan | | Brisbane, Australia
---------------------------------------------------------------
http://boredomandlaziness.skystorm.net
Nick Coghlan, Feb 26, 2005

3. ### JustGuest

In article <>,
Nick Coghlan <> wrote:

> Just wrote:
> > While googling for a non-linear equation solver, I found
> > Math:olynomial::Solve in CPAN. It seems a great little module, except
> > it's not Python... I'm especially looking for its poly_root()
> > functionality (which solves arbitrary polynomials). Does anyone know of
> > a Python module/package that implements that?
> >
> > Just

>
> Does SciPy do what you want? Specifically Scientific.Functions.FindRoot [1] &
> Scientific.Functions.Polynomial [2]
>
> Regards,
> Nick.
>
> [1]
> http://starship.python.net/~hinsen/ScientificPython/ScientificPythonManual/Sci
> entific_9.html
> [2]
> http://starship.python.net/~hinsen/ScientificPython/ScientificPythonManual/Sci
> entific_13.html

module.)

I had played with [1], but it "only" calculates one root, and I need all
roots (specifically, for a quintic equation). [2] doesn't seem to be a
solver?

Just
Just, Feb 26, 2005
4. ### Terry ReedyGuest

"Just" <> wrote in message
news:4all.nl...
>> Does SciPy do what you want? Specifically Scientific.Functions.FindRoot
>> [1] &
>> Scientific.Functions.Polynomial [2]
>> http://starship.python.net/~hinsen/ScientificPython/ScientificPythonManual/Sci
>> entific_9.html
>> [2]
>> http://starship.python.net/~hinsen/ScientificPython/ScientificPythonManual/Sci
>> entific_13.html

>
> (Hm, I had the impression that scipy != Konrad Hinsen's Scientific
> module.)

www.scipy.org (first hit on "Python SciPy" google search)

Terry J. Reedy
Terry Reedy, Feb 26, 2005
5. ### John M. GambleGuest

In article <4all.nl>,
Just <> wrote:
>While googling for a non-linear equation solver, I found
>Math:olynomial::Solve in CPAN. It seems a great little module, except

Thank you.

>it's not Python...

> I'm especially looking for its poly_root()
>functionality (which solves arbitrary polynomials). Does anyone know of
>a Python module/package that implements that?

Are you looking for that particular algorithm, or for any
source that will find the roots of the polynomial? The
original source for the algorithm used in the module is
from Hiroshi Murakami's Fortran source, and it shouldn't
be too difficult to repeat the translation process to python.

--
-john

February 28 1997: Last day libraries could order catalogue cards
from the Library of Congress.
John M. Gamble, Feb 26, 2005
6. ### JustGuest

In article <cvqgn5$n0b$>,
(John M. Gamble) wrote:

> In article <4all.nl>,
> Just <> wrote:
> >While googling for a non-linear equation solver, I found
> >Math:olynomial::Solve in CPAN. It seems a great little module, except

>
> Thank you.
>
> >it's not Python...

>

Heh, how big are the odds you find the author of an arbitrary Perl
module on c.l.py...

> > I'm especially looking for its poly_root()
> >functionality (which solves arbitrary polynomials). Does anyone know of
> >a Python module/package that implements that?

>
> Are you looking for that particular algorithm, or for any
> source that will find the roots of the polynomial?

Any will do. As I wrote in another post, I'm currently only looking for
a quintic equation solver, which your module does very nicely.

> The
> original source for the algorithm used in the module is
> from Hiroshi Murakami's Fortran source, and it shouldn't
> be too difficult to repeat the translation process to python.

Ah ok, I'll try to locate that (following the instruction in Solve.pm
didn't work for me ).

Just
Just, Feb 26, 2005
7. ### John M. GambleGuest

In article <4all.nl>,
Just <> wrote:
>
>Heh, how big are the odds you find the author of an arbitrary Perl
>module on c.l.py...
>

Hey, that's why it's called lurking.

>
>Any will do. As I wrote in another post, I'm currently only looking for
>a quintic equation solver, which your module does very nicely.
>
>> The
>> original source for the algorithm used in the module is
>> from Hiroshi Murakami's Fortran source, and it shouldn't
>> be too difficult to repeat the translation process to python.

>
>Ah ok, I'll try to locate that (following the instruction in Solve.pm
>didn't work for me ).
>

Ouch. I just did a quick search and found that that site has undergone
a few changes, and the code that i reference is missing. A few other
links in the docs are stale too. I need to update the documentation.

Anyway, doing a search for 'hqr' and Eispack got me a lot of sites.
In particular, this one is pretty friendly:

<http://netlib.enseeiht.fr/eispack/>

Look at the source for balanc.f (does the prep-work) and hqr.f
(does the solving). Minor annoyance: the real and imaginary
parts of the roots are in separate arrays. I combined them into
complex types in my perl source, in case you want to make a
comparison.

Of course, all this may be moot if the other suggestions
work out.

--
-john

February 28 1997: Last day libraries could order catalogue cards
from the Library of Congress.
John M. Gamble, Feb 26, 2005
8. ### JustGuest

In article <cvqmvt$2rm$>,
(John M. Gamble) wrote:

> >> The
> >> original source for the algorithm used in the module is
> >> from Hiroshi Murakami's Fortran source, and it shouldn't
> >> be too difficult to repeat the translation process to python.

> >
> >Ah ok, I'll try to locate that (following the instruction in Solve.pm
> >didn't work for me ).
> >

>
> Ouch. I just did a quick search and found that that site has undergone
> a few changes, and the code that i reference is missing. A few other
> links in the docs are stale too. I need to update the documentation.
>
> Anyway, doing a search for 'hqr' and Eispack got me a lot of sites.
> In particular, this one is pretty friendly:
>
> <http://netlib.enseeiht.fr/eispack/>
>
> Look at the source for balanc.f (does the prep-work) and hqr.f
> (does the solving). Minor annoyance: the real and imaginary
> parts of the roots are in separate arrays. I combined them into
> complex types in my perl source, in case you want to make a
> comparison.

Thanks! I'll check that out.

> Of course, all this may be moot if the other suggestions
> work out.

SciPy indeed appear to contain a solver, but I'm currently stuck in
trying to _get_ it for my platform (OSX). I'm definitely not going to
install a Fortran compiler just to evaluate it (even though my name is
not "Ilias" ;-). Also, SciPy is _huge_, so maybe a Python translation of
that Fortran code or your Perl code will turn out to be more attractive
after all...

Just
Just, Feb 26, 2005
9. ### Raymond L. BuvelGuest

Just wrote:
<snip>
>
> SciPy indeed appear to contain a solver, but I'm currently stuck in
> trying to _get_ it for my platform (OSX). I'm definitely not going to
> install a Fortran compiler just to evaluate it (even though my name is
> not "Ilias" ;-). Also, SciPy is _huge_, so maybe a Python translation of
> that Fortran code or your Perl code will turn out to be more attractive
> after all...
>
> Just

The GNU Scientific Library has a nice root finder for polynomials with
real coefficients. I have wrapped this with Pyrex to work with my
ratfun module see:

http://calcrpnpy.sourceforge.net/ratfun.html

If this will suit your needs, I can send you an alpha release of the
package with the root finder. It is not pure Python. I requires Pyrex
and a C compiler to install. My guess is that it will work on OSX as
well as it does on Linux. This functionality will be included in the
next release of the ratfun package but I still have to unit test a
number of components and update the documentation. Consequently, an
official release will not happen soon.

Ray
Raymond L. Buvel, Feb 27, 2005
10. ### JustGuest

In article <I68Ud.577\$-kc.rr.com>,
"Raymond L. Buvel" <> wrote:

> Just wrote:
> <snip>
> >
> > SciPy indeed appear to contain a solver, but I'm currently stuck in
> > trying to _get_ it for my platform (OSX). I'm definitely not going to
> > install a Fortran compiler just to evaluate it (even though my name is
> > not "Ilias" ;-). Also, SciPy is _huge_, so maybe a Python translation of
> > that Fortran code or your Perl code will turn out to be more attractive
> > after all...
> >
> > Just

>
> The GNU Scientific Library has a nice root finder for polynomials with
> real coefficients. I have wrapped this with Pyrex to work with my
> ratfun module see:
>
> http://calcrpnpy.sourceforge.net/ratfun.html
>
> If this will suit your needs, I can send you an alpha release of the
> package with the root finder. It is not pure Python. I requires Pyrex
> and a C compiler to install. My guess is that it will work on OSX as
> well as it does on Linux. This functionality will be included in the
> next release of the ratfun package but I still have to unit test a
> number of components and update the documentation. Consequently, an
> official release will not happen soon.

Thank you, I'll check this out. I had come across GSL, but not Python
bindings. (GPL is probably a problem for my project, but it's very good
to know anyway.)

On the other hand, I just finished translating the relevant portions of
Math:olynomial::Solve to Python, so I'm probably all set, at least for
now. Thanks everyone for the responses, especially to John Gamble!

Just
Just, Feb 27, 2005
11. ### Nick CoghlanGuest

Just wrote:
> (Hm, I had the impression that scipy != Konrad Hinsen's Scientific
> module.)

You're probably right

> I had played with [1], but it "only" calculates one root, and I need all
> roots (specifically, for a quintic equation). [2] doesn't seem to be a
> solver?

Actually, I was curious whether the 'zeros' method in [2] did the right thing.

Cheers,
Nick.

--
Nick Coghlan | | Brisbane, Australia
---------------------------------------------------------------
http://boredomandlaziness.skystorm.net
Nick Coghlan, Feb 27, 2005
12. ### Carl BanksGuest

Just wrote:
> While googling for a non-linear equation solver, I found
> Math:olynomial::Solve in CPAN. It seems a great little module,

except
> it's not Python... I'm especially looking for its poly_root()
> functionality (which solves arbitrary polynomials). Does anyone know

of
> a Python module/package that implements that?

If you don't have a great need for speed, you can accomplish this
easily with the linear algebra module of Numeric/numarray. Suppose
your quintic polynomial's in the form

a + b*x + c*x**2 + d*x**3 + e*x**4 + x**5

The roots of it are equal to the eigenvalues of the companion matrix:

0 1 0 0 0
0 0 1 0 0
0 0 0 1 0
0 0 0 0 1
-a -b -c -d -e

It should be pretty easy to set up a Numeric matrix and call
LinearAlgebra.eigenvalues. For example, here is a simple quintic
solver:

.. from Numeric import *
.. from LinearAlgebra import *
..
.. def quinticroots(p):
.. cm = zeros((5,5),Float32)
.. cm[0,1] = cm[1,2] = cm[2,3] = cm[3,4] = 1.0
.. cm[4,0] = -p[0]
.. cm[4,1] = -p[1]
.. cm[4,2] = -p[2]
.. cm[4,3] = -p[3]
.. cm[4,4] = -p[4]
.. return eigenvalues(cm)

now-you-can-find-all-five-Lagrange-points-ly yr's,

--
CARL BANKS
Carl Banks, Feb 27, 2005
13. ### JustGuest

In article <>,
"Carl Banks" <> wrote:

> Just wrote:
> > While googling for a non-linear equation solver, I found
> > Math:olynomial::Solve in CPAN. It seems a great little module,

> except
> > it's not Python... I'm especially looking for its poly_root()
> > functionality (which solves arbitrary polynomials). Does anyone know

> of
> > a Python module/package that implements that?

>
> If you don't have a great need for speed, you can accomplish this
> easily with the linear algebra module of Numeric/numarray. Suppose
> your quintic polynomial's in the form
>
> a + b*x + c*x**2 + d*x**3 + e*x**4 + x**5
>
> The roots of it are equal to the eigenvalues of the companion matrix:
>
> 0 1 0 0 0
> 0 0 1 0 0
> 0 0 0 1 0
> 0 0 0 0 1
> -a -b -c -d -e
>
> It should be pretty easy to set up a Numeric matrix and call
> LinearAlgebra.eigenvalues. For example, here is a simple quintic
> solver:
>
> . from Numeric import *
> . from LinearAlgebra import *
> .
> . def quinticroots(p):
> . cm = zeros((5,5),Float32)
> . cm[0,1] = cm[1,2] = cm[2,3] = cm[3,4] = 1.0
> . cm[4,0] = -p[0]
> . cm[4,1] = -p[1]
> . cm[4,2] = -p[2]
> . cm[4,3] = -p[3]
> . cm[4,4] = -p[4]
> . return eigenvalues(cm)
>
>
> now-you-can-find-all-five-Lagrange-points-ly yr's,

Wow, THANKS. This was the answer I was secretly hoping for... "Great
need for speed", no, not really, but this Numeric-based version is about
9 times faster than what I translated from Perl code yesterday, so from
where I'm standing your version is blazingly fast...

Thanks again,

Just
Just, Feb 27, 2005
14. ### Alex ReneltGuest

Just wrote:
> In article <>,
> "Carl Banks" <> wrote:

>>It should be pretty easy to set up a Numeric matrix and call
>>LinearAlgebra.eigenvalues. For example, here is a simple quintic
>>solver:
>>
>>. from Numeric import *
>>. from LinearAlgebra import *
>>.
>>. def quinticroots(p):
>>. cm = zeros((5,5),Float32)
>>. cm[0,1] = cm[1,2] = cm[2,3] = cm[3,4] = 1.0
>>. cm[4,0] = -p[0]
>>. cm[4,1] = -p[1]
>>. cm[4,2] = -p[2]
>>. cm[4,3] = -p[3]
>>. cm[4,4] = -p[4]
>>. return eigenvalues(cm)
>>
>>
>>now-you-can-find-all-five-Lagrange-points-ly yr's,

>
>
> Wow, THANKS. This was the answer I was secretly hoping for... "Great
> need for speed", no, not really, but this Numeric-based version is about
> 9 times faster than what I translated from Perl code yesterday, so from
> where I'm standing your version is blazingly fast...
>
> Thanks again,
>
> Just

I'm writing a class for polynomial manipulation. The generalization of
the above code is:

definitions:
1.) p = array([a_0, a_i, ..., a_n]) represents your polynomial
P(x) = \sum _{i=0} ^n a_i x^i

2.) deg(p) is its degree

3.) monic(p) makes P monic, i.e. monic(p) = p / p[:-1]

then you get:
from numarray import *
import numarray.linear_algebra as la

def roots(p):
p = monic(p); n = deg(p)
M = asarray(zeros((n,n)), typecode = 'f8')
# or 'c16' if you need complex coefficients
M[:-1,1:] = identity(n-1)
M[-1,:] = -p[:-1]
return la.eigenvalues(M)

Alex
Alex Renelt, Feb 27, 2005
15. ### Alex ReneltGuest

Alex Renelt wrote:

> I'm writing a class for polynomial manipulation. The generalization of
> the above code is:
>
> definitions:
> 1.) p = array([a_0, a_i, ..., a_n]) represents your polynomial
> P(x) = \sum _{i=0} ^n a_i x^i
>
> 2.) deg(p) is its degree
>
> 3.) monic(p) makes P monic, i.e. monic(p) = p / p[:-1]
>
> then you get:
> from numarray import *
> import numarray.linear_algebra as la
>
> def roots(p):
> p = monic(p); n = deg(p)
> M = asarray(zeros((n,n)), typecode = 'f8')
> # or 'c16' if you need complex coefficients
> M[:-1,1:] = identity(n-1)
> M[-1,:] = -p[:-1]
> return la.eigenvalues(M)
>
> Alex

under definitions, 3.)
its "monic(p) = p / p[-1]" of course

Alex
Alex Renelt, Feb 27, 2005
16. ### Raymond L. BuvelGuest

Alex Renelt wrote:
> Alex Renelt wrote:
>
>> I'm writing a class for polynomial manipulation. The generalization of
>> the above code is:
>>
>> definitions:
>> 1.) p = array([a_0, a_i, ..., a_n]) represents your polynomial
>> P(x) = \sum _{i=0} ^n a_i x^i
>>
>> 2.) deg(p) is its degree
>>
>> 3.) monic(p) makes P monic, i.e. monic(p) = p / p[:-1]
>>
>> then you get:
>> from numarray import *
>> import numarray.linear_algebra as la
>>
>> def roots(p):
>> p = monic(p); n = deg(p)
>> M = asarray(zeros((n,n)), typecode = 'f8')
>> # or 'c16' if you need complex coefficients
>> M[:-1,1:] = identity(n-1)
>> M[-1,:] = -p[:-1]
>> return la.eigenvalues(M)
>>
>> Alex

>
>
> uhh, I made a mistake:
> under definitions, 3.)
> its "monic(p) = p / p[-1]" of course
>
> Alex

Alex,

If you want a class for polynomial manipulation, you should check out my
ratfun module.

http://calcrpnpy.sourceforge.net/ratfun.html

Ray
Raymond L. Buvel, Feb 27, 2005
17. ### David E. Konerding DSD staffGuest

On 2005-02-26, Just <> wrote:
> While googling for a non-linear equation solver, I found
> Math:olynomial::Solve in CPAN. It seems a great little module, except
> it's not Python... I'm especially looking for its poly_root()
> functionality (which solves arbitrary polynomials). Does anyone know of
> a Python module/package that implements that?
>
> Just

Although I dont' really work on it any more, the Py-ML module which interfaces
Python to Mathematica would almost certain be able to do this, although you'd need an installation of
Mathematica.

http://sourceforge.net/projects/py-ml/
David E. Konerding DSD staff, Feb 27, 2005
18. ### Carl BanksGuest

Carl Banks wrote:
> . from Numeric import *
> . from LinearAlgebra import *
> .
> . def quinticroots(p):
> . cm = zeros((5,5),Float32)
> . cm[0,1] = cm[1,2] = cm[2,3] = cm[3,4] = 1.0
> . cm[4,0] = -p[0]
> . cm[4,1] = -p[1]
> . cm[4,2] = -p[2]
> . cm[4,3] = -p[3]
> . cm[4,4] = -p[4]
> . return eigenvalues(cm)

Here's an improved version. It uses 64-bit numbers (I had used type
Float32 because I often use a float32 type at work, not in Python,
unfortunately), and array assignment.

.. def quinticroots(p):
.. cm = zeros((5,5),Float)
.. cm[0,1] = cm[1,2] = cm[2,3] = cm[3,4] = 1.0
.. cm[4,:] = -array(p)
.. return eigenvalues(cm)

--
CARL BANKS
Carl Banks, Feb 27, 2005
19. ### Cousin StanleyGuest

Alex ....

Thanks for posting your generalized numarray
eigenvalue solution ....

It's been almost 30 years since I've looked at
any characteristic equation, eigenvalue, eignevector
type of processing and at this point I don't recall
many of the particulars ....

Not being sure about the nature of the monic( p ) function,
I implemented it as an element-wise division of each
of the coefficients ....

Is this anywhere near the correct interpretation
for monic( p ) ?

Using the version below, Python complained

.. M[ -1 , : ] = -p[ : -1 ]

I tried without the slicing on the right ....

.. . M[ -1 , : ] = -p[ -1 ]

The following code will run and produce results,
but I'm wondering if I've totally screwed it up
since the results I obtain are different from
those obtained from the specific 5th order Numeric
solution previously posted here ....

.. from numarray import *
..
.. import numarray.linear_algebra as LA
..
.. def monic( this_list ) :
..
.. m = [ ]
..
.. last_item = this_list[ -1 ]
..
.. for this_item in this_list :
..
.. m.append( this_item / last_item )
..
.. return m
..
..
.. def roots( p ) :
..
.. p = monic( p )
..
.. n = len( p ) # degree of polynomial
..
.. z = zeros( ( n , n ) )
..
.. M = asarray( z , typecode = 'f8' ) # typecode = c16, complex
..
.. M[ : -1 , 1 : ] = identity( n - 1 )
..
.. M[ -1 , : ] = -p[ -1 ] # removed : slicing on the right
..
.. return LA.eigenvalues( M )
..
..
.. coeff = [ 1. , 3. , 5. , 7. , 9. ]
..
.. print ' Coefficients ..'
.. print
.. print ' %s' % coeff
.. print
.. print ' Eigen Values .. '
.. print
..
.. eigen_values = roots( coeff )
..
.. for this_value in eigen_values :
..
.. print ' %s' % this_value
..

Any clues would be greatly appreciated ....

--
Stanley C. Kitching
Human Being
Phoenix, Arizona
Cousin Stanley, Mar 1, 2005
20. ### Pierre SchnizerGuest

In case you are still interested pygsl wraps the GSL solver.
<snip>
from pygsl import poly
pc = poly.poly_complex(3)
tmp, rs = pc.solve((2,3,1))
print rs
</snip>

You get pygsl at http://sourceforge.net/projects/pygsl/
Pierre
Pierre Schnizer, Mar 2, 2005