[SciPy-dev] roots, poly, comparison operators, etc.

[Pearu Peterson]

On Sat, 16 Feb 2002, eric wrote:

> Hey crew,
> 
> I spent some quality time with roots() and poly() tonight trying to find
> a solution to the rounding issue acceptable to all.  In the process,I've
> rewritten both of them.  As Heiko suggested, it is easy to test if the
> roots are complex congugates and reals in poly(), and, if so, ignore the
> imaginary part.  Things are never as easy as they seem...    

> Here is what I have learned:
> 
> roots() calculates eigenvalues under the covers.  This is done using eig
> which in turn calls some LAPACK (Atlas) routine xgeev, where x is the
> appropriate numeric type.  This is where the trouble starts.  The values
> returned by eig() that should be complex conjugates are actually off by
> as much as 2-3 digits of precision.  I've compared the output of SciPy
> to Matlab and Octave below.  

scipy.eig uses Fortran LAPACK (octave uses also lapack, I think(??)).

Actually the trouble starts from roots that forces scipy.eig to use
flapack.zgeev though flapack.dgeev would be more appropiate (real
input,more efficient). Namely, the function roots() uses complex array 

  A = diag(Numeric.ones((N-2,),'D'),-1)

for real input array

  p

If I make a change

  A = diag(Numeric.ones((N-2,),'d'),-1)

then complex conjugates will match exactly:

>>> r=scipy.roots([-1,1,2,1,2,3,4,5,6,7,8,9])
>>> r[0],r[1]
((-0.20505419046158363+1.0643579627494124j),
(-0.20505419046158363-1.0643579627494124j))

I believe that also Matlab and octave use the most appropiate routines in
their root calculations. So, scipy.roots should be fixed so that for real
and complex inputs scipy.eig will use flapack.dgeev and flapacl.zgeev,
respectively.

Regards,
	Pearu