Mercurial > octave
changeset 30668:965ead41658b
Avoid using M_PI in public header file (bug #61812).
* liboctave/util/oct-cmplx.h, liboctave/util/oct-cmplx.cc: Move implementation
of template from header to source file. It is using M_PI which is an XSI
extension that isn't defined in ISO C/C++.
* liboctave/util/module.mk: Add new file to build system.
| author | Markus Mützel <markus.muetzel@gmx.de> |
|---|---|
| date | Sun, 23 Jan 2022 21:20:13 +0100 |
| parents | f3e5f531e6b2 |
| children | 20b05c5c41a1 |
| files | liboctave/util/module.mk liboctave/util/oct-cmplx.cc liboctave/util/oct-cmplx.h |
| diffstat | 3 files changed, 280 insertions(+), 199 deletions(-) [+] |
line wrap: on
line diff
--- a/liboctave/util/module.mk Sun Jan 23 19:56:13 2022 -0800 +++ b/liboctave/util/module.mk Sun Jan 23 21:20:13 2022 +0100 @@ -73,6 +73,7 @@ %reldir%/quit.cc \ %reldir%/oct-atomic.c \ %reldir%/oct-base64.cc \ + %reldir%/oct-cmplx.cc \ %reldir%/oct-glob.cc \ %reldir%/oct-inttypes.cc \ %reldir%/oct-mutex.cc \
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/liboctave/util/oct-cmplx.cc Sun Jan 23 21:20:13 2022 +0100 @@ -0,0 +1,252 @@ +//////////////////////////////////////////////////////////////////////// +// +// Copyright (C) 1995-2022 The Octave Project Developers +// +// See the file COPYRIGHT.md in the top-level directory of this +// distribution or <https://octave.org/copyright/>. +// +// 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 +// <https://www.gnu.org/licenses/>. +// +//////////////////////////////////////////////////////////////////////// + +#if defined (HAVE_CONFIG_H) +# include "config.h" +#endif + +#include <complex> + +// 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. + +#if ! defined (__APPLE__) + + // General templated code for all (double, float) complex operators +# define DEF_COMPLEXR_COMP(OP, OPS) \ + template <typename T> \ + 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> \ + 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> \ + 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 in single precision +# define A_PI 3.1415925025939941f + +# define DEF_COMPLEXR_COMP(OP, OPS) \ + template <typename T> \ + 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> \ + 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> \ + 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 <> OCTAVE_API \ + /* specializations for float */ \ + 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 <> OCTAVE_API \ + 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 <> OCTAVE_API \ + 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 (<, <) +DEF_COMPLEXR_COMP (>=, >) +DEF_COMPLEXR_COMP (<=, <) + + +// Instantiate templates + +# define INST_COMPLEXR_COMP(OP, T) \ + template OCTAVE_API bool operator OP (const std::complex<T>& a, \ + const std::complex<T>& b); \ + template OCTAVE_API bool operator OP (const std::complex<T>& a, T b); \ + template OCTAVE_API bool operator OP (T a, const std::complex<T>& b); + +INST_COMPLEXR_COMP (>, double) +INST_COMPLEXR_COMP (<, double) +INST_COMPLEXR_COMP (>=, double) +INST_COMPLEXR_COMP (<=, double) + +#if ! defined (__APPLE__) +INST_COMPLEXR_COMP (>, float) +INST_COMPLEXR_COMP (<, float) +INST_COMPLEXR_COMP (>=, float) +INST_COMPLEXR_COMP (<=, float) +#endif
--- a/liboctave/util/oct-cmplx.h Sun Jan 23 19:56:13 2022 -0800 +++ b/liboctave/util/oct-cmplx.h Sun Jan 23 21:20:13 2022 +0100 @@ -33,206 +33,34 @@ typedef std::complex<double> Complex; typedef std::complex<float> FloatComplex; -// 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. - -#if ! defined (__APPLE__) +// template declaration +# define DECL_COMPLEXR_COMP(OP) \ + template <typename T> bool operator OP (const std::complex<T>& a, \ + const std::complex<T>& b); \ + template <typename T> bool operator OP (const std::complex<T>& a, T b); \ + template <typename T> bool operator OP (T a, const std::complex<T>& b); - // 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 +DECL_COMPLEXR_COMP(>) +DECL_COMPLEXR_COMP(<) +DECL_COMPLEXR_COMP(>=) +DECL_COMPLEXR_COMP(<=) -# 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; \ - } +// extern template instantiations +# define EXT_INST_COMPLEXR_COMP(OP, T) \ + extern template OCTAVE_API \ + bool operator OP (const std::complex<T>& a, const std::complex<T>& b); \ + extern template OCTAVE_API \ + bool operator OP (const std::complex<T>& a, T b); \ + extern template OCTAVE_API \ + bool operator OP (T a, const std::complex<T>& b); + +EXT_INST_COMPLEXR_COMP(>, double) +EXT_INST_COMPLEXR_COMP(>, float) +EXT_INST_COMPLEXR_COMP(<, double) +EXT_INST_COMPLEXR_COMP(<, float) +EXT_INST_COMPLEXR_COMP(>=, double) +EXT_INST_COMPLEXR_COMP(>=, float) +EXT_INST_COMPLEXR_COMP(<=, double) +EXT_INST_COMPLEXR_COMP(<=, float) #endif - -DEF_COMPLEXR_COMP (>, >) -DEF_COMPLEXR_COMP (<, <) -DEF_COMPLEXR_COMP (<=, <) -DEF_COMPLEXR_COMP (>=, >) - -#endif