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

[Alex Renelt]

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