[SciPy-user] Finding complex roots of complex polynomials in newscipy

Andrew D drewdemento at yahoo.com
Fri Nov 11 11:31:08 CST 2005


Hello,
I'm using an up-to-date build of newcore and newscipy (newcore
1467, newscipy 1431) on Red Hat Enterprise Linux WS release 3 (Taroon
Update 6) running on a Pentium 4.

Running the "roots" command with complex arguments seems to fail,
producing incorrect results even in the simplest cases:

********************************************************************
Python 2.4.2 (#1, Nov  7 2005, 13:09:14)
[GCC 3.2.3 20030502 (Red Hat Linux 3.2.3-53)] on linux2
Type "help", "copyright", "credits" or "license" for more
information.
>>> from scipy import *
Importing io to scipy
Importing fftpack to scipy
Importing special to scipy
Importing cluster to scipy
Importing sparse to scipy
Importing utils to scipy
Importing interpolate to scipy
Importing lib to scipy
Importing integrate to scipy
Importing signal to scipy
Importing optimize to scipy
Importing linalg to scipy
Importing stats to scipy
>>> roots([1+2.j,1])
array([-1.+0.j])

>>> roots([2.j,1])
Traceback (most recent call last):
  File "<stdin>", line 1, in ?
  File ".../lib/python2.4/site-packages/scipy/base/polynomial.py",
line 112, in roots
    roots = _eigvals(A)
  File ".../lib/python2.4/site-packages/scipy/base/polynomial.py",
line 30, in _eigvals
    return eigvals(arg)
  File ".../lib/python2.4/site-packages/scipy/linalg/decomp.py", line
172, in eigvals
    return eig(a,b=b,left=0,right=0,overwrite_a=overwrite_a)
  File ".../lib/python2.4/site-packages/scipy/linalg/decomp.py", line
114, in eig
    a1 = asarray_chkfinite(a)
  File ".../lib/python2.4/site-packages/scipy/base/function_base.py",
line 211, in asarray_chkfinite
    raise ValueError, "array must not contain infs or NaNs"
ValueError: array must not contain infs or NaNs
********************************************************************


SciPy v0.32 get's it right:
>>> roots([1.+2.j,1])
array([-0.2+0.4j])

>>> roots([2.j,1])
array([-0.+0.5j])


Is this a known problem? The culprit seems to be the lines 103-105 in
polynomial.py:
# casting: if incoming array isn't floating point, make it floating
point.
if not isinstance(p.dtype, (NX.floating, NX.complexfloating)):
    p = p.astype(float)

which force p to be real, if it is complex, even though it appears
that they should leave things be if they are complex.

Andrew


	
		
__________________________________ 
Yahoo! Mail - PC Magazine Editors' Choice 2005 
http://mail.yahoo.com



More information about the SciPy-user mailing list