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

[eric]

----- Original Message -----
From: "Pearu Peterson" <pearu at cens.ioc.ee>
To: <scipy-dev at scipy.org>
Sent: Saturday, February 16, 2002 5:22 AM
Subject: Re: [SciPy-dev] roots, poly, comparison operators, etc.


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

Thanks.  I've fixed this and checked it into the CVS.  So roots and poly are now
"rounding free" and seem to work.

eric