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

[Pearu Peterson]

Hi!

On Sat, 16 Feb 2002, eric wrote:

> Other things:
> 
> * Array printing is far superior in Matlab and Octave -- they
> generally always look nice.  We should clean up how arrays are output.  
> Also, the "format long" and "format short", etc options for specifying
> how arrays are printed are pretty nice.

I agree. May be ipython could be exploited here?

In fact, rounding of tiny numbers to zeros could be done only when
printing (though, personally I wouldn't prefer that either but I am
just looking for a compromise), not inside calculation routines. In this
way, no valuable information is lost when using these routines from
other calculation routines and computation will be even more efficient.

> * On Matlab/Octave, sort() sorts by magnitude of complex and then by
> angle.  On the other hand ==, >, <, etc. seem to only compare the real
> part of complex numbers.  
> These seem fine to me.  I know they aren't mathematically correct, but
> they seem to be pragmatically correct.  I'd like comments on these
> conventions and what others think.

There is no mathematically correct way to compare complex numbers,
they just cannot be ordered in an unique and sensible way.

However, in different applications different conventions may be useful or
reasonable for ordering complex numbers. Whatever is the convention, their
mathematical correctness is irrelevant and this cannot be used as an
argument for prefering one convention to another.

I would propose providing number of efficient comparison methods for
complex (or any) numbers that users may use in sort functions as an
optional argument. For example,

scipy.sort([2,1+2j],cmpmth='abs') -> [1+2j,2]  # sorts by abs value
scipy.sort([2,1+2j],cmpmth='real') -> [2,1+2j] # sorts by real part
scipy.sort([2,1+2j],cmpmth='realimag') # sorts by real then by imag
scipy.sort([2,1+2j],cmpmth='imagreal') # sorts by imag then by real
scipy.sort([2,1+2j],cmpmth='absangle') # sorts by abs then by angle
etc.
scipy.sort([2,1+2j],cmpfunc=<user defined comparison function>)

Note that

scipy.sort([-1,1],cmpmth='absangle') -> [1,-1]

which also demonstrates the arbitrariness of sorting complex numbers.

Btw, why do you want to sort the output of roots()? As far as I know,
there is no order defined for roots of polynomials. May be an optional
flag could be provided here?

> * Comparison of IEEE floats is hairy.  It looks to me like Matlab and
> Octave have chosen to limit precision to 15 digits.  This seems like a
> reasonable thing to do, for SciPy also, but we'd currently have to
> limit to 14 digits to deal with problems of imprecise LAPACK routines.  
> Pearu and Fernando are arguing against this, but maybe they wouldn't
> mind a 15 digit limit for double and 6 for float.  We'd have to modify
> Numeric to do this internally on comparisons.  There could be a flag
> that enables exact comparisons.

I hope this issue is solved by fixing roots() and also the accuracy 
of LAPACK routines can be rehabilitated now (I don't claim that their
output is always accurate, just floating numbers cannot be represented
accuratelly in computers memory, in principle. It has nothing to do with
the programs that manipulate these numbers, one can only improve the
algorithtms to minimize the computational errors).

The deep lesson to learn here is:
  Fix the computation algorithm.
  Do not fix the output of the computation.

Matlab and Octave are programs that are oriented to certain applications,
namely, to engineering applications where linear algebra is very
important.

SciPy needs not to choose the same orientation as Matlab, however, it can
certainly cover Matlab's orientations.

Python is a general purpose language and SciPy could also be a general
purpose package for scientific computing (whatever are the applications).

Regards,
	Pearu