Mercurial > octave

changeset 19635:bdf90710dddf

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Map -pi to pi for principal argument used in complex operators (bug #43313). * oct-cmplx.h: For all operator templates in file, use conditional tests on arg() output to effectively map -pi to pi. * Array-C.cc (nan_ascending_compare, nan_descending_compare) : Update sort routines used by issorted to map -pi to pi. Use intermediate variables for a 40% speed-up in calculation. * Array-fC.cc (nan_ascending_compare, nan_descending_compare) : Update sort routines used by issorted to map -pi to pi. Use intermediate variables for a 40% speed-up in calculation. * test/complex.tst: New test file to check behavior. * test/Makefile.am: Add complex.tst to build system.
author Daniel J Sebald <daniel.sebald@ieee.org>
date Sat, 04 Oct 2014 02:17:17 -0500
parents d9c0b08e3da6
children 3da4b083e0b8
files liboctave/array/Array-C.cc liboctave/array/Array-fC.cc liboctave/util/oct-cmplx.h test/Makefile.am test/complex.tst
diffstat 5 files changed, 143 insertions(+), 17 deletions(-) [+]
line wrap: on
line diff
--- a/liboctave/array/Array-C.cc	Mon Jan 26 22:36:58 2015 -0800
+++ b/liboctave/array/Array-C.cc	Sat Oct 04 02:17:17 2014 -0500
@@ -41,22 +41,54 @@
   return xisnan (x);
 }
 
+// Sort Criteria: 1) isnan, 2) magnitude of z, 3) phase of z in range (-pi, pi]
+
 static bool
 nan_ascending_compare (const Complex& x, const Complex& y)
 {
-  return (xisnan (y)
-          ? ! xisnan (x)
-          : ((std::abs (x) < std::abs (y))
-             || ((std::abs (x) == std::abs (y)) && (arg (x) < arg (y)))));
+  if (xisnan (y))
+    return (! xisnan (x));
+
+  double xabs = std::abs (x);
+  double yabs = std::abs (y);
+
+  if (xabs < yabs)
+    return true;
+  else if (xabs == yabs)
+    {
+      double xarg = arg (x);
+      double yarg = arg (y);
+      xarg = (xarg == -M_PI) ? M_PI : xarg; 
+      yarg = (yarg == -M_PI) ? M_PI : yarg; 
+
+      return xarg < yarg;
+    }
+  else
+    return false;
 }
 
 static bool
 nan_descending_compare (const Complex& x, const Complex& y)
 {
-  return (xisnan (x)
-          ? ! xisnan (y)
-          : ((std::abs (x) > std::abs (y))
-             || ((std::abs (x) == std::abs (y)) && (arg (x) > arg (y)))));
+  if (xisnan (x))
+    return (! xisnan (y));
+
+  double xabs = std::abs (x);
+  double yabs = std::abs (y);
+
+  if (xabs > yabs)
+    return true;
+  else if (xabs == yabs)
+    {
+      double xarg = arg (x);
+      double yarg = arg (y);
+      xarg = (xarg == -M_PI) ? M_PI : xarg; 
+      yarg = (yarg == -M_PI) ? M_PI : yarg; 
+
+      return xarg > yarg;
+    }
+  else
+    return false;
 }
 
 Array<Complex>::compare_fcn_type
--- a/liboctave/array/Array-fC.cc	Mon Jan 26 22:36:58 2015 -0800
+++ b/liboctave/array/Array-fC.cc	Sat Oct 04 02:17:17 2014 -0500
@@ -41,22 +41,54 @@
   return xisnan (x);
 }
 
+// Sort Criteria: 1) isnan, 2) magnitude of z, 3) phase of z in range (-pi, pi]
+
 static bool
 nan_ascending_compare (const FloatComplex& x, const FloatComplex& y)
 {
-  return (xisnan (y)
-          ? ! xisnan (x)
-          : ((std::abs (x) < std::abs (x))
-             || ((std::abs (x) == std::abs (y)) && (arg (x) < arg (y)))));
+  if (xisnan (y))
+    return (! xisnan (x));
+
+  float xabs = std::abs (x);
+  float yabs = std::abs (y);
+
+  if (xabs < yabs)
+    return true;
+  else if (xabs == yabs)
+    {
+      float xarg = arg (x);
+      float yarg = arg (y);
+      xarg = (xarg == -M_PI) ? M_PI : xarg; 
+      yarg = (yarg == -M_PI) ? M_PI : yarg; 
+
+      return xarg < yarg;
+    }
+  else
+    return false;
 }
 
 static bool
 nan_descending_compare (const FloatComplex& x, const FloatComplex& y)
 {
-  return (xisnan (x)
-          ? ! xisnan (y)
-          : ((std::abs (x) > std::abs (x))
-             || ((std::abs (x) == std::abs (y)) && (arg (x) > arg (y)))));
+  if (xisnan (x))
+    return (! xisnan (y));
+
+  float xabs = std::abs (x);
+  float yabs = std::abs (y);
+
+  if (xabs > yabs)
+    return true;
+  else if (xabs == yabs)
+    {
+      float xarg = arg (x);
+      float yarg = arg (y);
+      xarg = (xarg == -M_PI) ? M_PI : xarg; 
+      yarg = (yarg == -M_PI) ? M_PI : yarg; 
+
+      return xarg > yarg;
+    }
+  else
+    return false;
 }
 
 Array<FloatComplex>::compare_fcn_type
