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

[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::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