GNU Octave 7.1.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-2022 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
33typedef std::complex<double> Complex;
34typedef 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
233DEF_COMPLEXR_COMP (>, >)
234DEF_COMPLEXR_COMP (<, <)
235DEF_COMPLEXR_COMP (<=, <)
236DEF_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