--- a/liboctave/util/oct-cmplx.h	Mon Jan 26 22:36:58 2015 -0800
+++ b/liboctave/util/oct-cmplx.h	Sat Oct 04 02:17:17 2014 -0500
@@ -36,6 +36,10 @@
 // The abs/arg comparison is definitely more useful (the other one is emulated
 // rather trivially), so let's be consistent and use that all over.
 
+// The standard C library function arg() returns [-pi,pi], which creates a
+// non-unique representation for numbers along the negative real axis branch
+// cut.  Change this to principal value (-pi,pi] by mapping -pi to pi.
+
 #define DEF_COMPLEXR_COMP(OP, OPS) \
 template <class T> \
 inline bool operator OP (const std::complex<T>& a, const std::complex<T>& b) \
@@ -46,6 +50,15 @@
     { \
       FLOAT_TRUNCATE const T ay = std::arg (a); \
       FLOAT_TRUNCATE const T by = std::arg (b); \
+      if (ay == -M_PI) \
+        { \
+          if (by != -M_PI) \
+            return M_PI OP by; \
+        } \
+      else if (by == -M_PI) \
+        { \
+          return ay OP M_PI; \
+        } \
       return ay OP by; \
     } \
   else \
@@ -59,6 +72,8 @@
   if (ax == bx) \
     { \
       FLOAT_TRUNCATE const T ay = std::arg (a); \
+      if (ay == -M_PI) \
+        return M_PI OP 0; \
       return ay OP 0; \
     } \
   else \
@@ -72,6 +87,8 @@
   if (ax == bx) \
     { \
       FLOAT_TRUNCATE const T by = std::arg (b); \
+      if (by == -M_PI) \
+        return 0 OP M_PI; \
       return 0 OP by; \
     } \
   else \
--- a/test/Makefile.am	Mon Jan 26 22:36:58 2015 -0800
+++ b/test/Makefile.am	Sat Oct 04 02:17:17 2014 -0500
@@ -25,7 +25,8 @@
   args.tst \
   bug-31371.tst \
   bug-38576.tst \
- colormaps.tst \
+  colormaps.tst \
+  complex.tst \
   diag-perm.tst \
   error.tst \
   eval-catch.tst \
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/test/complex.tst	Sat Oct 04 02:17:17 2014 -0500
@@ -0,0 +1,44 @@
+## Copyright (C) 2015 Rik Wehbring
+##
+## This file is part of Octave.
+##
+## Octave is free software; you can redistribute it and/or modify it
+## under the terms of the GNU General Public License as published by
+## the Free Software Foundation; either version 3 of the License, or (at
+## your option) any later version.
+##
+## Octave is distributed in the hope that it will be useful, but
+## WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+## General Public License for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with Octave; see the file COPYING.  If not, see
+## <http://www.gnu.org/licenses/>.
+
+## Test ordering of complex values by magnitude and then by phase
+%!test
+%! x = [0 i 1+i 2 3i 3+4i];
+%! assert (sort (x, "descend"), fliplr (x));
+
+%!test
+%! x = [1, -1, i, -i];
+%! xs = [-i, 1, i, -1];
+%! assert (sort (x), xs);
+%! assert (sort (x, "descend"), fliplr (xs));
+
+## bug #44071, issorted incorrect because it uses different sort routine.
+%!assert (issorted ([1, -1, i, -i]), false)
+
+## bug #43313, -1 is both '>' and '==' to (-1 - 0i)
+%!test
+%! assert (complex(-1,0) == complex(-1,-0), true);
+%! assert (complex(-1,0) > complex(-1,-0), false);
+%! assert (complex(-1,0) < complex(-1,-0), false);
+
+## Test that sort and issorted both agree on boundary case
+%!test
+%! x = [complex(-1,0), complex(-1,-0), i, -i, 1];
+%! xs = sort (x);
+%! assert (issorted (xs));
+