Mercurial > octave-nkf
changeset 19669:bdf90710dddf
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)); +