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

[Just]

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

 

[Nick Coghlan]

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

 

[Just]

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

(Hm, I had the impression that scipy != Konrad Hinsen's Scientific
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

 

[Terry Reedy]

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

 

[John M. Gamble]

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...

Sorry about that.

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.

 

[Just]

Thank you.


Sorry about that.

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

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

 

[John M. Gamble]

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.


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.

 

[Just]

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

 

[Raymond L. Buvel]

Just wrote:

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

 

[Just]

Just wrote:


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

 

[Nick Coghlan]

(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.

 

[Carl Banks]

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,

 

[Just]

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

 

[Alex Renelt]

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

in addition:
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]

in addition:
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

 

[Raymond L. Buvel]

in addition:
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

 

[David E. Konerding DSD staff]

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/

 

[Carl Banks]

. 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)

 

[Cousin Stanley]

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
about the line ....

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

So, in view of you comments about slicing in you follow-up,
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 ....

 

[Pierre Schnizer]

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