468,168 Members | 1,461 Online
Bytes | Developer Community
New Post

Home Posts Topics Members FAQ

Post your question to a community of 468,168 developers. It's quick & easy.

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

While googling for a non-linear equation solver, I found
Math::Polynomial::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
Jul 18 '05 #1
17 3230
Just wrote:
While googling for a non-linear equation solver, I found
Math::Polynomial::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/S...entific_9.html
[2]
http://starship.python.net/~hinsen/S...ntific_13.html
--
Nick Coghlan | nc******@email.com | Brisbane, Australia
---------------------------------------------------------------
http://boredomandlaziness.skystorm.net
Jul 18 '05 #2
In article <ma***************************************@python. org>,
Nick Coghlan <nc******@iinet.net.au> wrote:
Just wrote:
While googling for a non-linear equation solver, I found
Math::Polynomial::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/S...thonManual/Sci
entific_9.html
[2]
http://starship.python.net/~hinsen/S...thonManual/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
Jul 18 '05 #3

"Just" <ju**@xs4all.nl> wrote in message
news:ju************************@news1.news.xs4all. nl...
Does SciPy do what you want? Specifically Scientific.Functions.FindRoot
[1] &
Scientific.Functions.Polynomial [2]
http://starship.python.net/~hinsen/S...thonManual/Sci
entific_9.html
[2]
http://starship.python.net/~hinsen/S...thonManual/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

Jul 18 '05 #4
In article <ju************************@news1.news.xs4all.nl >,
Just <ju**@xs4all.nl> wrote:
While googling for a non-linear equation solver, I found
Math::Polynomial::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.

--
-john

February 28 1997: Last day libraries could order catalogue cards
from the Library of Congress.
Jul 18 '05 #5
In article <cv**********@e250.ripco.com>,
jg*****@ripco.com (John M. Gamble) wrote:
In article <ju************************@news1.news.xs4all.nl >,
Just <ju**@xs4all.nl> wrote:
While googling for a non-linear equation solver, I found
Math::Polynomial::Solve in CPAN. It seems a great little module, except
Thank you.
it's not Python...


Sorry about that.


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
Jul 18 '05 #6
In article <ju************************@news1.news.xs4all.nl >,
Just <ju**@xs4all.nl> 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.
Jul 18 '05 #7
In article <cv**********@e250.ripco.com>,
jg*****@ripco.com (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
Jul 18 '05 #8
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
Jul 18 '05 #9
In article <I6***************@twister.rdc-kc.rr.com>,
"Raymond L. Buvel" <le******@wi.rr.com> 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::Polynomial::Solve to Python, so I'm probably all set, at least for
now. Thanks everyone for the responses, especially to John Gamble!

Just
Jul 18 '05 #10
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 ....
--
Stanley C. Kitching
Human Being
Phoenix, Arizona
Jul 18 '05 #11

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
Jul 18 '05 #12
| In case you are still interested pygsl wraps the GSL solver.
| ....
| from pygsl import poly
|
| pc = poly.poly_complex( 3 )
|
| tmp, rs = pc.solve( ( 2 , 3 , 1 ) )
|
| print rs
| ....
|
| You get pygsl at http://sourceforge.net/projects/pygsl/

Pierre ....

I am still interested and have downloaded the PyGSL source
from SourceForge and will attempt to build under Debian Linux ....

Thanks for the information ....
--
Stanley C. Kitching
Human Being
Phoenix, Arizona
Jul 18 '05 #13
Carl Banks wrote:
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


The method "zeros()" in Scientific.Functions.Polynomial uses exactly
that trick for finding the zeros of a general polynomial. If you need
to do more with polynomials than just finding the zeros, the Polynomial
class is probably better than an on-the-spot solution.

Root finding through eigenvalues is not the fastest method, but it's
simple and stable, and not terribly bad either.

Sorry for not making that comment earlier, I don't have the time to
follow this list at the moment (to my great regret), but I was made
aware of this thread through PythonURL.

Konrad.

Jul 18 '05 #14
Cousin Stanley wrote:
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 ....
This is correct. The aim is that p has a leading coefficient that equals 1.

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

Using the version below, Python complained
about the line ....

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

Works for me. Did you perhaps use a list (type(p) == type([])) for p?
Then python does not know what -p means (numeric or numarray does).
So, in view of you comments about slicing in you follow-up,
I tried without the slicing on the right ....

. . M[ -1 , : ] = -p[ -1 ]
That's wrong because you don't set a slice but a single item!
Old code should work. We need all coefficients but the leading coeff. So
take the slice 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
your function equals the following one:

def monic(p): return p / p[-1]

But you have to ensure that p is an array object. It does element-wise
operations per default.

Remember we need future division or take float(p[-1]) or the denominator. .
.
. def roots( p ) :
.
. p = monic( p )
.
. n = len( p ) # degree of polynomial The degree is len(p) -1 or something smaller (some people are calling
len(p) -1 a "degree bound" instead). It is smaller if p contains leading
zeros which should be deleted, e.g. P(x) = x^2 + 4 x + 4 could be entered as
p = array([4, 4, 1, 0, 0])
which would produce a deg of 4 instead of 2. .
. 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. ] remember to use an array or convert it on-the-fly inside your roots
function:

