GNU Octave  6.2.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
oct-cmplx.h
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 //
3 // Copyright (C) 1995-2021 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 (octave_oct_cmplx_h)
27 #define octave_oct_cmplx_h 1
28 
29 #include "octave-config.h"
30 
31 #include <complex>
32 
33 typedef std::complex<double> Complex;
34 typedef std::complex<float> FloatComplex;
35 
36 // For complex-complex and complex-real comparisons, Octave uses the following
37 // ordering: compare absolute values first; if they match, compare phase
38 // angles. This is partially inconsistent with M*b, which compares complex
39 // numbers only by their real parts; OTOH, it uses the same definition for
40 // max/min and sort. The abs/arg comparison is definitely more useful (the
41 // other one is emulated rather trivially), so let's be consistent and use that
42 // all over.
43 
44 // The standard C library function arg() returns [-pi,pi], which creates a
45 // non-unique representation for numbers along the negative real axis branch
46 // cut. Change this to principal value (-pi,pi] by mapping -pi to pi.
47 
48 #if ! defined (__APPLE__)
49 
50  // General templated code for all (double, float) complex operators
51 # define DEF_COMPLEXR_COMP(OP, OPS) \
52  template <typename T> \
53  inline bool operator OP (const std::complex<T>& a, \
54  const std::complex<T>& b) \
55  { \
56  OCTAVE_FLOAT_TRUNCATE const T ax = std::abs (a); \
57  OCTAVE_FLOAT_TRUNCATE const T bx = std::abs (b); \
58  if (ax == bx) \
59  { \
60  OCTAVE_FLOAT_TRUNCATE const T ay = std::arg (a); \
61  OCTAVE_FLOAT_TRUNCATE const T by = std::arg (b); \
62  if (ay == static_cast<T> (-M_PI)) \
63  { \
64  if (by != static_cast<T> (-M_PI)) \
65  return static_cast<T> (M_PI) OP by; \
66  } \
67  else if (by == static_cast<T> (-M_PI)) \
68  { \
69  return ay OP static_cast<T> (M_PI); \
70  } \
71  return ay OP by; \
72  } \
73  else \
74  return ax OPS bx; \
75  } \
76  template <typename T> \
77  inline bool operator OP (const std::complex<T>& a, T b) \
78  { \
79  OCTAVE_FLOAT_TRUNCATE const T ax = std::abs (a); \
80  OCTAVE_FLOAT_TRUNCATE const T bx = std::abs (b); \
81  if (ax == bx) \
82  { \
83  OCTAVE_FLOAT_TRUNCATE const T ay = std::arg (a); \
84  if (ay == static_cast<T> (-M_PI)) \
85  return static_cast<T> (M_PI) OP 0; \
86  return ay OP 0; \
87  } \
88  else \
89  return ax OPS bx; \
90  } \
91  template <typename T> \
92  inline bool operator OP (T a, const std::complex<T>& b) \
93  { \
94  OCTAVE_FLOAT_TRUNCATE const T ax = std::abs (a); \
95  OCTAVE_FLOAT_TRUNCATE const T bx = std::abs (b); \
96  if (ax == bx) \
97  { \
98  OCTAVE_FLOAT_TRUNCATE const T by = std::arg (b); \
99  if (by == static_cast<T> (-M_PI)) \
100  return 0 OP static_cast<T> (M_PI); \
101  return 0 OP by; \
102  } \
103  else \
104  return ax OPS bx; \
105  }
106 
107 #else
108  // Apple specialization.
109 
110  // FIXME: Apple's math library chooses to return a different value for
111  // std::arg with floats than the constant static_cast<float> (M_PI).
112  // Work around this obtuse behavior by providing the constant A_PI which
113  // is Apple's definition of the constant PI for float variables.
114  // The template code for doubles is the same as that for UNIX platforms.
115  // Use C++ template specialization to add specific functions for float
116  // values that make use of A_PI.
117 
118  // Apple version of PI
119 # define A_PI 3.1415925025939941f
120 
121 # define DEF_COMPLEXR_COMP(OP, OPS) \
122  template <typename T> \
123  inline bool operator OP (const std::complex<T>& a, \
124  const std::complex<T>& b) \
125  { \
126  OCTAVE_FLOAT_TRUNCATE const T ax = std::abs (a); \
127  OCTAVE_FLOAT_TRUNCATE const T bx = std::abs (b); \
128  if (ax == bx) \
129  { \
130  OCTAVE_FLOAT_TRUNCATE const T ay = std::arg (a); \
131  OCTAVE_FLOAT_TRUNCATE const T by = std::arg (b); \
132  if (ay == static_cast<T> (-M_PI)) \
133  { \
134  if (by != static_cast<T> (-M_PI)) \
135  return static_cast<T> (M_PI) OP by; \
136  } \
137  else if (by == static_cast<T> (-M_PI)) \
138  { \
139  return ay OP static_cast<T> (M_PI); \
140  } \
141  return ay OP by; \
142  } \
143  else \
144  return ax OPS bx; \
145  } \
146  template <typename T> \
147  inline bool operator OP (const std::complex<T>& a, T b) \
148  { \
149  OCTAVE_FLOAT_TRUNCATE const T ax = std::abs (a); \
150  OCTAVE_FLOAT_TRUNCATE const T bx = std::abs (b); \
151  if (ax == bx) \
152  { \
153  OCTAVE_FLOAT_TRUNCATE const T ay = std::arg (a); \
154  if (ay == static_cast<T> (-M_PI)) \
155  return static_cast<T> (M_PI) OP 0; \
156  return ay OP 0; \
157  } \
158  else \
159  return ax OPS bx; \
160  } \
161  template <typename T> \
162  inline bool operator OP (T a, const std::complex<T>& b) \
163  { \
164  OCTAVE_FLOAT_TRUNCATE const T ax = std::abs (a); \
165  OCTAVE_FLOAT_TRUNCATE const T bx = std::abs (b); \
166  if (ax == bx) \
167  { \
168  OCTAVE_FLOAT_TRUNCATE const T by = std::arg (b); \
169  if (by == static_cast<T> (-M_PI)) \
170  return 0 OP static_cast<T> (M_PI); \
171  return 0 OP by; \
172  } \
173  else \
174  return ax OPS bx; \
175  } \
176  template <> \
177  inline bool operator OP<float> (const std::complex<float>& a, \
178  const std::complex<float>& b) \
179  { \
180  OCTAVE_FLOAT_TRUNCATE const float ax = std::abs (a); \
181  OCTAVE_FLOAT_TRUNCATE const float bx = std::abs (b); \
182  if (ax == bx) \
183  { \
184  OCTAVE_FLOAT_TRUNCATE const float ay = std::arg (a); \
185  OCTAVE_FLOAT_TRUNCATE const float by = std::arg (b); \
186  if (ay == -A_PI) \
187  { \
188  if (by != -A_PI) \
189  return static_cast<float> (M_PI) OP by; \
190  } \
191  else if (by == -A_PI) \
192  { \
193  return ay OP static_cast<float> (M_PI); \
194  } \
195  return ay OP by; \
196  } \
197  else \
198  return ax OPS bx; \
199  } \
200  template <> \
201  inline bool operator OP<float> (const std::complex<float>& a, float b) \
202  { \
203  OCTAVE_FLOAT_TRUNCATE const float ax = std::abs (a); \
204  OCTAVE_FLOAT_TRUNCATE const float bx = std::abs (b); \
205  if (ax == bx) \
206  { \
207  OCTAVE_FLOAT_TRUNCATE const float ay = std::arg (a); \
208  if (ay == -A_PI) \
209  return static_cast<float> (M_PI) OP 0; \
210  return ay OP 0; \
211  } \
212  else \
213  return ax OPS bx; \
214  } \
215  template <> \
216  inline bool operator OP<float> (float a, const std::complex<float>& b) \
217  { \
218  OCTAVE_FLOAT_TRUNCATE const float ax = std::abs (a); \
219  OCTAVE_FLOAT_TRUNCATE const float bx = std::abs (b); \
220  if (ax == bx) \
221  { \
222  OCTAVE_FLOAT_TRUNCATE const float by = std::arg (b); \
223  if (by == -A_PI) \
224  return 0 OP static_cast<float> (M_PI); \
225  return 0 OP by; \
226  } \
227  else \
228  return ax OPS bx; \
229  }
230 
231 #endif
232 
233 DEF_COMPLEXR_COMP (>, >)
234 DEF_COMPLEXR_COMP (<, <)
235 DEF_COMPLEXR_COMP (<=, <)
236 DEF_COMPLEXR_COMP (>=, >)
237 
238 #endif
Complex
std::complex< double > Complex
Definition: oct-cmplx.h:33
DEF_COMPLEXR_COMP
#define DEF_COMPLEXR_COMP(OP, OPS)
Definition: oct-cmplx.h:51
FloatComplex
std::complex< float > FloatComplex
Definition: oct-cmplx.h:34