While googling for a nonlinear 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 17 3317
In article <ma***************************************@python. org>,
Nick Coghlan <nc******@iinet.net.au> wrote: Just wrote: While googling for a nonlinear 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
In article <ju************************@news1.news.xs4all.nl >,
Just <ju**@xs4all.nl> wrote: While googling for a nonlinear 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.
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 nonlinear 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
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 prepwork) 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.
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 prepwork) 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 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
In article <I6***************@twister.rdckc.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
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 elementwise 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 followup,
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
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
 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
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 onthespot 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.
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 elementwise 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 followup, 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 elementwise
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 onthefly 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
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(n1) 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
 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
 ....
 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 offbyone 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 ......... 20050227
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 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
 
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
           