GNU Octave  9.1.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
oct-cmplx.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 //
3 // Copyright (C) 1995-2024 The Octave Project Developers
4 //
5 // See the file COPYRIGHT.md in the top-level directory of this
6 // distribution or <https://octave.org/copyright/>.
7 //
8 // This file is part of Octave.
9 //
10 // Octave is free software: you can redistribute it and/or modify it
11 // under the terms of the GNU General Public License as published by
12 // the Free Software Foundation, either version 3 of the License, or
13 // (at your option) any later version.
14 //
15 // Octave is distributed in the hope that it will be useful, but
16 // WITHOUT ANY WARRANTY; without even the implied warranty of
17 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18 // GNU General Public License for more details.
19 //
20 // You should have received a copy of the GNU General Public License
21 // along with Octave; see the file COPYING. If not, see
22 // <https://www.gnu.org/licenses/>.
23 //
24 ////////////////////////////////////////////////////////////////////////
25 
26 #if defined (HAVE_CONFIG_H)
27 # include "config.h"
28 #endif
29 
30 #include <complex>
31 
32 // For complex-complex and complex-real comparisons, Octave uses the following
33 // ordering: compare absolute values first; if they match, compare phase
34 // angles. This is partially inconsistent with M*b, which compares complex
35 // numbers only by their real parts; OTOH, it uses the same definition for
36 // max/min and sort. The abs/arg comparison is definitely more useful (the
37 // other one is emulated rather trivially), so let's be consistent and use that
38 // all over.
39 
40 // The standard C library function arg() returns [-pi,pi], which creates a
41 // non-unique representation for numbers along the negative real axis branch
42 // cut. Change this to principal value (-pi,pi] by mapping -pi to pi.
43 
44 // General templated code for all (double, float) complex operators
45 #define DEF_COMPLEXR_COMP(OP, OPS) \
46  template <typename T> \
47  bool operator OP (const std::complex<T>& a, const std::complex<T>& b) \
48  { \
49  OCTAVE_FLOAT_TRUNCATE const T ax = std::abs (a); \
50  OCTAVE_FLOAT_TRUNCATE const T bx = std::abs (b); \
51  if (ax == bx) \
52  { \
53  OCTAVE_FLOAT_TRUNCATE const T ay = std::arg (a); \
54  OCTAVE_FLOAT_TRUNCATE const T by = std::arg (b); \
55  if (ay == static_cast<T> (-M_PI)) \
56  { \
57  if (by != static_cast<T> (-M_PI)) \
58  return static_cast<T> (M_PI) OP by; \
59  } \
60  else if (by == static_cast<T> (-M_PI)) \
61  { \
62  return ay OP static_cast<T> (M_PI); \
63  } \
64  return ay OP by; \
65  } \
66  else \
67  return ax OPS bx; \
68  } \
69  template <typename T> \
70  bool operator OP (const std::complex<T>& a, T b) \
71  { \
72  OCTAVE_FLOAT_TRUNCATE const T ax = std::abs (a); \
73  OCTAVE_FLOAT_TRUNCATE const T bx = std::abs (b); \
74  if (ax == bx) \
75  { \
76  OCTAVE_FLOAT_TRUNCATE const T ay = std::arg (a); \
77  if (ay == static_cast<T> (-M_PI)) \
78  return static_cast<T> (M_PI) OP 0; \
79  return ay OP 0; \
80  } \
81  else \
82  return ax OPS bx; \
83  } \
84  template <typename T> \
85  bool operator OP (T a, const std::complex<T>& b) \
86  { \
87  OCTAVE_FLOAT_TRUNCATE const T ax = std::abs (a); \
88  OCTAVE_FLOAT_TRUNCATE const T bx = std::abs (b); \
89  if (ax == bx) \
90  { \
91  OCTAVE_FLOAT_TRUNCATE const T by = std::arg (b); \
92  if (by == static_cast<T> (-M_PI)) \
93  return 0 OP static_cast<T> (M_PI); \
94  return 0 OP by; \
95  } \
96  else \
97  return ax OPS bx; \
98  }
99 
100 #if defined (__APPLE__)
101  // Apple specializations
102 
103  // FIXME: Apple's math library chooses to return a different value for
104  // std::arg with floats than the constant static_cast<float> (M_PI).
105  // Work around this obtuse behavior by providing the constant A_PI which
106  // is Apple's definition of the constant PI for float variables.
107  // The template code for doubles is the same as that for UNIX platforms.
108  // Use C++ template specialization to add specific functions for float
109  // values that make use of A_PI.
110 
111  // Apple version of PI in single precision
112 # define A_PI 3.1415925025939941f
113 
114 # define INST_SPECIAL_COMPLEXR_COMP(OP, OPS) \
115  template <> OCTAVE_API \
116  bool operator OP<float> (const std::complex<float>& a, \
117  const std::complex<float>& b) \
118  { \
119  OCTAVE_FLOAT_TRUNCATE const float ax = std::abs (a); \
120  OCTAVE_FLOAT_TRUNCATE const float bx = std::abs (b); \
121  if (ax == bx) \
122  { \
123  OCTAVE_FLOAT_TRUNCATE const float ay = std::arg (a); \
124  OCTAVE_FLOAT_TRUNCATE const float by = std::arg (b); \
125  if (ay == -A_PI) \
126  { \
127  if (by != -A_PI) \
128  return static_cast<float> (M_PI) OP by; \
129  } \
130  else if (by == -A_PI) \
131  { \
132  return ay OP static_cast<float> (M_PI); \
133  } \
134  return ay OP by; \
135  } \
136  else \
137  return ax OPS bx; \
138  } \
139  template <> OCTAVE_API \
140  bool operator OP<float> (const std::complex<float>& a, float b) \
141  { \
142  OCTAVE_FLOAT_TRUNCATE const float ax = std::abs (a); \
143  OCTAVE_FLOAT_TRUNCATE const float bx = std::abs (b); \
144  if (ax == bx) \
145  { \
146  OCTAVE_FLOAT_TRUNCATE const float ay = std::arg (a); \
147  if (ay == -A_PI) \
148  return static_cast<float> (M_PI) OP 0; \
149  return ay OP 0; \
150  } \
151  else \
152  return ax OPS bx; \
153  } \
154  template <> OCTAVE_API \
155  bool operator OP<float> (float a, const std::complex<float>& b) \
156  { \
157  OCTAVE_FLOAT_TRUNCATE const float ax = std::abs (a); \
158  OCTAVE_FLOAT_TRUNCATE const float bx = std::abs (b); \
159  if (ax == bx) \
160  { \
161  OCTAVE_FLOAT_TRUNCATE const float by = std::arg (b); \
162  if (by == -A_PI) \
163  return 0 OP static_cast<float> (M_PI); \
164  return 0 OP by; \
165  } \
166  else \
167  return ax OPS bx; \
168  }
169 
170 #endif
171 
172 DEF_COMPLEXR_COMP (>, >)
173 DEF_COMPLEXR_COMP (<, <)
174 DEF_COMPLEXR_COMP (>=, >)
175 DEF_COMPLEXR_COMP (<=, <)
176 
177 
178 // Instantiate templates
179 
180 # define INST_COMPLEXR_COMP(OP, T) \
181  template OCTAVE_API bool operator OP (const std::complex<T>& a, \
182  const std::complex<T>& b); \
183  template OCTAVE_API bool operator OP (const std::complex<T>& a, T b); \
184  template OCTAVE_API bool operator OP (T a, const std::complex<T>& b);
185 
186 INST_COMPLEXR_COMP (>, double)
187 INST_COMPLEXR_COMP (<, double)
188 INST_COMPLEXR_COMP (>=, double)
189 INST_COMPLEXR_COMP (<=, double)
190 
191 #if defined (__APPLE__)
192 INST_SPECIAL_COMPLEXR_COMP (>, >)
193 INST_SPECIAL_COMPLEXR_COMP (<, <)
194 INST_SPECIAL_COMPLEXR_COMP (>=, >)
195 INST_SPECIAL_COMPLEXR_COMP (<=, <)
196 #else
197 INST_COMPLEXR_COMP (>, float)
198 INST_COMPLEXR_COMP (<, float)
199 INST_COMPLEXR_COMP (>=, float)
200 INST_COMPLEXR_COMP (<=, float)
201 #endif
DEF_COMPLEXR_COMP
#define DEF_COMPLEXR_COMP(OP, OPS)
Definition: oct-cmplx.cc:45
INST_COMPLEXR_COMP
#define INST_COMPLEXR_COMP(OP, T)
Definition: oct-cmplx.cc:180