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