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

[Cousin Stanley]

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

 

[konrad.hinsen]

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.

 

[Alex Renelt]

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

 

[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

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

 

[Cousin Stanley]

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

 

[Cousin Stanley]

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