| ....
| 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
--
http://mail.python.org/mailman/listinfo/python-list