M[-1,:] = - asarray(p)[:-1] .
. 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 ....

Hope that helps.

Alex
Jul 18 '05 #15
Raymond L. Buvel wrote:
Alex Renelt wrote:
Alex Renelt wrote:
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


Ray,

thanks a lot for your hint but I'm writing it for a students paper in a
german math class so I believe I should better do some work alone ;-)

In addition I only need a class for polynomials and not for rational
functions and I'm testing different iterative polynomial solvers. So I'm
happy to have my own small class which I understand 100%.

Generally I'm against rediscovering the wheel again and again!

Alex
Jul 18 '05 #16
| In case you are still interested pygsl wraps the GSL solver.
| ....
|
| from pygsl import poly
|
| pc = poly.poly_complex( 3 )
|
| tmp , rs = pc.solve( ( 2 , 3 , 1 ) )
|
| print rs
| ....
|
| You get pygsl at http://sourceforge.net/projects/pygsl/

Pierre ....

Thanks again for the link to the Python GNU Scientific Library ....

I managed to build it under Debian GNU/Linux
and the solution you posted above works as is ....

The PyGSL results agree with the numarray eigenvalue
solution posted by Alex Renelt and a few timing tests
for small problems seem to indicate that it runs a bit quicker ....

I'm looking forward to further use of this package
in the future ....
--
Stanley C. Kitching
Human Being
Phoenix, Arizona
Jul 18 '05 #17
| ....
| Did you perhaps use a list (type(p) == type([])) for p?
| ....

Alex ....

Using the coefficients in an array instead of a list
was the key in the solution to my problems ....

Your other suggestions regarding floating p
and the off-by-one error that I had with the
polynomial degree were also included ....

The results agree with solutions from PyGSL
as suggested by Pierre Schnizer, but seem
to run just a bit slower on my machine ....

Thanks again for your assistance ....

The version that I tested with follows ....

Stanley C. Kitching
# -----------------------------------------------------------------

#!/usr/bin/env python

'''
NewsGroup .... comp.lang.python
Date ......... 2005-02-27
Subject ...... any Python equivalent of Math::Polynomial::Solver
Reply_By ..... Alex Renelt
Edited_By .... Stanley c. Kitching

I'm writing a class for polynomial manipulation.

The generalization of the above code
for providing eigenvalues of a polynomial
is ....

definitions ....

1.) p = array( [ a_0 , a_i , .... , a_n ] )

P( x ) = \sum _{ i = 0 } ^n a_i x^i

2.) deg( p ) is its degree

3.) monic( p ) makes P monic

monic( p ) = p / p[ -1 ]

'''

import sys
import time

from numarray import *

import numarray.linear_algebra as LA
print '\n ' , sys.argv[ 0 ] , '\n'

def report( n , this_data ) :

print ' Coefficients ......' , list_coeff[ n ]
print
print ' Roots .........'
print

dt , roots = this_data

for this_item in roots :

print ' %s' % this_item

print
print ' Elapsed Time .... %.6f Seconds' % dt
print
print
def roots( p ) :

p = p / float( p[ -1 ] ) # monic( p )

n = len( p ) - 1 # 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 ]

return LA.eigenvalues( M )
list_coeff = [ array( ( 2. , 3. , 1. ) ) ,
array( ( 1. , 3. , 5. , 7. , 9. ) ) ,
array( ( 10. , 8. , 6. , 4. , 2. , 1. , 2. , 4. , 6. ) ) ]

list_roots = [ ]

for these_coeff in list_coeff :

beg = time.time()

rs = roots( these_coeff )

end = time.time()

dt = end - beg

list_roots.append( [ dt , list( rs ) ] )
i = 0

for this_data in list_roots :

report( i , this_data )

i += 1

print

Jul 18 '05 #18

This discussion thread is closed

Replies have been disabled for this discussion.

Similar topics

699 posts views Thread by mike420 | last post: by
21 posts views Thread by Glen Wheeler | last post: by
4 posts views Thread by Kyler Laird | last post: by
9 posts views Thread by Joseph Aldred | last post: by
2 posts views Thread by sekitoleko | last post: by
5 posts views Thread by MathNewbie | last post: by
116 posts views Thread by Dilip | last post: by
270 posts views Thread by Jordan | last post: by
reply views Thread by kamranasdasdas | last post: by
reply views Thread by gcreed | last post: by
By using this site, you agree to our Privacy Policy and Terms of Use.