Mercurial > octave
changeset 24448:138c2380521c
Fix failing sort of complex values on Mac platforms (bug #52627).
* oct-cmplx.h: Add explanation of the unusual fact that Apple's float
version of the constant PI is different from other platforms. Use
#if defined with __APPLE__ to determine whether compilation is on a Mac.
For Mac platforms, provide a template specialization which overrides the
default template code in the case of floats.
| author | Rik <rik@octave.org> |
|---|---|
| date | Thu, 21 Dec 2017 14:30:38 -0800 |
| parents | b431e2719000 |
| children | 9a9b2e0596f5 |
| files | liboctave/util/oct-cmplx.h |
| diffstat | 1 files changed, 191 insertions(+), 60 deletions(-) [+] |
line wrap: on
line diff
--- a/liboctave/util/oct-cmplx.h Thu Dec 21 14:36:06 2017 -0500 +++ b/liboctave/util/oct-cmplx.h Thu Dec 21 14:30:38 2017 -0800 @@ -31,71 +31,202 @@ typedef std::complex<double> Complex; typedef std::complex<float> FloatComplex; -// For complex-complex and complex-real comparisons, we use the following -// ordering: compare absolute values first; if they match, compare phase angles. -// This is partially inconsistent with M*b, which compares complex numbers only -// by their real parts; OTOH, it uses the same definition for max/min and sort. -// 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. +// For complex-complex and complex-real comparisons, Octave uses the following +// ordering: compare absolute values first; if they match, compare phase +// angles. This is partially inconsistent with M*b, which compares complex +// numbers only by their real parts; OTOH, it uses the same definition for +// max/min and sort. 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 <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; \ - } \ - template <typename T> \ - inline bool operator OP (const std::complex<T>& a, 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); \ - if (ay == static_cast<T> (-M_PI)) \ - return static_cast<T> (M_PI) OP 0; \ - return ay OP 0; \ - } \ - else \ - return ax OPS bx; \ - } \ - template <typename T> \ - inline bool operator OP (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 by = std::arg (b); \ - if (by == static_cast<T> (-M_PI)) \ - return 0 OP static_cast<T> (M_PI); \ - return 0 OP by; \ - } \ - else \ - return ax OPS bx; \ - } +#if ! defined (__APPLE__) + + // General templated code for all (double, float) complex operators +# define DEF_COMPLEXR_COMP(OP, OPS) \ + 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; \ + } \ + template <typename T> \ + inline bool operator OP (const std::complex<T>& a, 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); \ + if (ay == static_cast<T> (-M_PI)) \ + return static_cast<T> (M_PI) OP 0; \ + return ay OP 0; \ + } \ + else \ + return ax OPS bx; \ + } \ + template <typename T> \ + inline bool operator OP (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 by = std::arg (b); \ + if (by == static_cast<T> (-M_PI)) \ + return 0 OP static_cast<T> (M_PI); \ + return 0 OP by; \ + } \ + else \ + return ax OPS bx; \ + } + +#else + // Apple specialization. + + // FIXME: Apple's math library chooses to return a different value for + // std::arg with floats than the constant static_cast<float> (M_PI). + // Work around this obtuse behavior by providing the constant A_PI which + // is Apple's definition of the constant PI for float variables. + // The template code for doubles is the same as that for UNIX platforms. + // Use C++ template specialization to add specific functions for float + // values that make use of A_PI. + + // Apple version of PI +# define A_PI 3.1415925025939941f + +# define DEF_COMPLEXR_COMP(OP, OPS) \ + 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; \ + } \ + template <typename T> \ + inline bool operator OP (const std::complex<T>& a, 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); \ + if (ay == static_cast<T> (-M_PI)) \ + return static_cast<T> (M_PI) OP 0; \ + return ay OP 0; \ + } \ + else \ + return ax OPS bx; \ + } \ + template <typename T> \ + inline bool operator OP (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 by = std::arg (b); \ + if (by == static_cast<T> (-M_PI)) \ + return 0 OP static_cast<T> (M_PI); \ + return 0 OP by; \ + } \ + else \ + return ax OPS bx; \ + } \ + template <> \ + inline bool operator OP<float> (const std::complex<float>& a, \ + const std::complex<float>& b) \ + { \ + OCTAVE_FLOAT_TRUNCATE const float ax = std::abs (a); \ + OCTAVE_FLOAT_TRUNCATE const float bx = std::abs (b); \ + if (ax == bx) \ + { \ + OCTAVE_FLOAT_TRUNCATE const float ay = std::arg (a); \ + OCTAVE_FLOAT_TRUNCATE const float by = std::arg (b); \ + if (ay == -A_PI) \ + { \ + if (by != -A_PI) \ + return static_cast<float> (M_PI) OP by; \ + } \ + else if (by == -A_PI) \ + { \ + return ay OP static_cast<float> (M_PI); \ + } \ + return ay OP by; \ + } \ + else \ + return ax OPS bx; \ + } \ + template <> \ + inline bool operator OP<float> (const std::complex<float>& a, float b) \ + { \ + OCTAVE_FLOAT_TRUNCATE const float ax = std::abs (a); \ + OCTAVE_FLOAT_TRUNCATE const float bx = std::abs (b); \ + if (ax == bx) \ + { \ + OCTAVE_FLOAT_TRUNCATE const float ay = std::arg (a); \ + if (ay == -A_PI) \ + return static_cast<float> (M_PI) OP 0; \ + return ay OP 0; \ + } \ + else \ + return ax OPS bx; \ + } \ + template <> \ + inline bool operator OP<float> (float a, const std::complex<float>& b) \ + { \ + OCTAVE_FLOAT_TRUNCATE const float ax = std::abs (a); \ + OCTAVE_FLOAT_TRUNCATE const float bx = std::abs (b); \ + if (ax == bx) \ + { \ + OCTAVE_FLOAT_TRUNCATE const float by = std::arg (b); \ + if (by == -A_PI) \ + return 0 OP static_cast<float> (M_PI); \ + return 0 OP by; \ + } \ + else \ + return ax OPS bx; \ + } + +#endif DEF_COMPLEXR_COMP (>, >) DEF_COMPLEXR_COMP (<, <)