GNU Octave 10.1.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
 
Loading...
Searching...
No Matches
oct-cmplx.cc
Go to the documentation of this file.
1////////////////////////////////////////////////////////////////////////
2//
3// Copyright (C) 1995-2025 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 const T ax = std::abs (a); \
50 const T bx = std::abs (b); \
51 if (ax == bx) \
52 { \
53 const T ay = std::arg (a); \
54 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 const T ax = std::abs (a); \
73 const T bx = std::abs (b); \
74 if (ax == bx) \
75 { \
76 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 const T ax = std::abs (a); \
88 const T bx = std::abs (b); \
89 if (ax == bx) \
90 { \
91 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 const float ax = std::abs (a); \
120 const float bx = std::abs (b); \
121 if (ax == bx) \
122 { \
123 const float ay = std::arg (a); \
124 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 const float ax = std::abs (a); \
143 const float bx = std::abs (b); \
144 if (ax == bx) \
145 { \
146 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 const float ax = std::abs (a); \
158 const float bx = std::abs (b); \
159 if (ax == bx) \
160 { \
161 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
172DEF_COMPLEXR_COMP (>, >)
173DEF_COMPLEXR_COMP (<, <)
174DEF_COMPLEXR_COMP (>=, >)
175DEF_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
186INST_COMPLEXR_COMP (>, double)
187INST_COMPLEXR_COMP (<, double)
188INST_COMPLEXR_COMP (>=, double)
189INST_COMPLEXR_COMP (<=, double)
190
191#if defined (__APPLE__)
192INST_SPECIAL_COMPLEXR_COMP (>, >)
193INST_SPECIAL_COMPLEXR_COMP (<, <)
194INST_SPECIAL_COMPLEXR_COMP (>=, >)
195INST_SPECIAL_COMPLEXR_COMP (<=, <)
196#else
197INST_COMPLEXR_COMP (>, float)
198INST_COMPLEXR_COMP (<, float)
199INST_COMPLEXR_COMP (>=, float)
200INST_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