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

Discussion in 'Python' started by Just, Feb 26, 2005.

  1. Just

    Just Guest

    While googling for a non-linear 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
    Just, Feb 26, 2005
    #1
    1. Advertising

  2. Just

    Nick Coghlan Guest

    Just wrote:
    > While googling for a non-linear 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/ScientificPython/ScientificPythonManual/Scientific_9.html
    [2]
    http://starship.python.net/~hinsen/ScientificPython/ScientificPythonManual/Scientific_13.html


    --
    Nick Coghlan | | Brisbane, Australia
    ---------------------------------------------------------------
    http://boredomandlaziness.skystorm.net
    Nick Coghlan, Feb 26, 2005
    #2
    1. Advertising

  3. Just

    Just Guest

    In article <>,
    Nick Coghlan <> wrote:

    > Just wrote:
    > > While googling for a non-linear 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/ScientificPython/ScientificPythonManual/Sci
    > entific_9.html
    > [2]
    > http://starship.python.net/~hinsen/ScientificPython/ScientificPythonManual/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
    Just, Feb 26, 2005
    #3
  4. Just

    Terry Reedy Guest

    "Just" <> wrote in message
    news:4all.nl...
    >> Does SciPy do what you want? Specifically Scientific.Functions.FindRoot
    >> [1] &
    >> Scientific.Functions.Polynomial [2]
    >> http://starship.python.net/~hinsen/ScientificPython/ScientificPythonManual/Sci
    >> entific_9.html
    >> [2]
    >> http://starship.python.net/~hinsen/ScientificPython/ScientificPythonManual/Sci
    >> entific_13.html

    >
    > (Hm, I had the impression that scipy != Konrad Hinsen's Scientific
    > module.)


    www.scipy.org (first hit on "Python SciPy" google search)

    Terry J. Reedy
    Terry Reedy, Feb 26, 2005
    #4
  5. In article <4all.nl>,
    Just <> wrote:
    >While googling for a non-linear 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.
    John M. Gamble, Feb 26, 2005
    #5
  6. Just

    Just Guest

    In article <cvqgn5$n0b$>,
    (John M. Gamble) wrote:

    > In article <4all.nl>,
    > Just <> wrote:
    > >While googling for a non-linear 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
    Just, Feb 26, 2005
    #6
  7. In article <4all.nl>,
    Just <> 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 prep-work) 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.
    John M. Gamble, Feb 26, 2005
    #7
  8. Just

    Just Guest

    In article <cvqmvt$2rm$>,
    (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 prep-work) 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, Feb 26, 2005
    #8
  9. 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
    Raymond L. Buvel, Feb 27, 2005
    #9
  10. Just

    Just Guest

    In article <I68Ud.577$-kc.rr.com>,
    "Raymond L. Buvel" <> 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
    Just, Feb 27, 2005
    #10
  11. Just

    Nick Coghlan Guest

    Just wrote:
    > (Hm, I had the impression that scipy != Konrad Hinsen's Scientific
    > module.)


    You're probably right :)

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


    Actually, I was curious whether the 'zeros' method in [2] did the right thing.

    Cheers,
    Nick.

    --
    Nick Coghlan | | Brisbane, Australia
    ---------------------------------------------------------------
    http://boredomandlaziness.skystorm.net
    Nick Coghlan, Feb 27, 2005
    #11
  12. Just

    Carl Banks Guest

    Just wrote:
    > While googling for a non-linear 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?


    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

    It should be pretty easy to set up a Numeric matrix and call
    LinearAlgebra.eigenvalues. For example, here is a simple quintic
    solver:

    .. from Numeric import *
    .. from LinearAlgebra import *
    ..
    .. def quinticroots(p):
    .. cm = zeros((5,5),Float32)
    .. cm[0,1] = cm[1,2] = cm[2,3] = cm[3,4] = 1.0
    .. cm[4,0] = -p[0]
    .. cm[4,1] = -p[1]
    .. cm[4,2] = -p[2]
    .. cm[4,3] = -p[3]
    .. cm[4,4] = -p[4]
    .. return eigenvalues(cm)


    now-you-can-find-all-five-Lagrange-points-ly yr's,

    --
    CARL BANKS
    Carl Banks, Feb 27, 2005
    #12
  13. Just

    Just Guest

    In article <>,
    "Carl Banks" <> wrote:

    > Just wrote:
    > > While googling for a non-linear 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?

    >
    > 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
    >
    > It should be pretty easy to set up a Numeric matrix and call
    > LinearAlgebra.eigenvalues. For example, here is a simple quintic
    > solver:
    >
    > . from Numeric import *
    > . from LinearAlgebra import *
    > .
    > . def quinticroots(p):
    > . cm = zeros((5,5),Float32)
    > . cm[0,1] = cm[1,2] = cm[2,3] = cm[3,4] = 1.0
    > . cm[4,0] = -p[0]
    > . cm[4,1] = -p[1]
    > . cm[4,2] = -p[2]
    > . cm[4,3] = -p[3]
    > . cm[4,4] = -p[4]
    > . return eigenvalues(cm)
    >
    >
    > now-you-can-find-all-five-Lagrange-points-ly yr's,


    Wow, THANKS. This was the answer I was secretly hoping for... "Great
    need for speed", no, not really, but this Numeric-based version is about
    9 times faster than what I translated from Perl code yesterday, so from
    where I'm standing your version is blazingly fast...

    Thanks again,

    Just
    Just, Feb 27, 2005
    #13
  14. Just

    Alex Renelt Guest

    Just wrote:
    > In article <>,
    > "Carl Banks" <> wrote:


    >>It should be pretty easy to set up a Numeric matrix and call
    >>LinearAlgebra.eigenvalues. For example, here is a simple quintic
    >>solver:
    >>
    >>. from Numeric import *
    >>. from LinearAlgebra import *
    >>.
    >>. def quinticroots(p):
    >>. cm = zeros((5,5),Float32)
    >>. cm[0,1] = cm[1,2] = cm[2,3] = cm[3,4] = 1.0
    >>. cm[4,0] = -p[0]
    >>. cm[4,1] = -p[1]
    >>. cm[4,2] = -p[2]
    >>. cm[4,3] = -p[3]
    >>. cm[4,4] = -p[4]
    >>. return eigenvalues(cm)
    >>
    >>
    >>now-you-can-find-all-five-Lagrange-points-ly yr's,

    >
    >
    > Wow, THANKS. This was the answer I was secretly hoping for... "Great
    > need for speed", no, not really, but this Numeric-based version is about
    > 9 times faster than what I translated from Perl code yesterday, so from
    > where I'm standing your version is blazingly fast...
    >
    > Thanks again,
    >
    > Just


    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
    Alex Renelt, Feb 27, 2005
    #14
  15. Just

    Alex Renelt Guest

    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(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 Renelt, Feb 27, 2005
    #15
  16. 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(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
    Raymond L. Buvel, Feb 27, 2005
    #16
  17. On 2005-02-26, Just <> wrote:
    > While googling for a non-linear 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


    Although I dont' really work on it any more, the Py-ML module which interfaces
    Python to Mathematica would almost certain be able to do this, although you'd need an installation of
    Mathematica.

    http://sourceforge.net/projects/py-ml/
    David E. Konerding DSD staff, Feb 27, 2005
    #17
  18. Just

    Carl Banks Guest

    Carl Banks wrote:
    > . from Numeric import *
    > . from LinearAlgebra import *
    > .
    > . def quinticroots(p):
    > . cm = zeros((5,5),Float32)
    > . cm[0,1] = cm[1,2] = cm[2,3] = cm[3,4] = 1.0
    > . cm[4,0] = -p[0]
    > . cm[4,1] = -p[1]
    > . cm[4,2] = -p[2]
    > . cm[4,3] = -p[3]
    > . cm[4,4] = -p[4]
    > . return eigenvalues(cm)



    Here's an improved version. It uses 64-bit numbers (I had used type
    Float32 because I often use a float32 type at work, not in Python,
    unfortunately), and array assignment.

    .. def quinticroots(p):
    .. cm = zeros((5,5),Float)
    .. cm[0,1] = cm[1,2] = cm[2,3] = cm[3,4] = 1.0
    .. cm[4,:] = -array(p)
    .. return eigenvalues(cm)


    --
    CARL BANKS
    Carl Banks, Feb 27, 2005
    #18
  19. 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 ....

    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 follow-up,
    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
    Cousin Stanley, Mar 1, 2005
    #19
  20. 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
    Pierre Schnizer, Mar 2, 2005
    #20
    1. Advertising

Want to reply to this thread or ask your own question?

It takes just 2 minutes to sign up (and it's free!). Just click the sign up button to choose a username and then you can ask your own questions on the forum.
Similar Threads
  1. chirs
    Replies:
    18
    Views:
    743
    Chris Uppal
    Mar 2, 2004
  2. Chen L.
    Replies:
    2
    Views:
    861
    Dik T. Winter
    Jul 6, 2004
  3. jitasi
    Replies:
    1
    Views:
    732
    Terry Reedy
    Mar 4, 2007
  4. Kelly Jones
    Replies:
    4
    Views:
    162
    Kelly Jones
    Apr 1, 2010
  5. VK
    Replies:
    15
    Views:
    1,111
    Dr J R Stockton
    May 2, 2010
Loading...

Share This Page