[SciPy-dev] roots, poly, comparison operators, etc.
eric at scipy.org
Sun Feb 17 01:08:59 CST 2002
----- 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
> 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,r
> 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,
Thanks. I've fixed this and checked it into the CVS. So roots and poly are now
"rounding free" and seem to work.
More information about the Scipy-dev