[Octave-bug-tracker] [bug #52919] comparison of complex values

[Rik]

Follow-up Comment #10, bug #52919 (project octave):

This is not an easy thing to address.  It should probably be discussed back on
the Octave Maintainer's list if there is going to be a major change to
ordering.

First, I agree with principle of least surprise.  To begin with that means
operators and max, min, sort should all behave the same way.  Octave is not
quite there yet.  The operators and sort are self-consistent.  And max/min are
consistent when called with one argument.  When called with a two argument
form such as max (x, 1+1i) the ordering is not consistent.  But we already
have a bug report for that.  See bug #51307.

One thing that we are working against is Matlab's definition of complex. 
iscomplex() does not mean the imaginary part is non-zero, it merely indicates
that Matlab has reserved storage for a complex portion of the number.  I know
that it may be tempting to think of real numbers as the limiting case of
complex numbers as the imaginary part goes to zero.


xreal =   limit  (x + epsilon*i)
        epsilion -> 0


but this isn't what Matlab, and hence Octave, have implemented.  There is
always a categorical distinction between a real number and a complex number.

So Octave has two categories of numbers and applies a different ordering
system for each number type.  Incidentally, I believe Matlab used to sort real
numbers by value and complex numbers by abs() without checking the phase angle
for ties.  This is a horrible idea because the same program can produce
different results based on the input data.  In the case of sort, equal
magnitude numbers might appear in the output in the order they were in the
input, or frankly any other permutation.  Not only does this create programs
which fail based on differences in input data, but given that the output could
depend on the internals of the sort algorithm, it also introduced a dependence
on which version of Matlab you were using.  At least Matlab now uses all
information (magnitude and phase angle) to unambiguously order things. 
According to the documentation for sort


If A is complex, then by default, sort sorts the elements by magnitude. If
more than one element has equal magnitude, then the elements are sorted by
phase angle on the interval (−π, π].


This already seems to be what Octave does.


x = [-1+1i, 0, 1+1i]
x =

  -1 + 1i   0 + 0i   1 + 1i

angle (x)
ans =

   2.35619   0.00000   0.78540

sort (x)
ans =

   0 + 0i   1 + 1i  -1 + 1i


Hooray, Octave is Matlab-compatible as the number -1+1i with phase angle 2.356
is placed after 1+1i which only has a phase angle of 0.7854.

Octave could just stop there since we are compatible, but we do try to be
better than Matlab where we can.

I agree that it is somewhat galling to see 1+1i ahead of -1+1i.  Intuitively,
it feels like the number with a negative real part should be smaller.  In a
perfect world, -1 + 1i is less than 1+1i.  In addition, I understand that it
is counter-intuitive that adding a complex number to an existing complex
number can reverse the sort order.


x = -1 + 1i;
y = 1+1i;
z = -2i;
x < y
ans = 0
(x+z) < (y+z)
ans = 1


One possibility would be to do a slightly more sophisticated ordering on the
phase angle rather than just using the value.  If we first resolved whether
the angles were in the left half of the complex plane (corresponding to
negative real numbers) versus the right (corresponding to positive real
numbers) then at least we would have an ordering that would be consistent with
the real numbers as the imaginary part of the number was shrunk to zero.  For
the test example above, the result of the '<' operator would also not flip
based on the addition of an imaginary number.  That would be another nice
gain.

Here is a sample of 4 relevant numbers with equal magnitude, but different
phase angles.


x = [-1-1i, -1+1i, 1-1i, 1+1i]
x =

  -1 - 1i  -1 + 1i   1 - 1i   1 + 1i

octave:51> angle (x)
ans =

  -2.35619   2.35619  -0.78540   0.78540


That seems very distinguishable.  There may be an efficient way to calculate
this without even calling the angle() function and relying on the signbits and
magnitudes of the real and imaginary parts directly.

For reference, the code that implements operators is in
liboctave/util/oct-cmplx.h


     template <typename T>                                                  
\
     inline bool operator OP (const std::complex<T>& a,                     
\
                              const std::complex<T>& b)                     
\
     {                                                                      
\
       OCTAVE_FLOAT_TRUNCATE const T ax = std::abs (a);                     
\
       OCTAVE_FLOAT_TRUNCATE const T bx = std::abs (b);                     
\
       if (ax == bx)                                                        
\
         {                                                                  
\
           OCTAVE_FLOAT_TRUNCATE const T ay = std::arg (a);                 
\
           OCTAVE_FLOAT_TRUNCATE const T by = std::arg (b);                 
\
           if (ay == static_cast<T> (-M_PI))                                
\
             {                                                              
\
               if (by != static_cast<T> (-M_PI))                            
\
                 return static_cast<T> (M_PI) OP by;                        
\
             }                                                              
\
           else if (by == static_cast<T> (-M_PI))                           
\
             {                                                              
\
               return ay OP static_cast<T> (M_PI);                          
\
             }                                                              
\
           return ay OP by;                                                 
\
         }                                                                  
\
       else                                                                 
\
         return ax OPS bx;                                                  
\
     }                                                                      
\


I don't know how expensive the std::arg call is, but Octave is already doing a
bit of work to re-map -PI to +PI that we might be able to do away with.

In pseduocode, after determining that the magnitudes are equal.


  areal = real (a);
  breal = real (b);
  if (areal != breal)
    return areal OP breal;   # could also just use signbit for comparison

  aimag = imag (a);
  bimag = imag (b);
  if (aimag != bimag)
    return aimag OP bimag;   # could also just use signbit for comparison
  
  ## In the same quadrant
  return (aimag / areal) OP (bimag / breal);




    _______________________________________________________

Reply to this item at:

  <http://savannah.gnu.org/bugs/?52919>

_______________________________________________
  Message sent via/by Savannah
  http://savannah.gnu.org/