Mercurial > hg > octave-nkf
annotate src/xpow.cc @ 8978:0a58c4cd1405
optimize element-wise power
author | Jaroslav Hajek <highegg@gmail.com> |
---|---|
date | Fri, 13 Mar 2009 22:33:26 +0100 |
parents | 93f18f166aba |
children | a7c00773a089 |
rev | line source |
---|---|
1 | 1 /* |
2 | |
7017 | 3 Copyright (C) 1993, 1994, 1995, 1996, 1997, 1998, 2000, 2002, 2003, |
8920 | 4 2004, 2005, 2006, 2007, 2008 John W. Eaton |
1 | 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 the | |
7016 | 10 Free Software Foundation; either version 3 of the License, or (at your |
11 option) any later version. | |
1 | 12 |
13 Octave is distributed in the hope that it will be useful, but WITHOUT | |
14 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or | |
15 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License | |
16 for more details. | |
17 | |
18 You should have received a copy of the GNU General Public License | |
7016 | 19 along with Octave; see the file COPYING. If not, see |
20 <http://www.gnu.org/licenses/>. | |
1 | 21 |
22 */ | |
23 | |
240 | 24 #ifdef HAVE_CONFIG_H |
1192 | 25 #include <config.h> |
1 | 26 #endif |
27 | |
1343 | 28 #include <cassert> |
1580 | 29 #include <climits> |
1343 | 30 |
4669 | 31 #include "Array-util.h" |
1352 | 32 #include "CColVector.h" |
453 | 33 #include "CDiagMatrix.h" |
8382
9b20a4847056
implement scalar powers of diag matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
34 #include "fCDiagMatrix.h" |
1352 | 35 #include "CMatrix.h" |
453 | 36 #include "EIG.h" |
7789
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
37 #include "fEIG.h" |
1352 | 38 #include "dDiagMatrix.h" |
8382
9b20a4847056
implement scalar powers of diag matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
39 #include "fDiagMatrix.h" |
1352 | 40 #include "dMatrix.h" |
8958
6ccc12cc65ef
implement raising a permutation matrix to integer power
Jaroslav Hajek <highegg@gmail.com>
parents:
8920
diff
changeset
|
41 #include "PermMatrix.h" |
3585 | 42 #include "mx-cm-cdm.h" |
1651 | 43 #include "oct-cmplx.h" |
4153 | 44 #include "quit.h" |
1352 | 45 |
46 #include "error.h" | |
4055 | 47 #include "oct-obj.h" |
1567 | 48 #include "utils.h" |
1352 | 49 #include "xpow.h" |
1 | 50 |
5275 | 51 #ifdef _OPENMP |
52 #include <omp.h> | |
53 #endif | |
54 | |
1567 | 55 static inline int |
56 xisint (double x) | |
57 { | |
58 return (D_NINT (x) == x | |
59 && ((x >= 0 && x < INT_MAX) | |
60 || (x <= 0 && x > INT_MIN))); | |
61 } | |
62 | |
767 | 63 // Safer pow functions. |
64 // | |
65 // op2 \ op1: s m cs cm | |
66 // +-- +---+---+----+----+ | |
67 // scalar | | 1 | 5 | 7 | 11 | | |
68 // +---+---+----+----+ | |
2365 | 69 // matrix | 2 | * | 8 | * | |
767 | 70 // +---+---+----+----+ |
71 // complex_scalar | 3 | 6 | 9 | 12 | | |
72 // +---+---+----+----+ | |
2365 | 73 // complex_matrix | 4 | * | 10 | * | |
767 | 74 // +---+---+----+----+ |
1 | 75 |
767 | 76 // -*- 1 -*- |
2086 | 77 octave_value |
1 | 78 xpow (double a, double b) |
79 { | |
7543
b84c5cbc0812
print.m: handle gif and jpg devices
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
80 double retval; |
b84c5cbc0812
print.m: handle gif and jpg devices
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
81 |
8978
0a58c4cd1405
optimize element-wise power
Jaroslav Hajek <highegg@gmail.com>
parents:
8960
diff
changeset
|
82 if (a < 0.0 && ! xisint (b)) |
1 | 83 { |
5260 | 84 Complex atmp (a); |
4682 | 85 |
5260 | 86 return std::pow (atmp, b); |
1 | 87 } |
88 else | |
7543
b84c5cbc0812
print.m: handle gif and jpg devices
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
89 retval = std::pow (a, b); |
b84c5cbc0812
print.m: handle gif and jpg devices
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
90 |
b84c5cbc0812
print.m: handle gif and jpg devices
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
91 return retval; |
1 | 92 } |
93 | |
767 | 94 // -*- 2 -*- |
2086 | 95 octave_value |
164 | 96 xpow (double a, const Matrix& b) |
1 | 97 { |
2086 | 98 octave_value retval; |
1 | 99 |
5275 | 100 octave_idx_type nr = b.rows (); |
101 octave_idx_type nc = b.cols (); | |
1 | 102 |
103 if (nr == 0 || nc == 0 || nr != nc) | |
104 error ("for x^A, A must be square"); | |
105 else | |
106 { | |
107 EIG b_eig (b); | |
108 | |
7578
91f8446ce4ae
handle possible error from EIG
John W. Eaton <jwe@octave.org>
parents:
7543
diff
changeset
|
109 if (! error_state) |
1 | 110 { |
7578
91f8446ce4ae
handle possible error from EIG
John W. Eaton <jwe@octave.org>
parents:
7543
diff
changeset
|
111 ComplexColumnVector lambda (b_eig.eigenvalues ()); |
91f8446ce4ae
handle possible error from EIG
John W. Eaton <jwe@octave.org>
parents:
7543
diff
changeset
|
112 ComplexMatrix Q (b_eig.eigenvectors ()); |
91f8446ce4ae
handle possible error from EIG
John W. Eaton <jwe@octave.org>
parents:
7543
diff
changeset
|
113 |
91f8446ce4ae
handle possible error from EIG
John W. Eaton <jwe@octave.org>
parents:
7543
diff
changeset
|
114 for (octave_idx_type i = 0; i < nr; i++) |
91f8446ce4ae
handle possible error from EIG
John W. Eaton <jwe@octave.org>
parents:
7543
diff
changeset
|
115 { |
91f8446ce4ae
handle possible error from EIG
John W. Eaton <jwe@octave.org>
parents:
7543
diff
changeset
|
116 Complex elt = lambda(i); |
91f8446ce4ae
handle possible error from EIG
John W. Eaton <jwe@octave.org>
parents:
7543
diff
changeset
|
117 if (std::imag (elt) == 0.0) |
91f8446ce4ae
handle possible error from EIG
John W. Eaton <jwe@octave.org>
parents:
7543
diff
changeset
|
118 lambda(i) = std::pow (a, std::real (elt)); |
91f8446ce4ae
handle possible error from EIG
John W. Eaton <jwe@octave.org>
parents:
7543
diff
changeset
|
119 else |
91f8446ce4ae
handle possible error from EIG
John W. Eaton <jwe@octave.org>
parents:
7543
diff
changeset
|
120 lambda(i) = std::pow (a, elt); |
91f8446ce4ae
handle possible error from EIG
John W. Eaton <jwe@octave.org>
parents:
7543
diff
changeset
|
121 } |
91f8446ce4ae
handle possible error from EIG
John W. Eaton <jwe@octave.org>
parents:
7543
diff
changeset
|
122 ComplexDiagMatrix D (lambda); |
91f8446ce4ae
handle possible error from EIG
John W. Eaton <jwe@octave.org>
parents:
7543
diff
changeset
|
123 |
91f8446ce4ae
handle possible error from EIG
John W. Eaton <jwe@octave.org>
parents:
7543
diff
changeset
|
124 retval = ComplexMatrix (Q * D * Q.inverse ()); |
1 | 125 } |
7578
91f8446ce4ae
handle possible error from EIG
John W. Eaton <jwe@octave.org>
parents:
7543
diff
changeset
|
126 else |
91f8446ce4ae
handle possible error from EIG
John W. Eaton <jwe@octave.org>
parents:
7543
diff
changeset
|
127 error ("xpow: matrix diagonalization failed"); |
1 | 128 } |
129 | |
130 return retval; | |
131 } | |
132 | |
767 | 133 // -*- 3 -*- |
2086 | 134 octave_value |
164 | 135 xpow (double a, const Complex& b) |
1 | 136 { |
137 Complex result; | |
138 Complex atmp (a); | |
5260 | 139 result = std::pow (atmp, b); |
1567 | 140 return result; |
1 | 141 } |
142 | |
767 | 143 // -*- 4 -*- |
2086 | 144 octave_value |
164 | 145 xpow (double a, const ComplexMatrix& b) |
1 | 146 { |
2086 | 147 octave_value retval; |
1 | 148 |
5275 | 149 octave_idx_type nr = b.rows (); |
150 octave_idx_type nc = b.cols (); | |
1 | 151 |
152 if (nr == 0 || nc == 0 || nr != nc) | |
2365 | 153 error ("for x^A, A must be square"); |
1 | 154 else |
155 { | |
156 EIG b_eig (b); | |
157 | |
7578
91f8446ce4ae
handle possible error from EIG
John W. Eaton <jwe@octave.org>
parents:
7543
diff
changeset
|
158 if (! error_state) |
1 | 159 { |
7578
91f8446ce4ae
handle possible error from EIG
John W. Eaton <jwe@octave.org>
parents:
7543
diff
changeset
|
160 ComplexColumnVector lambda (b_eig.eigenvalues ()); |
91f8446ce4ae
handle possible error from EIG
John W. Eaton <jwe@octave.org>
parents:
7543
diff
changeset
|
161 ComplexMatrix Q (b_eig.eigenvectors ()); |
91f8446ce4ae
handle possible error from EIG
John W. Eaton <jwe@octave.org>
parents:
7543
diff
changeset
|
162 |
91f8446ce4ae
handle possible error from EIG
John W. Eaton <jwe@octave.org>
parents:
7543
diff
changeset
|
163 for (octave_idx_type i = 0; i < nr; i++) |
91f8446ce4ae
handle possible error from EIG
John W. Eaton <jwe@octave.org>
parents:
7543
diff
changeset
|
164 { |
91f8446ce4ae
handle possible error from EIG
John W. Eaton <jwe@octave.org>
parents:
7543
diff
changeset
|
165 Complex elt = lambda(i); |
91f8446ce4ae
handle possible error from EIG
John W. Eaton <jwe@octave.org>
parents:
7543
diff
changeset
|
166 if (std::imag (elt) == 0.0) |
91f8446ce4ae
handle possible error from EIG
John W. Eaton <jwe@octave.org>
parents:
7543
diff
changeset
|
167 lambda(i) = std::pow (a, std::real (elt)); |
91f8446ce4ae
handle possible error from EIG
John W. Eaton <jwe@octave.org>
parents:
7543
diff
changeset
|
168 else |
91f8446ce4ae
handle possible error from EIG
John W. Eaton <jwe@octave.org>
parents:
7543
diff
changeset
|
169 lambda(i) = std::pow (a, elt); |
91f8446ce4ae
handle possible error from EIG
John W. Eaton <jwe@octave.org>
parents:
7543
diff
changeset
|
170 } |
91f8446ce4ae
handle possible error from EIG
John W. Eaton <jwe@octave.org>
parents:
7543
diff
changeset
|
171 ComplexDiagMatrix D (lambda); |
91f8446ce4ae
handle possible error from EIG
John W. Eaton <jwe@octave.org>
parents:
7543
diff
changeset
|
172 |
91f8446ce4ae
handle possible error from EIG
John W. Eaton <jwe@octave.org>
parents:
7543
diff
changeset
|
173 retval = ComplexMatrix (Q * D * Q.inverse ()); |
1 | 174 } |
7578
91f8446ce4ae
handle possible error from EIG
John W. Eaton <jwe@octave.org>
parents:
7543
diff
changeset
|
175 else |
91f8446ce4ae
handle possible error from EIG
John W. Eaton <jwe@octave.org>
parents:
7543
diff
changeset
|
176 error ("xpow: matrix diagonalization failed"); |
1 | 177 } |
178 | |
179 return retval; | |
180 } | |
181 | |
767 | 182 // -*- 5 -*- |
2086 | 183 octave_value |
164 | 184 xpow (const Matrix& a, double b) |
1 | 185 { |
2086 | 186 octave_value retval; |
1 | 187 |
5275 | 188 octave_idx_type nr = a.rows (); |
189 octave_idx_type nc = a.cols (); | |
1 | 190 |
191 if (nr == 0 || nc == 0 || nr != nc) | |
2365 | 192 error ("for A^b, A must be square"); |
1567 | 193 else |
1 | 194 { |
2800 | 195 if (static_cast<int> (b) == b) |
1 | 196 { |
2804 | 197 int btmp = static_cast<int> (b); |
1567 | 198 if (btmp == 0) |
199 { | |
200 retval = DiagMatrix (nr, nr, 1.0); | |
201 } | |
202 else | |
203 { | |
204 // Too much copying? | |
5775 | 205 // FIXME -- we shouldn't do this if the exponent is |
1567 | 206 // large... |
207 | |
208 Matrix atmp; | |
209 if (btmp < 0) | |
210 { | |
211 btmp = -btmp; | |
1655 | 212 |
5275 | 213 octave_idx_type info; |
1655 | 214 double rcond = 0.0; |
6207 | 215 MatrixType mattype (a); |
1655 | 216 |
6207 | 217 atmp = a.inverse (mattype, info, rcond, 1); |
1655 | 218 |
219 if (info == -1) | |
220 warning ("inverse: matrix singular to machine\ | |
221 precision, rcond = %g", rcond); | |
1567 | 222 } |
223 else | |
224 atmp = a; | |
225 | |
226 Matrix result (atmp); | |
3178 | 227 |
228 btmp--; | |
229 | |
230 while (btmp > 0) | |
231 { | |
232 if (btmp & 1) | |
233 result = result * atmp; | |
234 | |
235 btmp >>= 1; | |
236 | |
237 if (btmp > 0) | |
238 atmp = atmp * atmp; | |
239 } | |
1567 | 240 |
241 retval = result; | |
242 } | |
1 | 243 } |
244 else | |
245 { | |
1567 | 246 EIG a_eig (a); |
7578
91f8446ce4ae
handle possible error from EIG
John W. Eaton <jwe@octave.org>
parents:
7543
diff
changeset
|
247 |
91f8446ce4ae
handle possible error from EIG
John W. Eaton <jwe@octave.org>
parents:
7543
diff
changeset
|
248 if (! error_state) |
91f8446ce4ae
handle possible error from EIG
John W. Eaton <jwe@octave.org>
parents:
7543
diff
changeset
|
249 { |
91f8446ce4ae
handle possible error from EIG
John W. Eaton <jwe@octave.org>
parents:
7543
diff
changeset
|
250 ComplexColumnVector lambda (a_eig.eigenvalues ()); |
91f8446ce4ae
handle possible error from EIG
John W. Eaton <jwe@octave.org>
parents:
7543
diff
changeset
|
251 ComplexMatrix Q (a_eig.eigenvectors ()); |
1567 | 252 |
7578
91f8446ce4ae
handle possible error from EIG
John W. Eaton <jwe@octave.org>
parents:
7543
diff
changeset
|
253 for (octave_idx_type i = 0; i < nr; i++) |
91f8446ce4ae
handle possible error from EIG
John W. Eaton <jwe@octave.org>
parents:
7543
diff
changeset
|
254 lambda(i) = std::pow (lambda(i), b); |
91f8446ce4ae
handle possible error from EIG
John W. Eaton <jwe@octave.org>
parents:
7543
diff
changeset
|
255 |
91f8446ce4ae
handle possible error from EIG
John W. Eaton <jwe@octave.org>
parents:
7543
diff
changeset
|
256 ComplexDiagMatrix D (lambda); |
1567 | 257 |
7578
91f8446ce4ae
handle possible error from EIG
John W. Eaton <jwe@octave.org>
parents:
7543
diff
changeset
|
258 retval = ComplexMatrix (Q * D * Q.inverse ()); |
91f8446ce4ae
handle possible error from EIG
John W. Eaton <jwe@octave.org>
parents:
7543
diff
changeset
|
259 } |
91f8446ce4ae
handle possible error from EIG
John W. Eaton <jwe@octave.org>
parents:
7543
diff
changeset
|
260 else |
91f8446ce4ae
handle possible error from EIG
John W. Eaton <jwe@octave.org>
parents:
7543
diff
changeset
|
261 error ("xpow: matrix diagonalization failed"); |
1567 | 262 } |
263 } | |
1358 | 264 |
1567 | 265 return retval; |
266 } | |
1 | 267 |
8382
9b20a4847056
implement scalar powers of diag matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
268 // -*- 5d -*- |
9b20a4847056
implement scalar powers of diag matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
269 octave_value |
9b20a4847056
implement scalar powers of diag matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
270 xpow (const DiagMatrix& a, double b) |
9b20a4847056
implement scalar powers of diag matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
271 { |
9b20a4847056
implement scalar powers of diag matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
272 octave_value retval; |
9b20a4847056
implement scalar powers of diag matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
273 |
9b20a4847056
implement scalar powers of diag matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
274 octave_idx_type nr = a.rows (); |
9b20a4847056
implement scalar powers of diag matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
275 octave_idx_type nc = a.cols (); |
9b20a4847056
implement scalar powers of diag matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
276 |
9b20a4847056
implement scalar powers of diag matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
277 if (nr == 0 || nc == 0 || nr != nc) |
9b20a4847056
implement scalar powers of diag matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
278 error ("for A^b, A must be square"); |
9b20a4847056
implement scalar powers of diag matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
279 else |
9b20a4847056
implement scalar powers of diag matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
280 { |
9b20a4847056
implement scalar powers of diag matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
281 if (static_cast<int> (b) == b) |
9b20a4847056
implement scalar powers of diag matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
282 { |
9b20a4847056
implement scalar powers of diag matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
283 DiagMatrix r (nr, nc); |
9b20a4847056
implement scalar powers of diag matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
284 for (octave_idx_type i = 0; i < nc; i++) |
9b20a4847056
implement scalar powers of diag matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
285 r(i, i) = std::pow (a(i, i), b); |
9b20a4847056
implement scalar powers of diag matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
286 retval = r; |
9b20a4847056
implement scalar powers of diag matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
287 } |
9b20a4847056
implement scalar powers of diag matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
288 else |
9b20a4847056
implement scalar powers of diag matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
289 { |
9b20a4847056
implement scalar powers of diag matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
290 ComplexDiagMatrix r (nr, nc); |
9b20a4847056
implement scalar powers of diag matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
291 for (octave_idx_type i = 0; i < nc; i++) |
9b20a4847056
implement scalar powers of diag matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
292 r(i, i) = std::pow (static_cast<Complex> (a(i, i)), b); |
9b20a4847056
implement scalar powers of diag matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
293 retval = r; |
9b20a4847056
implement scalar powers of diag matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
294 } |
9b20a4847056
implement scalar powers of diag matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
295 } |
9b20a4847056
implement scalar powers of diag matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
296 |
9b20a4847056
implement scalar powers of diag matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
297 return retval; |
9b20a4847056
implement scalar powers of diag matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
298 } |
9b20a4847056
implement scalar powers of diag matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
299 |
8958
6ccc12cc65ef
implement raising a permutation matrix to integer power
Jaroslav Hajek <highegg@gmail.com>
parents:
8920
diff
changeset
|
300 // -*- 5p -*- |
6ccc12cc65ef
implement raising a permutation matrix to integer power
Jaroslav Hajek <highegg@gmail.com>
parents:
8920
diff
changeset
|
301 octave_value |
6ccc12cc65ef
implement raising a permutation matrix to integer power
Jaroslav Hajek <highegg@gmail.com>
parents:
8920
diff
changeset
|
302 xpow (const PermMatrix& a, double b) |
6ccc12cc65ef
implement raising a permutation matrix to integer power
Jaroslav Hajek <highegg@gmail.com>
parents:
8920
diff
changeset
|
303 { |
6ccc12cc65ef
implement raising a permutation matrix to integer power
Jaroslav Hajek <highegg@gmail.com>
parents:
8920
diff
changeset
|
304 octave_value retval; |
6ccc12cc65ef
implement raising a permutation matrix to integer power
Jaroslav Hajek <highegg@gmail.com>
parents:
8920
diff
changeset
|
305 int btmp = static_cast<int> (b); |
6ccc12cc65ef
implement raising a permutation matrix to integer power
Jaroslav Hajek <highegg@gmail.com>
parents:
8920
diff
changeset
|
306 if (btmp == b) |
6ccc12cc65ef
implement raising a permutation matrix to integer power
Jaroslav Hajek <highegg@gmail.com>
parents:
8920
diff
changeset
|
307 return a.power (btmp); |
6ccc12cc65ef
implement raising a permutation matrix to integer power
Jaroslav Hajek <highegg@gmail.com>
parents:
8920
diff
changeset
|
308 else |
6ccc12cc65ef
implement raising a permutation matrix to integer power
Jaroslav Hajek <highegg@gmail.com>
parents:
8920
diff
changeset
|
309 return xpow (Matrix (a), b); |
6ccc12cc65ef
implement raising a permutation matrix to integer power
Jaroslav Hajek <highegg@gmail.com>
parents:
8920
diff
changeset
|
310 } |
6ccc12cc65ef
implement raising a permutation matrix to integer power
Jaroslav Hajek <highegg@gmail.com>
parents:
8920
diff
changeset
|
311 |
1567 | 312 // -*- 6 -*- |
2086 | 313 octave_value |
1567 | 314 xpow (const Matrix& a, const Complex& b) |
315 { | |
2086 | 316 octave_value retval; |
1 | 317 |
5275 | 318 octave_idx_type nr = a.rows (); |
319 octave_idx_type nc = a.cols (); | |
1567 | 320 |
321 if (nr == 0 || nc == 0 || nr != nc) | |
2365 | 322 error ("for A^b, A must be square"); |
1 | 323 else |
324 { | |
325 EIG a_eig (a); | |
7578
91f8446ce4ae
handle possible error from EIG
John W. Eaton <jwe@octave.org>
parents:
7543
diff
changeset
|
326 |
91f8446ce4ae
handle possible error from EIG
John W. Eaton <jwe@octave.org>
parents:
7543
diff
changeset
|
327 if (! error_state) |
91f8446ce4ae
handle possible error from EIG
John W. Eaton <jwe@octave.org>
parents:
7543
diff
changeset
|
328 { |
91f8446ce4ae
handle possible error from EIG
John W. Eaton <jwe@octave.org>
parents:
7543
diff
changeset
|
329 ComplexColumnVector lambda (a_eig.eigenvalues ()); |
91f8446ce4ae
handle possible error from EIG
John W. Eaton <jwe@octave.org>
parents:
7543
diff
changeset
|
330 ComplexMatrix Q (a_eig.eigenvectors ()); |
1 | 331 |
7578
91f8446ce4ae
handle possible error from EIG
John W. Eaton <jwe@octave.org>
parents:
7543
diff
changeset
|
332 for (octave_idx_type i = 0; i < nr; i++) |
91f8446ce4ae
handle possible error from EIG
John W. Eaton <jwe@octave.org>
parents:
7543
diff
changeset
|
333 lambda(i) = std::pow (lambda(i), b); |
91f8446ce4ae
handle possible error from EIG
John W. Eaton <jwe@octave.org>
parents:
7543
diff
changeset
|
334 |
91f8446ce4ae
handle possible error from EIG
John W. Eaton <jwe@octave.org>
parents:
7543
diff
changeset
|
335 ComplexDiagMatrix D (lambda); |
1 | 336 |
7578
91f8446ce4ae
handle possible error from EIG
John W. Eaton <jwe@octave.org>
parents:
7543
diff
changeset
|
337 retval = ComplexMatrix (Q * D * Q.inverse ()); |
91f8446ce4ae
handle possible error from EIG
John W. Eaton <jwe@octave.org>
parents:
7543
diff
changeset
|
338 } |
91f8446ce4ae
handle possible error from EIG
John W. Eaton <jwe@octave.org>
parents:
7543
diff
changeset
|
339 else |
91f8446ce4ae
handle possible error from EIG
John W. Eaton <jwe@octave.org>
parents:
7543
diff
changeset
|
340 error ("xpow: matrix diagonalization failed"); |
1 | 341 } |
342 | |
343 return retval; | |
344 } | |
345 | |
767 | 346 // -*- 7 -*- |
2086 | 347 octave_value |
164 | 348 xpow (const Complex& a, double b) |
1 | 349 { |
350 Complex result; | |
1567 | 351 |
352 if (xisint (b)) | |
5260 | 353 result = std::pow (a, static_cast<int> (b)); |
1567 | 354 else |
5260 | 355 result = std::pow (a, b); |
1567 | 356 |
357 return result; | |
1 | 358 } |
359 | |
767 | 360 // -*- 8 -*- |
2086 | 361 octave_value |
164 | 362 xpow (const Complex& a, const Matrix& b) |
1 | 363 { |
2086 | 364 octave_value retval; |
1 | 365 |
5275 | 366 octave_idx_type nr = b.rows (); |
367 octave_idx_type nc = b.cols (); | |
1 | 368 |
369 if (nr == 0 || nc == 0 || nr != nc) | |
2365 | 370 error ("for x^A, A must be square"); |
1 | 371 else |
372 { | |
373 EIG b_eig (b); | |
374 | |
7578
91f8446ce4ae
handle possible error from EIG
John W. Eaton <jwe@octave.org>
parents:
7543
diff
changeset
|
375 if (! error_state) |
1 | 376 { |
7578
91f8446ce4ae
handle possible error from EIG
John W. Eaton <jwe@octave.org>
parents:
7543
diff
changeset
|
377 ComplexColumnVector lambda (b_eig.eigenvalues ()); |
91f8446ce4ae
handle possible error from EIG
John W. Eaton <jwe@octave.org>
parents:
7543
diff
changeset
|
378 ComplexMatrix Q (b_eig.eigenvectors ()); |
91f8446ce4ae
handle possible error from EIG
John W. Eaton <jwe@octave.org>
parents:
7543
diff
changeset
|
379 |
91f8446ce4ae
handle possible error from EIG
John W. Eaton <jwe@octave.org>
parents:
7543
diff
changeset
|
380 for (octave_idx_type i = 0; i < nr; i++) |
91f8446ce4ae
handle possible error from EIG
John W. Eaton <jwe@octave.org>
parents:
7543
diff
changeset
|
381 { |
91f8446ce4ae
handle possible error from EIG
John W. Eaton <jwe@octave.org>
parents:
7543
diff
changeset
|
382 Complex elt = lambda(i); |
91f8446ce4ae
handle possible error from EIG
John W. Eaton <jwe@octave.org>
parents:
7543
diff
changeset
|
383 if (std::imag (elt) == 0.0) |
91f8446ce4ae
handle possible error from EIG
John W. Eaton <jwe@octave.org>
parents:
7543
diff
changeset
|
384 lambda(i) = std::pow (a, std::real (elt)); |
91f8446ce4ae
handle possible error from EIG
John W. Eaton <jwe@octave.org>
parents:
7543
diff
changeset
|
385 else |
91f8446ce4ae
handle possible error from EIG
John W. Eaton <jwe@octave.org>
parents:
7543
diff
changeset
|
386 lambda(i) = std::pow (a, elt); |
91f8446ce4ae
handle possible error from EIG
John W. Eaton <jwe@octave.org>
parents:
7543
diff
changeset
|
387 } |
91f8446ce4ae
handle possible error from EIG
John W. Eaton <jwe@octave.org>
parents:
7543
diff
changeset
|
388 ComplexDiagMatrix D (lambda); |
91f8446ce4ae
handle possible error from EIG
John W. Eaton <jwe@octave.org>
parents:
7543
diff
changeset
|
389 |
91f8446ce4ae
handle possible error from EIG
John W. Eaton <jwe@octave.org>
parents:
7543
diff
changeset
|
390 retval = ComplexMatrix (Q * D * Q.inverse ()); |
1 | 391 } |
7578
91f8446ce4ae
handle possible error from EIG
John W. Eaton <jwe@octave.org>
parents:
7543
diff
changeset
|
392 else |
91f8446ce4ae
handle possible error from EIG
John W. Eaton <jwe@octave.org>
parents:
7543
diff
changeset
|
393 error ("xpow: matrix diagonalization failed"); |
1 | 394 } |
395 | |
396 return retval; | |
397 } | |
398 | |
767 | 399 // -*- 9 -*- |
2086 | 400 octave_value |
164 | 401 xpow (const Complex& a, const Complex& b) |
1 | 402 { |
403 Complex result; | |
5260 | 404 result = std::pow (a, b); |
1567 | 405 return result; |
1 | 406 } |
407 | |
767 | 408 // -*- 10 -*- |
2086 | 409 octave_value |
164 | 410 xpow (const Complex& a, const ComplexMatrix& b) |
1 | 411 { |
2086 | 412 octave_value retval; |
1 | 413 |
5275 | 414 octave_idx_type nr = b.rows (); |
415 octave_idx_type nc = b.cols (); | |
1 | 416 |
417 if (nr == 0 || nc == 0 || nr != nc) | |
2365 | 418 error ("for x^A, A must be square"); |
1 | 419 else |
420 { | |
421 EIG b_eig (b); | |
422 | |
7578
91f8446ce4ae
handle possible error from EIG
John W. Eaton <jwe@octave.org>
parents:
7543
diff
changeset
|
423 if (! error_state) |
1 | 424 { |
7578
91f8446ce4ae
handle possible error from EIG
John W. Eaton <jwe@octave.org>
parents:
7543
diff
changeset
|
425 ComplexColumnVector lambda (b_eig.eigenvalues ()); |
91f8446ce4ae
handle possible error from EIG
John W. Eaton <jwe@octave.org>
parents:
7543
diff
changeset
|
426 ComplexMatrix Q (b_eig.eigenvectors ()); |
91f8446ce4ae
handle possible error from EIG
John W. Eaton <jwe@octave.org>
parents:
7543
diff
changeset
|
427 |
91f8446ce4ae
handle possible error from EIG
John W. Eaton <jwe@octave.org>
parents:
7543
diff
changeset
|
428 for (octave_idx_type i = 0; i < nr; i++) |
91f8446ce4ae
handle possible error from EIG
John W. Eaton <jwe@octave.org>
parents:
7543
diff
changeset
|
429 { |
91f8446ce4ae
handle possible error from EIG
John W. Eaton <jwe@octave.org>
parents:
7543
diff
changeset
|
430 Complex elt = lambda(i); |
91f8446ce4ae
handle possible error from EIG
John W. Eaton <jwe@octave.org>
parents:
7543
diff
changeset
|
431 if (std::imag (elt) == 0.0) |
91f8446ce4ae
handle possible error from EIG
John W. Eaton <jwe@octave.org>
parents:
7543
diff
changeset
|
432 lambda(i) = std::pow (a, std::real (elt)); |
91f8446ce4ae
handle possible error from EIG
John W. Eaton <jwe@octave.org>
parents:
7543
diff
changeset
|
433 else |
91f8446ce4ae
handle possible error from EIG
John W. Eaton <jwe@octave.org>
parents:
7543
diff
changeset
|
434 lambda(i) = std::pow (a, elt); |
91f8446ce4ae
handle possible error from EIG
John W. Eaton <jwe@octave.org>
parents:
7543
diff
changeset
|
435 } |
91f8446ce4ae
handle possible error from EIG
John W. Eaton <jwe@octave.org>
parents:
7543
diff
changeset
|
436 ComplexDiagMatrix D (lambda); |
91f8446ce4ae
handle possible error from EIG
John W. Eaton <jwe@octave.org>
parents:
7543
diff
changeset
|
437 |
91f8446ce4ae
handle possible error from EIG
John W. Eaton <jwe@octave.org>
parents:
7543
diff
changeset
|
438 retval = ComplexMatrix (Q * D * Q.inverse ()); |
1 | 439 } |
7578
91f8446ce4ae
handle possible error from EIG
John W. Eaton <jwe@octave.org>
parents:
7543
diff
changeset
|
440 else |
91f8446ce4ae
handle possible error from EIG
John W. Eaton <jwe@octave.org>
parents:
7543
diff
changeset
|
441 error ("xpow: matrix diagonalization failed"); |
1 | 442 } |
443 | |
444 return retval; | |
445 } | |
446 | |
767 | 447 // -*- 11 -*- |
2086 | 448 octave_value |
164 | 449 xpow (const ComplexMatrix& a, double b) |
1 | 450 { |
2086 | 451 octave_value retval; |
1 | 452 |
5275 | 453 octave_idx_type nr = a.rows (); |
454 octave_idx_type nc = a.cols (); | |
1 | 455 |
456 if (nr == 0 || nc == 0 || nr != nc) | |
2365 | 457 error ("for A^b, A must be square"); |
1567 | 458 else |
1 | 459 { |
2800 | 460 if (static_cast<int> (b) == b) |
1 | 461 { |
2804 | 462 int btmp = static_cast<int> (b); |
1567 | 463 if (btmp == 0) |
464 { | |
465 retval = DiagMatrix (nr, nr, 1.0); | |
466 } | |
467 else | |
468 { | |
469 // Too much copying? | |
5775 | 470 // FIXME -- we shouldn't do this if the exponent is |
1567 | 471 // large... |
472 | |
473 ComplexMatrix atmp; | |
474 if (btmp < 0) | |
475 { | |
476 btmp = -btmp; | |
1655 | 477 |
5275 | 478 octave_idx_type info; |
1655 | 479 double rcond = 0.0; |
6207 | 480 MatrixType mattype (a); |
1655 | 481 |
6207 | 482 atmp = a.inverse (mattype, info, rcond, 1); |
1655 | 483 |
484 if (info == -1) | |
485 warning ("inverse: matrix singular to machine\ | |
486 precision, rcond = %g", rcond); | |
1567 | 487 } |
488 else | |
489 atmp = a; | |
490 | |
491 ComplexMatrix result (atmp); | |
3178 | 492 |
493 btmp--; | |
494 | |
495 while (btmp > 0) | |
496 { | |
497 if (btmp & 1) | |
498 result = result * atmp; | |
499 | |
500 btmp >>= 1; | |
501 | |
502 if (btmp > 0) | |
503 atmp = atmp * atmp; | |
504 } | |
1567 | 505 |
506 retval = result; | |
507 } | |
1 | 508 } |
509 else | |
510 { | |
1567 | 511 EIG a_eig (a); |
7578
91f8446ce4ae
handle possible error from EIG
John W. Eaton <jwe@octave.org>
parents:
7543
diff
changeset
|
512 |
91f8446ce4ae
handle possible error from EIG
John W. Eaton <jwe@octave.org>
parents:
7543
diff
changeset
|
513 if (! error_state) |
91f8446ce4ae
handle possible error from EIG
John W. Eaton <jwe@octave.org>
parents:
7543
diff
changeset
|
514 { |
91f8446ce4ae
handle possible error from EIG
John W. Eaton <jwe@octave.org>
parents:
7543
diff
changeset
|
515 ComplexColumnVector lambda (a_eig.eigenvalues ()); |
91f8446ce4ae
handle possible error from EIG
John W. Eaton <jwe@octave.org>
parents:
7543
diff
changeset
|
516 ComplexMatrix Q (a_eig.eigenvectors ()); |
1567 | 517 |
7578
91f8446ce4ae
handle possible error from EIG
John W. Eaton <jwe@octave.org>
parents:
7543
diff
changeset
|
518 for (octave_idx_type i = 0; i < nr; i++) |
91f8446ce4ae
handle possible error from EIG
John W. Eaton <jwe@octave.org>
parents:
7543
diff
changeset
|
519 lambda(i) = std::pow (lambda(i), b); |
91f8446ce4ae
handle possible error from EIG
John W. Eaton <jwe@octave.org>
parents:
7543
diff
changeset
|
520 |
91f8446ce4ae
handle possible error from EIG
John W. Eaton <jwe@octave.org>
parents:
7543
diff
changeset
|
521 ComplexDiagMatrix D (lambda); |
1567 | 522 |
7578
91f8446ce4ae
handle possible error from EIG
John W. Eaton <jwe@octave.org>
parents:
7543
diff
changeset
|
523 retval = ComplexMatrix (Q * D * Q.inverse ()); |
91f8446ce4ae
handle possible error from EIG
John W. Eaton <jwe@octave.org>
parents:
7543
diff
changeset
|
524 } |
91f8446ce4ae
handle possible error from EIG
John W. Eaton <jwe@octave.org>
parents:
7543
diff
changeset
|
525 else |
91f8446ce4ae
handle possible error from EIG
John W. Eaton <jwe@octave.org>
parents:
7543
diff
changeset
|
526 error ("xpow: matrix diagonalization failed"); |
1567 | 527 } |
528 } | |
1358 | 529 |
1567 | 530 return retval; |
531 } | |
1 | 532 |
1567 | 533 // -*- 12 -*- |
2086 | 534 octave_value |
1567 | 535 xpow (const ComplexMatrix& a, const Complex& b) |
536 { | |
2086 | 537 octave_value retval; |
1 | 538 |
5275 | 539 octave_idx_type nr = a.rows (); |
540 octave_idx_type nc = a.cols (); | |
1567 | 541 |
542 if (nr == 0 || nc == 0 || nr != nc) | |
2365 | 543 error ("for A^b, A must be square"); |
1 | 544 else |
545 { | |
546 EIG a_eig (a); | |
7578
91f8446ce4ae
handle possible error from EIG
John W. Eaton <jwe@octave.org>
parents:
7543
diff
changeset
|
547 |
91f8446ce4ae
handle possible error from EIG
John W. Eaton <jwe@octave.org>
parents:
7543
diff
changeset
|
548 if (! error_state) |
91f8446ce4ae
handle possible error from EIG
John W. Eaton <jwe@octave.org>
parents:
7543
diff
changeset
|
549 { |
91f8446ce4ae
handle possible error from EIG
John W. Eaton <jwe@octave.org>
parents:
7543
diff
changeset
|
550 ComplexColumnVector lambda (a_eig.eigenvalues ()); |
91f8446ce4ae
handle possible error from EIG
John W. Eaton <jwe@octave.org>
parents:
7543
diff
changeset
|
551 ComplexMatrix Q (a_eig.eigenvectors ()); |
1 | 552 |
7578
91f8446ce4ae
handle possible error from EIG
John W. Eaton <jwe@octave.org>
parents:
7543
diff
changeset
|
553 for (octave_idx_type i = 0; i < nr; i++) |
91f8446ce4ae
handle possible error from EIG
John W. Eaton <jwe@octave.org>
parents:
7543
diff
changeset
|
554 lambda(i) = std::pow (lambda(i), b); |
91f8446ce4ae
handle possible error from EIG
John W. Eaton <jwe@octave.org>
parents:
7543
diff
changeset
|
555 |
91f8446ce4ae
handle possible error from EIG
John W. Eaton <jwe@octave.org>
parents:
7543
diff
changeset
|
556 ComplexDiagMatrix D (lambda); |
1 | 557 |
7578
91f8446ce4ae
handle possible error from EIG
John W. Eaton <jwe@octave.org>
parents:
7543
diff
changeset
|
558 retval = ComplexMatrix (Q * D * Q.inverse ()); |
91f8446ce4ae
handle possible error from EIG
John W. Eaton <jwe@octave.org>
parents:
7543
diff
changeset
|
559 } |
91f8446ce4ae
handle possible error from EIG
John W. Eaton <jwe@octave.org>
parents:
7543
diff
changeset
|
560 else |
91f8446ce4ae
handle possible error from EIG
John W. Eaton <jwe@octave.org>
parents:
7543
diff
changeset
|
561 error ("xpow: matrix diagonalization failed"); |
1 | 562 } |
563 | |
564 return retval; | |
565 } | |
566 | |
8382
9b20a4847056
implement scalar powers of diag matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
567 // -*- 12d -*- |
9b20a4847056
implement scalar powers of diag matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
568 octave_value |
9b20a4847056
implement scalar powers of diag matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
569 xpow (const ComplexDiagMatrix& a, const Complex& b) |
9b20a4847056
implement scalar powers of diag matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
570 { |
9b20a4847056
implement scalar powers of diag matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
571 octave_value retval; |
9b20a4847056
implement scalar powers of diag matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
572 |
9b20a4847056
implement scalar powers of diag matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
573 octave_idx_type nr = a.rows (); |
9b20a4847056
implement scalar powers of diag matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
574 octave_idx_type nc = a.cols (); |
9b20a4847056
implement scalar powers of diag matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
575 |
9b20a4847056
implement scalar powers of diag matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
576 if (nr == 0 || nc == 0 || nr != nc) |
9b20a4847056
implement scalar powers of diag matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
577 error ("for A^b, A must be square"); |
9b20a4847056
implement scalar powers of diag matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
578 else |
9b20a4847056
implement scalar powers of diag matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
579 { |
9b20a4847056
implement scalar powers of diag matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
580 ComplexDiagMatrix r (nr, nc); |
9b20a4847056
implement scalar powers of diag matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
581 for (octave_idx_type i = 0; i < nc; i++) |
9b20a4847056
implement scalar powers of diag matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
582 r(i, i) = std::pow (a(i, i), b); |
9b20a4847056
implement scalar powers of diag matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
583 retval = r; |
9b20a4847056
implement scalar powers of diag matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
584 } |
9b20a4847056
implement scalar powers of diag matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
585 |
9b20a4847056
implement scalar powers of diag matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
586 return retval; |
9b20a4847056
implement scalar powers of diag matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
587 } |
9b20a4847056
implement scalar powers of diag matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
588 |
9b20a4847056
implement scalar powers of diag matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
589 // mixed |
9b20a4847056
implement scalar powers of diag matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
590 octave_value |
9b20a4847056
implement scalar powers of diag matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
591 xpow (const ComplexDiagMatrix& a, double b) |
9b20a4847056
implement scalar powers of diag matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
592 { |
9b20a4847056
implement scalar powers of diag matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
593 return xpow (a, static_cast<Complex> (b)); |
9b20a4847056
implement scalar powers of diag matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
594 } |
9b20a4847056
implement scalar powers of diag matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
595 |
9b20a4847056
implement scalar powers of diag matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
596 octave_value |
9b20a4847056
implement scalar powers of diag matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
597 xpow (const DiagMatrix& a, const Complex& b) |
9b20a4847056
implement scalar powers of diag matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
598 { |
9b20a4847056
implement scalar powers of diag matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
599 return xpow (ComplexDiagMatrix (a), b); |
9b20a4847056
implement scalar powers of diag matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
600 } |
9b20a4847056
implement scalar powers of diag matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
601 |
9b20a4847056
implement scalar powers of diag matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
602 |
767 | 603 // Safer pow functions that work elementwise for matrices. |
604 // | |
605 // op2 \ op1: s m cs cm | |
606 // +-- +---+---+----+----+ | |
607 // scalar | | * | 3 | * | 9 | | |
608 // +---+---+----+----+ | |
609 // matrix | 1 | 4 | 7 | 10 | | |
610 // +---+---+----+----+ | |
611 // complex_scalar | * | 5 | * | 11 | | |
612 // +---+---+----+----+ | |
613 // complex_matrix | 2 | 6 | 8 | 12 | | |
614 // +---+---+----+----+ | |
615 // | |
616 // * -> not needed. | |
1 | 617 |
5775 | 618 // FIXME -- these functions need to be fixed so that things |
3162 | 619 // like |
620 // | |
621 // a = -1; b = [ 0, 0.5, 1 ]; r = a .^ b | |
622 // | |
623 // and | |
624 // | |
625 // a = -1; b = [ 0, 0.5, 1 ]; for i = 1:3, r(i) = a .^ b(i), end | |
626 // | |
627 // produce identical results. Also, it would be nice if -1^0.5 | |
628 // produced a pure imaginary result instead of a complex number with a | |
629 // small real part. But perhaps that's really a problem with the math | |
630 // library... | |
631 | |
767 | 632 // -*- 1 -*- |
2086 | 633 octave_value |
164 | 634 elem_xpow (double a, const Matrix& b) |
1 | 635 { |
2086 | 636 octave_value retval; |
1 | 637 |
5275 | 638 octave_idx_type nr = b.rows (); |
639 octave_idx_type nc = b.cols (); | |
1 | 640 |
3162 | 641 double d1, d2; |
1358 | 642 |
3162 | 643 if (a < 0.0 && ! b.all_integers (d1, d2)) |
1 | 644 { |
645 Complex atmp (a); | |
646 ComplexMatrix result (nr, nc); | |
5275 | 647 |
648 for (octave_idx_type j = 0; j < nc; j++) | |
649 for (octave_idx_type i = 0; i < nr; i++) | |
4153 | 650 { |
651 OCTAVE_QUIT; | |
5260 | 652 result (i, j) = std::pow (atmp, b (i, j)); |
4153 | 653 } |
1 | 654 |
1567 | 655 retval = result; |
1 | 656 } |
657 else | |
658 { | |
659 Matrix result (nr, nc); | |
5275 | 660 |
661 for (octave_idx_type j = 0; j < nc; j++) | |
662 for (octave_idx_type i = 0; i < nr; i++) | |
4153 | 663 { |
664 OCTAVE_QUIT; | |
5260 | 665 result (i, j) = std::pow (a, b (i, j)); |
4153 | 666 } |
1 | 667 |
1567 | 668 retval = result; |
1 | 669 } |
670 | |
671 return retval; | |
672 } | |
673 | |
767 | 674 // -*- 2 -*- |
2086 | 675 octave_value |
164 | 676 elem_xpow (double a, const ComplexMatrix& b) |
1 | 677 { |
5275 | 678 octave_idx_type nr = b.rows (); |
679 octave_idx_type nc = b.cols (); | |
1 | 680 |
681 ComplexMatrix result (nr, nc); | |
3125 | 682 Complex atmp (a); |
5275 | 683 |
684 for (octave_idx_type j = 0; j < nc; j++) | |
685 for (octave_idx_type i = 0; i < nr; i++) | |
4153 | 686 { |
687 OCTAVE_QUIT; | |
5260 | 688 result (i, j) = std::pow (atmp, b (i, j)); |
4153 | 689 } |
1 | 690 |
1567 | 691 return result; |
1 | 692 } |
693 | |
767 | 694 // -*- 3 -*- |
2086 | 695 octave_value |
164 | 696 elem_xpow (const Matrix& a, double b) |
1 | 697 { |
2086 | 698 octave_value retval; |
1 | 699 |
5275 | 700 octave_idx_type nr = a.rows (); |
701 octave_idx_type nc = a.cols (); | |
1 | 702 |
8978
0a58c4cd1405
optimize element-wise power
Jaroslav Hajek <highegg@gmail.com>
parents:
8960
diff
changeset
|
703 if (! xisint (b) && a.any_element_is_negative ()) |
1 | 704 { |
705 ComplexMatrix result (nr, nc); | |
5275 | 706 |
707 for (octave_idx_type j = 0; j < nc; j++) | |
708 for (octave_idx_type i = 0; i < nr; i++) | |
1 | 709 { |
5665 | 710 OCTAVE_QUIT; |
711 | |
712 Complex atmp (a (i, j)); | |
713 | |
714 result (i, j) = std::pow (atmp, b); | |
1 | 715 } |
716 | |
1567 | 717 retval = result; |
1 | 718 } |
719 else | |
720 { | |
721 Matrix result (nr, nc); | |
5275 | 722 |
723 for (octave_idx_type j = 0; j < nc; j++) | |
724 for (octave_idx_type i = 0; i < nr; i++) | |
4153 | 725 { |
726 OCTAVE_QUIT; | |
5260 | 727 result (i, j) = std::pow (a (i, j), b); |
4153 | 728 } |
1 | 729 |
1567 | 730 retval = result; |
1 | 731 } |
732 | |
733 return retval; | |
734 } | |
735 | |
767 | 736 // -*- 4 -*- |
2086 | 737 octave_value |
164 | 738 elem_xpow (const Matrix& a, const Matrix& b) |
1 | 739 { |
2086 | 740 octave_value retval; |
1567 | 741 |
5275 | 742 octave_idx_type nr = a.rows (); |
743 octave_idx_type nc = a.cols (); | |
2365 | 744 |
5275 | 745 octave_idx_type b_nr = b.rows (); |
746 octave_idx_type b_nc = b.cols (); | |
1 | 747 |
2365 | 748 if (nr != b_nr || nc != b_nc) |
749 { | |
750 gripe_nonconformant ("operator .^", nr, nc, b_nr, b_nc); | |
751 return octave_value (); | |
752 } | |
1 | 753 |
754 int convert_to_complex = 0; | |
5275 | 755 for (octave_idx_type j = 0; j < nc; j++) |
756 for (octave_idx_type i = 0; i < nr; i++) | |
1 | 757 { |
4153 | 758 OCTAVE_QUIT; |
2305 | 759 double atmp = a (i, j); |
760 double btmp = b (i, j); | |
2800 | 761 if (atmp < 0.0 && static_cast<int> (btmp) != btmp) |
1 | 762 { |
763 convert_to_complex = 1; | |
764 goto done; | |
765 } | |
766 } | |
767 | |
2365 | 768 done: |
1 | 769 |
770 if (convert_to_complex) | |
771 { | |
772 ComplexMatrix complex_result (nr, nc); | |
773 | |
5275 | 774 for (octave_idx_type j = 0; j < nc; j++) |
775 for (octave_idx_type i = 0; i < nr; i++) | |
1 | 776 { |
4153 | 777 OCTAVE_QUIT; |
5665 | 778 Complex atmp (a (i, j)); |
779 Complex btmp (b (i, j)); | |
780 complex_result (i, j) = std::pow (atmp, btmp); | |
1 | 781 } |
1567 | 782 |
783 retval = complex_result; | |
1 | 784 } |
785 else | |
786 { | |
787 Matrix result (nr, nc); | |
788 | |
5275 | 789 for (octave_idx_type j = 0; j < nc; j++) |
790 for (octave_idx_type i = 0; i < nr; i++) | |
4153 | 791 { |
792 OCTAVE_QUIT; | |
5260 | 793 result (i, j) = std::pow (a (i, j), b (i, j)); |
4153 | 794 } |
1 | 795 |
1567 | 796 retval = result; |
1 | 797 } |
1567 | 798 |
799 return retval; | |
1 | 800 } |
801 | |
767 | 802 // -*- 5 -*- |
2086 | 803 octave_value |
164 | 804 elem_xpow (const Matrix& a, const Complex& b) |
1 | 805 { |
5275 | 806 octave_idx_type nr = a.rows (); |
807 octave_idx_type nc = a.cols (); | |
1 | 808 |
809 ComplexMatrix result (nr, nc); | |
5275 | 810 |
811 for (octave_idx_type j = 0; j < nc; j++) | |
812 for (octave_idx_type i = 0; i < nr; i++) | |
4153 | 813 { |
814 OCTAVE_QUIT; | |
5260 | 815 result (i, j) = std::pow (Complex (a (i, j)), b); |
4153 | 816 } |
1 | 817 |
1567 | 818 return result; |
1 | 819 } |
820 | |
767 | 821 // -*- 6 -*- |
2086 | 822 octave_value |
164 | 823 elem_xpow (const Matrix& a, const ComplexMatrix& b) |
1 | 824 { |
5275 | 825 octave_idx_type nr = a.rows (); |
826 octave_idx_type nc = a.cols (); | |
2365 | 827 |
5275 | 828 octave_idx_type b_nr = b.rows (); |
829 octave_idx_type b_nc = b.cols (); | |
1 | 830 |
2365 | 831 if (nr != b_nr || nc != b_nc) |
832 { | |
833 gripe_nonconformant ("operator .^", nr, nc, b_nr, b_nc); | |
834 return octave_value (); | |
835 } | |
1 | 836 |
837 ComplexMatrix result (nr, nc); | |
5275 | 838 |
839 for (octave_idx_type j = 0; j < nc; j++) | |
840 for (octave_idx_type i = 0; i < nr; i++) | |
4153 | 841 { |
842 OCTAVE_QUIT; | |
5260 | 843 result (i, j) = std::pow (Complex (a (i, j)), b (i, j)); |
4153 | 844 } |
1 | 845 |
1567 | 846 return result; |
1 | 847 } |
848 | |
767 | 849 // -*- 7 -*- |
2086 | 850 octave_value |
164 | 851 elem_xpow (const Complex& a, const Matrix& b) |
1 | 852 { |
5275 | 853 octave_idx_type nr = b.rows (); |
854 octave_idx_type nc = b.cols (); | |
1 | 855 |
856 ComplexMatrix result (nr, nc); | |
5275 | 857 |
858 for (octave_idx_type j = 0; j < nc; j++) | |
859 for (octave_idx_type i = 0; i < nr; i++) | |
1567 | 860 { |
4153 | 861 OCTAVE_QUIT; |
2305 | 862 double btmp = b (i, j); |
1567 | 863 if (xisint (btmp)) |
5260 | 864 result (i, j) = std::pow (a, static_cast<int> (btmp)); |
1567 | 865 else |
5260 | 866 result (i, j) = std::pow (a, btmp); |
1567 | 867 } |
1 | 868 |
1567 | 869 return result; |
1 | 870 } |
871 | |
767 | 872 // -*- 8 -*- |
2086 | 873 octave_value |
164 | 874 elem_xpow (const Complex& a, const ComplexMatrix& b) |
1 | 875 { |
5275 | 876 octave_idx_type nr = b.rows (); |
877 octave_idx_type nc = b.cols (); | |
1 | 878 |
879 ComplexMatrix result (nr, nc); | |
5275 | 880 |
881 for (octave_idx_type j = 0; j < nc; j++) | |
882 for (octave_idx_type i = 0; i < nr; i++) | |
4153 | 883 { |
884 OCTAVE_QUIT; | |
5260 | 885 result (i, j) = std::pow (a, b (i, j)); |
4153 | 886 } |
1 | 887 |
1567 | 888 return result; |
1 | 889 } |
890 | |
767 | 891 // -*- 9 -*- |
2086 | 892 octave_value |
164 | 893 elem_xpow (const ComplexMatrix& a, double b) |
1 | 894 { |
5275 | 895 octave_idx_type nr = a.rows (); |
896 octave_idx_type nc = a.cols (); | |
1 | 897 |
898 ComplexMatrix result (nr, nc); | |
899 | |
1567 | 900 if (xisint (b)) |
901 { | |
5275 | 902 for (octave_idx_type j = 0; j < nc; j++) |
903 for (octave_idx_type i = 0; i < nr; i++) | |
4153 | 904 { |
905 OCTAVE_QUIT; | |
5260 | 906 result (i, j) = std::pow (a (i, j), static_cast<int> (b)); |
4153 | 907 } |
1567 | 908 } |
909 else | |
910 { | |
5275 | 911 for (octave_idx_type j = 0; j < nc; j++) |
912 for (octave_idx_type i = 0; i < nr; i++) | |
4153 | 913 { |
914 OCTAVE_QUIT; | |
5260 | 915 result (i, j) = std::pow (a (i, j), b); |
4153 | 916 } |
1567 | 917 } |
918 | |
919 return result; | |
1 | 920 } |
921 | |
767 | 922 // -*- 10 -*- |
2086 | 923 octave_value |
164 | 924 elem_xpow (const ComplexMatrix& a, const Matrix& b) |
1 | 925 { |
5275 | 926 octave_idx_type nr = a.rows (); |
927 octave_idx_type nc = a.cols (); | |
2365 | 928 |
5275 | 929 octave_idx_type b_nr = b.rows (); |
930 octave_idx_type b_nc = b.cols (); | |
1 | 931 |
2365 | 932 if (nr != b_nr || nc != b_nc) |
933 { | |
934 gripe_nonconformant ("operator .^", nr, nc, b_nr, b_nc); | |
935 return octave_value (); | |
936 } | |
1 | 937 |
938 ComplexMatrix result (nr, nc); | |
5275 | 939 |
940 for (octave_idx_type j = 0; j < nc; j++) | |
941 for (octave_idx_type i = 0; i < nr; i++) | |
1567 | 942 { |
4153 | 943 OCTAVE_QUIT; |
2305 | 944 double btmp = b (i, j); |
1567 | 945 if (xisint (btmp)) |
5260 | 946 result (i, j) = std::pow (a (i, j), static_cast<int> (btmp)); |
1567 | 947 else |
5260 | 948 result (i, j) = std::pow (a (i, j), btmp); |
1567 | 949 } |
1 | 950 |
1567 | 951 return result; |
1 | 952 } |
953 | |
767 | 954 // -*- 11 -*- |
2086 | 955 octave_value |
164 | 956 elem_xpow (const ComplexMatrix& a, const Complex& b) |
1 | 957 { |
5275 | 958 octave_idx_type nr = a.rows (); |
959 octave_idx_type nc = a.cols (); | |
1 | 960 |
961 ComplexMatrix result (nr, nc); | |
5275 | 962 |
963 for (octave_idx_type j = 0; j < nc; j++) | |
964 for (octave_idx_type i = 0; i < nr; i++) | |
4153 | 965 { |
966 OCTAVE_QUIT; | |
5260 | 967 result (i, j) = std::pow (a (i, j), b); |
4153 | 968 } |
1 | 969 |
1567 | 970 return result; |
1 | 971 } |
972 | |
767 | 973 // -*- 12 -*- |
2086 | 974 octave_value |
164 | 975 elem_xpow (const ComplexMatrix& a, const ComplexMatrix& b) |
1 | 976 { |
5275 | 977 octave_idx_type nr = a.rows (); |
978 octave_idx_type nc = a.cols (); | |
2365 | 979 |
5275 | 980 octave_idx_type b_nr = b.rows (); |
981 octave_idx_type b_nc = b.cols (); | |
2365 | 982 |
983 if (nr != b_nr || nc != b_nc) | |
984 { | |
985 gripe_nonconformant ("operator .^", nr, nc, b_nr, b_nc); | |
986 return octave_value (); | |
987 } | |
1 | 988 |
989 ComplexMatrix result (nr, nc); | |
5275 | 990 |
991 for (octave_idx_type j = 0; j < nc; j++) | |
992 for (octave_idx_type i = 0; i < nr; i++) | |
4153 | 993 { |
994 OCTAVE_QUIT; | |
5260 | 995 result (i, j) = std::pow (a (i, j), b (i, j)); |
4153 | 996 } |
1 | 997 |
1567 | 998 return result; |
1 | 999 } |
1000 | |
4543 | 1001 // Safer pow functions that work elementwise for N-d arrays. |
1002 // | |
1003 // op2 \ op1: s nd cs cnd | |
1004 // +-- +---+---+----+----+ | |
1005 // scalar | | * | 3 | * | 9 | | |
1006 // +---+---+----+----+ | |
1007 // N_d | 1 | 4 | 7 | 10 | | |
1008 // +---+---+----+----+ | |
1009 // complex_scalar | * | 5 | * | 11 | | |
1010 // +---+---+----+----+ | |
1011 // complex_N_d | 2 | 6 | 8 | 12 | | |
1012 // +---+---+----+----+ | |
1013 // | |
1014 // * -> not needed. | |
1015 | |
5775 | 1016 // FIXME -- these functions need to be fixed so that things |
4543 | 1017 // like |
1018 // | |
1019 // a = -1; b = [ 0, 0.5, 1 ]; r = a .^ b | |
1020 // | |
1021 // and | |
1022 // | |
1023 // a = -1; b = [ 0, 0.5, 1 ]; for i = 1:3, r(i) = a .^ b(i), end | |
1024 // | |
1025 // produce identical results. Also, it would be nice if -1^0.5 | |
1026 // produced a pure imaginary result instead of a complex number with a | |
1027 // small real part. But perhaps that's really a problem with the math | |
1028 // library... | |
1029 | |
1030 // -*- 1 -*- | |
1031 octave_value | |
1032 elem_xpow (double a, const NDArray& b) | |
1033 { | |
1034 octave_value retval; | |
1035 | |
1036 double d1, d2; | |
1037 | |
1038 if (a < 0.0 && ! b.all_integers (d1, d2)) | |
1039 { | |
1040 Complex atmp (a); | |
1041 ComplexNDArray result (b.dims ()); | |
5275 | 1042 for (octave_idx_type i = 0; i < b.length (); i++) |
4543 | 1043 { |
1044 OCTAVE_QUIT; | |
5260 | 1045 result(i) = std::pow (atmp, b(i)); |
4543 | 1046 } |
1047 | |
1048 retval = result; | |
1049 } | |
1050 else | |
1051 { | |
1052 NDArray result (b.dims ()); | |
5275 | 1053 for (octave_idx_type i = 0; i < b.length (); i++) |
4543 | 1054 { |
1055 OCTAVE_QUIT; | |
5260 | 1056 result (i) = std::pow (a, b(i)); |
4543 | 1057 } |
1058 | |
1059 retval = result; | |
1060 } | |
1061 | |
1062 return retval; | |
1063 } | |
1064 | |
1065 // -*- 2 -*- | |
1066 octave_value | |
1067 elem_xpow (double a, const ComplexNDArray& b) | |
1068 { | |
1069 ComplexNDArray result (b.dims ()); | |
1070 Complex atmp (a); | |
5275 | 1071 |
1072 for (octave_idx_type i = 0; i < b.length (); i++) | |
4543 | 1073 { |
1074 OCTAVE_QUIT; | |
5260 | 1075 result(i) = std::pow (atmp, b(i)); |
4543 | 1076 } |
1077 | |
1078 return result; | |
1079 } | |
1080 | |
1081 // -*- 3 -*- | |
1082 octave_value | |
1083 elem_xpow (const NDArray& a, double b) | |
1084 { | |
1085 octave_value retval; | |
1086 | |
8978
0a58c4cd1405
optimize element-wise power
Jaroslav Hajek <highegg@gmail.com>
parents:
8960
diff
changeset
|
1087 if (! xisint (b)) |
4543 | 1088 { |
8978
0a58c4cd1405
optimize element-wise power
Jaroslav Hajek <highegg@gmail.com>
parents:
8960
diff
changeset
|
1089 if (a.any_element_is_negative ()) |
0a58c4cd1405
optimize element-wise power
Jaroslav Hajek <highegg@gmail.com>
parents:
8960
diff
changeset
|
1090 { |
0a58c4cd1405
optimize element-wise power
Jaroslav Hajek <highegg@gmail.com>
parents:
8960
diff
changeset
|
1091 ComplexNDArray result (a.dims ()); |
0a58c4cd1405
optimize element-wise power
Jaroslav Hajek <highegg@gmail.com>
parents:
8960
diff
changeset
|
1092 |
0a58c4cd1405
optimize element-wise power
Jaroslav Hajek <highegg@gmail.com>
parents:
8960
diff
changeset
|
1093 for (octave_idx_type i = 0; i < a.length (); i++) |
0a58c4cd1405
optimize element-wise power
Jaroslav Hajek <highegg@gmail.com>
parents:
8960
diff
changeset
|
1094 { |
0a58c4cd1405
optimize element-wise power
Jaroslav Hajek <highegg@gmail.com>
parents:
8960
diff
changeset
|
1095 OCTAVE_QUIT; |
0a58c4cd1405
optimize element-wise power
Jaroslav Hajek <highegg@gmail.com>
parents:
8960
diff
changeset
|
1096 |
0a58c4cd1405
optimize element-wise power
Jaroslav Hajek <highegg@gmail.com>
parents:
8960
diff
changeset
|
1097 Complex atmp (a (i)); |
0a58c4cd1405
optimize element-wise power
Jaroslav Hajek <highegg@gmail.com>
parents:
8960
diff
changeset
|
1098 |
0a58c4cd1405
optimize element-wise power
Jaroslav Hajek <highegg@gmail.com>
parents:
8960
diff
changeset
|
1099 result(i) = std::pow (atmp, b); |
0a58c4cd1405
optimize element-wise power
Jaroslav Hajek <highegg@gmail.com>
parents:
8960
diff
changeset
|
1100 } |
0a58c4cd1405
optimize element-wise power
Jaroslav Hajek <highegg@gmail.com>
parents:
8960
diff
changeset
|
1101 |
0a58c4cd1405
optimize element-wise power
Jaroslav Hajek <highegg@gmail.com>
parents:
8960
diff
changeset
|
1102 retval = result; |
0a58c4cd1405
optimize element-wise power
Jaroslav Hajek <highegg@gmail.com>
parents:
8960
diff
changeset
|
1103 } |
0a58c4cd1405
optimize element-wise power
Jaroslav Hajek <highegg@gmail.com>
parents:
8960
diff
changeset
|
1104 else |
0a58c4cd1405
optimize element-wise power
Jaroslav Hajek <highegg@gmail.com>
parents:
8960
diff
changeset
|
1105 { |
0a58c4cd1405
optimize element-wise power
Jaroslav Hajek <highegg@gmail.com>
parents:
8960
diff
changeset
|
1106 NDArray result (a.dims ()); |
0a58c4cd1405
optimize element-wise power
Jaroslav Hajek <highegg@gmail.com>
parents:
8960
diff
changeset
|
1107 for (octave_idx_type i = 0; i < a.length (); i++) |
0a58c4cd1405
optimize element-wise power
Jaroslav Hajek <highegg@gmail.com>
parents:
8960
diff
changeset
|
1108 { |
0a58c4cd1405
optimize element-wise power
Jaroslav Hajek <highegg@gmail.com>
parents:
8960
diff
changeset
|
1109 OCTAVE_QUIT; |
0a58c4cd1405
optimize element-wise power
Jaroslav Hajek <highegg@gmail.com>
parents:
8960
diff
changeset
|
1110 result(i) = std::pow (a(i), b); |
0a58c4cd1405
optimize element-wise power
Jaroslav Hajek <highegg@gmail.com>
parents:
8960
diff
changeset
|
1111 } |
0a58c4cd1405
optimize element-wise power
Jaroslav Hajek <highegg@gmail.com>
parents:
8960
diff
changeset
|
1112 |
0a58c4cd1405
optimize element-wise power
Jaroslav Hajek <highegg@gmail.com>
parents:
8960
diff
changeset
|
1113 retval = result; |
0a58c4cd1405
optimize element-wise power
Jaroslav Hajek <highegg@gmail.com>
parents:
8960
diff
changeset
|
1114 } |
4543 | 1115 } |
1116 else | |
1117 { | |
1118 NDArray result (a.dims ()); | |
1119 | |
8978
0a58c4cd1405
optimize element-wise power
Jaroslav Hajek <highegg@gmail.com>
parents:
8960
diff
changeset
|
1120 int ib = static_cast<int> (b); |
0a58c4cd1405
optimize element-wise power
Jaroslav Hajek <highegg@gmail.com>
parents:
8960
diff
changeset
|
1121 if (ib == 2) |
0a58c4cd1405
optimize element-wise power
Jaroslav Hajek <highegg@gmail.com>
parents:
8960
diff
changeset
|
1122 { |
0a58c4cd1405
optimize element-wise power
Jaroslav Hajek <highegg@gmail.com>
parents:
8960
diff
changeset
|
1123 for (octave_idx_type i = 0; i < a.length (); i++) |
0a58c4cd1405
optimize element-wise power
Jaroslav Hajek <highegg@gmail.com>
parents:
8960
diff
changeset
|
1124 result.xelem (i) = a(i) * a(i); |
0a58c4cd1405
optimize element-wise power
Jaroslav Hajek <highegg@gmail.com>
parents:
8960
diff
changeset
|
1125 } |
0a58c4cd1405
optimize element-wise power
Jaroslav Hajek <highegg@gmail.com>
parents:
8960
diff
changeset
|
1126 else |
0a58c4cd1405
optimize element-wise power
Jaroslav Hajek <highegg@gmail.com>
parents:
8960
diff
changeset
|
1127 { |
0a58c4cd1405
optimize element-wise power
Jaroslav Hajek <highegg@gmail.com>
parents:
8960
diff
changeset
|
1128 for (octave_idx_type i = 0; i < a.length (); i++) |
0a58c4cd1405
optimize element-wise power
Jaroslav Hajek <highegg@gmail.com>
parents:
8960
diff
changeset
|
1129 { |
0a58c4cd1405
optimize element-wise power
Jaroslav Hajek <highegg@gmail.com>
parents:
8960
diff
changeset
|
1130 OCTAVE_QUIT; |
0a58c4cd1405
optimize element-wise power
Jaroslav Hajek <highegg@gmail.com>
parents:
8960
diff
changeset
|
1131 result(i) = std::pow (a(i), ib); |
0a58c4cd1405
optimize element-wise power
Jaroslav Hajek <highegg@gmail.com>
parents:
8960
diff
changeset
|
1132 } |
0a58c4cd1405
optimize element-wise power
Jaroslav Hajek <highegg@gmail.com>
parents:
8960
diff
changeset
|
1133 } |
4543 | 1134 |
1135 retval = result; | |
1136 } | |
1137 | |
1138 return retval; | |
1139 } | |
1140 | |
1141 // -*- 4 -*- | |
1142 octave_value | |
1143 elem_xpow (const NDArray& a, const NDArray& b) | |
1144 { | |
1145 octave_value retval; | |
1146 | |
1147 dim_vector a_dims = a.dims (); | |
1148 dim_vector b_dims = b.dims (); | |
1149 | |
1150 if (a_dims != b_dims) | |
1151 { | |
1152 gripe_nonconformant ("operator .^", a_dims, b_dims); | |
1153 return octave_value (); | |
1154 } | |
1155 | |
1156 int len = a.length (); | |
1157 | |
1158 bool convert_to_complex = false; | |
1159 | |
5275 | 1160 for (octave_idx_type i = 0; i < len; i++) |
4543 | 1161 { |
1162 OCTAVE_QUIT; | |
1163 double atmp = a(i); | |
1164 double btmp = b(i); | |
1165 if (atmp < 0.0 && static_cast<int> (btmp) != btmp) | |
1166 { | |
1167 convert_to_complex = true; | |
1168 goto done; | |
1169 } | |
1170 } | |
1171 | |
1172 done: | |
1173 | |
1174 if (convert_to_complex) | |
1175 { | |
1176 ComplexNDArray complex_result (a_dims); | |
1177 | |
5275 | 1178 for (octave_idx_type i = 0; i < len; i++) |
4543 | 1179 { |
1180 OCTAVE_QUIT; | |
5665 | 1181 Complex atmp (a(i)); |
1182 Complex btmp (b(i)); | |
1183 complex_result(i) = std::pow (atmp, btmp); | |
4543 | 1184 } |
1185 | |
1186 retval = complex_result; | |
1187 } | |
1188 else | |
1189 { | |
1190 NDArray result (a_dims); | |
1191 | |
5275 | 1192 for (octave_idx_type i = 0; i < len; i++) |
4543 | 1193 { |
1194 OCTAVE_QUIT; | |
5260 | 1195 result(i) = std::pow (a(i), b(i)); |
4543 | 1196 } |
1197 | |
1198 retval = result; | |
1199 } | |
1200 | |
1201 return retval; | |
1202 } | |
1203 | |
1204 // -*- 5 -*- | |
1205 octave_value | |
1206 elem_xpow (const NDArray& a, const Complex& b) | |
1207 { | |
1208 ComplexNDArray result (a.dims ()); | |
1209 | |
5275 | 1210 for (octave_idx_type i = 0; i < a.length (); i++) |
4543 | 1211 { |
1212 OCTAVE_QUIT; | |
5260 | 1213 result(i) = std::pow (Complex (a(i)), b); |
4543 | 1214 } |
1215 | |
1216 return result; | |
1217 } | |
1218 | |
1219 // -*- 6 -*- | |
1220 octave_value | |
1221 elem_xpow (const NDArray& a, const ComplexNDArray& b) | |
1222 { | |
1223 dim_vector a_dims = a.dims (); | |
1224 dim_vector b_dims = b.dims (); | |
1225 | |
1226 if (a_dims != b_dims) | |
1227 { | |
1228 gripe_nonconformant ("operator .^", a_dims, b_dims); | |
1229 return octave_value (); | |
1230 } | |
1231 | |
1232 ComplexNDArray result (a_dims); | |
5275 | 1233 |
1234 for (octave_idx_type i = 0; i < a.length (); i++) | |
4543 | 1235 { |
1236 OCTAVE_QUIT; | |
5260 | 1237 result(i) = std::pow (Complex (a(i)), b(i)); |
4543 | 1238 } |
1239 | |
1240 return result; | |
1241 } | |
1242 | |
1243 // -*- 7 -*- | |
1244 octave_value | |
1245 elem_xpow (const Complex& a, const NDArray& b) | |
1246 { | |
1247 ComplexNDArray result (b.dims ()); | |
5275 | 1248 |
1249 for (octave_idx_type i = 0; i < b.length (); i++) | |
4543 | 1250 { |
1251 OCTAVE_QUIT; | |
1252 double btmp = b(i); | |
1253 if (xisint (btmp)) | |
5260 | 1254 result(i) = std::pow (a, static_cast<int> (btmp)); |
4543 | 1255 else |
5260 | 1256 result(i) = std::pow (a, btmp); |
4543 | 1257 } |
1258 | |
1259 return result; | |
1260 } | |
1261 | |
1262 // -*- 8 -*- | |
1263 octave_value | |
1264 elem_xpow (const Complex& a, const ComplexNDArray& b) | |
1265 { | |
1266 ComplexNDArray result (b.dims ()); | |
5275 | 1267 |
1268 for (octave_idx_type i = 0; i < b.length (); i++) | |
4543 | 1269 { |
1270 OCTAVE_QUIT; | |
5260 | 1271 result(i) = std::pow (a, b(i)); |
4543 | 1272 } |
1273 | |
1274 return result; | |
1275 } | |
1276 | |
1277 // -*- 9 -*- | |
1278 octave_value | |
1279 elem_xpow (const ComplexNDArray& a, double b) | |
1280 { | |
1281 ComplexNDArray result (a.dims ()); | |
1282 | |
1283 if (xisint (b)) | |
1284 { | |
5275 | 1285 for (octave_idx_type i = 0; i < a.length (); i++) |
4543 | 1286 { |
1287 OCTAVE_QUIT; | |
5260 | 1288 result(i) = std::pow (a(i), static_cast<int> (b)); |
4543 | 1289 } |
1290 } | |
1291 else | |
1292 { | |
5275 | 1293 for (octave_idx_type i = 0; i < a.length (); i++) |
4543 | 1294 { |
1295 OCTAVE_QUIT; | |
5260 | 1296 result(i) = std::pow (a(i), b); |
4543 | 1297 } |
1298 } | |
1299 | |
1300 return result; | |
1301 } | |
1302 | |
1303 // -*- 10 -*- | |
1304 octave_value | |
1305 elem_xpow (const ComplexNDArray& a, const NDArray& b) | |
1306 { | |
1307 dim_vector a_dims = a.dims (); | |
1308 dim_vector b_dims = b.dims (); | |
1309 | |
1310 if (a_dims != b_dims) | |
1311 { | |
1312 gripe_nonconformant ("operator .^", a_dims, b_dims); | |
1313 return octave_value (); | |
1314 } | |
1315 | |
1316 ComplexNDArray result (a_dims); | |
5275 | 1317 |
1318 for (octave_idx_type i = 0; i < a.length (); i++) | |
4543 | 1319 { |
1320 OCTAVE_QUIT; | |
1321 double btmp = b(i); | |
1322 if (xisint (btmp)) | |
5260 | 1323 result(i) = std::pow (a(i), static_cast<int> (btmp)); |
4543 | 1324 else |
5260 | 1325 result(i) = std::pow (a(i), btmp); |
4543 | 1326 } |
1327 | |
1328 return result; | |
1329 } | |
1330 | |
1331 // -*- 11 -*- | |
1332 octave_value | |
1333 elem_xpow (const ComplexNDArray& a, const Complex& b) | |
1334 { | |
1335 ComplexNDArray result (a.dims ()); | |
5275 | 1336 |
1337 for (octave_idx_type i = 0; i < a.length (); i++) | |
4543 | 1338 { |
1339 OCTAVE_QUIT; | |
5260 | 1340 result(i) = std::pow (a(i), b); |
4543 | 1341 } |
1342 | |
1343 return result; | |
1344 } | |
1345 | |
1346 // -*- 12 -*- | |
1347 octave_value | |
1348 elem_xpow (const ComplexNDArray& a, const ComplexNDArray& b) | |
1349 { | |
1350 dim_vector a_dims = a.dims (); | |
1351 dim_vector b_dims = b.dims (); | |
1352 | |
1353 if (a_dims != b_dims) | |
1354 { | |
1355 gripe_nonconformant ("operator .^", a_dims, b_dims); | |
1356 return octave_value (); | |
1357 } | |
1358 | |
1359 ComplexNDArray result (a_dims); | |
5275 | 1360 |
1361 for (octave_idx_type i = 0; i < a.length (); i++) | |
4543 | 1362 { |
1363 OCTAVE_QUIT; | |
5260 | 1364 result(i) = std::pow (a(i), b(i)); |
4543 | 1365 } |
1366 | |
1367 return result; | |
1368 } | |
1369 | |
7789
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1370 static inline int |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1371 xisint (float x) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1372 { |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1373 return (D_NINT (x) == x |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1374 && ((x >= 0 && x < INT_MAX) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1375 || (x <= 0 && x > INT_MIN))); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1376 } |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1377 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1378 // Safer pow functions. |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1379 // |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1380 // op2 \ op1: s m cs cm |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1381 // +-- +---+---+----+----+ |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1382 // scalar | | 1 | 5 | 7 | 11 | |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1383 // +---+---+----+----+ |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1384 // matrix | 2 | * | 8 | * | |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1385 // +---+---+----+----+ |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1386 // complex_scalar | 3 | 6 | 9 | 12 | |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1387 // +---+---+----+----+ |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1388 // complex_matrix | 4 | * | 10 | * | |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1389 // +---+---+----+----+ |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1390 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1391 // -*- 1 -*- |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1392 octave_value |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1393 xpow (float a, float b) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1394 { |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1395 float retval; |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1396 |
8978
0a58c4cd1405
optimize element-wise power
Jaroslav Hajek <highegg@gmail.com>
parents:
8960
diff
changeset
|
1397 if (a < 0.0 && ! xisint (b)) |
7789
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1398 { |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1399 FloatComplex atmp (a); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1400 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1401 return std::pow (atmp, b); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1402 } |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1403 else |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1404 retval = std::pow (a, b); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1405 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1406 return retval; |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1407 } |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1408 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1409 // -*- 2 -*- |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1410 octave_value |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1411 xpow (float a, const FloatMatrix& b) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1412 { |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1413 octave_value retval; |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1414 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1415 octave_idx_type nr = b.rows (); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1416 octave_idx_type nc = b.cols (); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1417 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1418 if (nr == 0 || nc == 0 || nr != nc) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1419 error ("for x^A, A must be square"); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1420 else |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1421 { |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1422 FloatEIG b_eig (b); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1423 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1424 if (! error_state) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1425 { |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1426 FloatComplexColumnVector lambda (b_eig.eigenvalues ()); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1427 FloatComplexMatrix Q (b_eig.eigenvectors ()); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1428 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1429 for (octave_idx_type i = 0; i < nr; i++) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1430 { |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1431 FloatComplex elt = lambda(i); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1432 if (std::imag (elt) == 0.0) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1433 lambda(i) = std::pow (a, std::real (elt)); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1434 else |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1435 lambda(i) = std::pow (a, elt); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1436 } |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1437 FloatComplexDiagMatrix D (lambda); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1438 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1439 retval = FloatComplexMatrix (Q * D * Q.inverse ()); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1440 } |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1441 else |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1442 error ("xpow: matrix diagonalization failed"); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1443 } |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1444 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1445 return retval; |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1446 } |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1447 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1448 // -*- 3 -*- |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1449 octave_value |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1450 xpow (float a, const FloatComplex& b) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1451 { |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1452 FloatComplex result; |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1453 FloatComplex atmp (a); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1454 result = std::pow (atmp, b); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1455 return result; |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1456 } |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1457 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1458 // -*- 4 -*- |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1459 octave_value |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1460 xpow (float a, const FloatComplexMatrix& b) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1461 { |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1462 octave_value retval; |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1463 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1464 octave_idx_type nr = b.rows (); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1465 octave_idx_type nc = b.cols (); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1466 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1467 if (nr == 0 || nc == 0 || nr != nc) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1468 error ("for x^A, A must be square"); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1469 else |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1470 { |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1471 FloatEIG b_eig (b); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1472 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1473 if (! error_state) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1474 { |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1475 FloatComplexColumnVector lambda (b_eig.eigenvalues ()); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1476 FloatComplexMatrix Q (b_eig.eigenvectors ()); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1477 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1478 for (octave_idx_type i = 0; i < nr; i++) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1479 { |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1480 FloatComplex elt = lambda(i); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1481 if (std::imag (elt) == 0.0) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1482 lambda(i) = std::pow (a, std::real (elt)); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1483 else |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1484 lambda(i) = std::pow (a, elt); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1485 } |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1486 FloatComplexDiagMatrix D (lambda); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1487 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1488 retval = FloatComplexMatrix (Q * D * Q.inverse ()); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1489 } |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1490 else |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1491 error ("xpow: matrix diagonalization failed"); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1492 } |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1493 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1494 return retval; |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1495 } |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1496 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1497 // -*- 5 -*- |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1498 octave_value |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1499 xpow (const FloatMatrix& a, float b) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1500 { |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1501 octave_value retval; |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1502 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1503 octave_idx_type nr = a.rows (); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1504 octave_idx_type nc = a.cols (); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1505 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1506 if (nr == 0 || nc == 0 || nr != nc) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1507 error ("for A^b, A must be square"); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1508 else |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1509 { |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1510 if (static_cast<int> (b) == b) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1511 { |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1512 int btmp = static_cast<int> (b); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1513 if (btmp == 0) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1514 { |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1515 retval = FloatDiagMatrix (nr, nr, 1.0); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1516 } |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1517 else |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1518 { |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1519 // Too much copying? |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1520 // FIXME -- we shouldn't do this if the exponent is |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1521 // large... |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1522 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1523 FloatMatrix atmp; |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1524 if (btmp < 0) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1525 { |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1526 btmp = -btmp; |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1527 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1528 octave_idx_type info; |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1529 float rcond = 0.0; |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1530 MatrixType mattype (a); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1531 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1532 atmp = a.inverse (mattype, info, rcond, 1); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1533 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1534 if (info == -1) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1535 warning ("inverse: matrix singular to machine\ |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1536 precision, rcond = %g", rcond); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1537 } |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1538 else |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1539 atmp = a; |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1540 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1541 FloatMatrix result (atmp); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1542 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1543 btmp--; |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1544 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1545 while (btmp > 0) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1546 { |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1547 if (btmp & 1) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1548 result = result * atmp; |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1549 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1550 btmp >>= 1; |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1551 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1552 if (btmp > 0) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1553 atmp = atmp * atmp; |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1554 } |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1555 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1556 retval = result; |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1557 } |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1558 } |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1559 else |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1560 { |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1561 FloatEIG a_eig (a); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1562 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1563 if (! error_state) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1564 { |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1565 FloatComplexColumnVector lambda (a_eig.eigenvalues ()); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1566 FloatComplexMatrix Q (a_eig.eigenvectors ()); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1567 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1568 for (octave_idx_type i = 0; i < nr; i++) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1569 lambda(i) = std::pow (lambda(i), b); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1570 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1571 FloatComplexDiagMatrix D (lambda); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1572 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1573 retval = FloatComplexMatrix (Q * D * Q.inverse ()); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1574 } |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1575 else |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1576 error ("xpow: matrix diagonalization failed"); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1577 } |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1578 } |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1579 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1580 return retval; |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1581 } |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1582 |
8382
9b20a4847056
implement scalar powers of diag matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
1583 // -*- 5d -*- |
9b20a4847056
implement scalar powers of diag matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
1584 octave_value |
9b20a4847056
implement scalar powers of diag matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
1585 xpow (const FloatDiagMatrix& a, float b) |
9b20a4847056
implement scalar powers of diag matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
1586 { |
9b20a4847056
implement scalar powers of diag matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
1587 octave_value retval; |
9b20a4847056
implement scalar powers of diag matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
1588 |
9b20a4847056
implement scalar powers of diag matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
1589 octave_idx_type nr = a.rows (); |
9b20a4847056
implement scalar powers of diag matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
1590 octave_idx_type nc = a.cols (); |
9b20a4847056
implement scalar powers of diag matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
1591 |
9b20a4847056
implement scalar powers of diag matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
1592 if (nr == 0 || nc == 0 || nr != nc) |
9b20a4847056
implement scalar powers of diag matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
1593 error ("for A^b, A must be square"); |
9b20a4847056
implement scalar powers of diag matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
1594 else |
9b20a4847056
implement scalar powers of diag matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
1595 { |
9b20a4847056
implement scalar powers of diag matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
1596 if (static_cast<int> (b) == b) |
9b20a4847056
implement scalar powers of diag matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
1597 { |
9b20a4847056
implement scalar powers of diag matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
1598 FloatDiagMatrix r (nr, nc); |
9b20a4847056
implement scalar powers of diag matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
1599 for (octave_idx_type i = 0; i < nc; i++) |
9b20a4847056
implement scalar powers of diag matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
1600 r(i, i) = std::pow (a(i, i), b); |
9b20a4847056
implement scalar powers of diag matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
1601 retval = r; |
9b20a4847056
implement scalar powers of diag matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
1602 } |
9b20a4847056
implement scalar powers of diag matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
1603 else |
9b20a4847056
implement scalar powers of diag matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
1604 { |
9b20a4847056
implement scalar powers of diag matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
1605 FloatComplexDiagMatrix r (nr, nc); |
9b20a4847056
implement scalar powers of diag matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
1606 for (octave_idx_type i = 0; i < nc; i++) |
9b20a4847056
implement scalar powers of diag matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
1607 r(i, i) = std::pow (static_cast<FloatComplex> (a(i, i)), b); |
9b20a4847056
implement scalar powers of diag matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
1608 retval = r; |
9b20a4847056
implement scalar powers of diag matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
1609 } |
9b20a4847056
implement scalar powers of diag matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
1610 } |
9b20a4847056
implement scalar powers of diag matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
1611 |
9b20a4847056
implement scalar powers of diag matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
1612 return retval; |
9b20a4847056
implement scalar powers of diag matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
1613 } |
9b20a4847056
implement scalar powers of diag matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
1614 |
7789
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1615 // -*- 6 -*- |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1616 octave_value |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1617 xpow (const FloatMatrix& a, const FloatComplex& b) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1618 { |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1619 octave_value retval; |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1620 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1621 octave_idx_type nr = a.rows (); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1622 octave_idx_type nc = a.cols (); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1623 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1624 if (nr == 0 || nc == 0 || nr != nc) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1625 error ("for A^b, A must be square"); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1626 else |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1627 { |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1628 FloatEIG a_eig (a); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1629 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1630 if (! error_state) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1631 { |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1632 FloatComplexColumnVector lambda (a_eig.eigenvalues ()); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1633 FloatComplexMatrix Q (a_eig.eigenvectors ()); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1634 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1635 for (octave_idx_type i = 0; i < nr; i++) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1636 lambda(i) = std::pow (lambda(i), b); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1637 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1638 FloatComplexDiagMatrix D (lambda); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1639 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1640 retval = FloatComplexMatrix (Q * D * Q.inverse ()); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1641 } |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1642 else |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1643 error ("xpow: matrix diagonalization failed"); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1644 } |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1645 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1646 return retval; |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1647 } |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1648 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1649 // -*- 7 -*- |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1650 octave_value |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1651 xpow (const FloatComplex& a, float b) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1652 { |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1653 FloatComplex result; |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1654 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1655 if (xisint (b)) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1656 result = std::pow (a, static_cast<int> (b)); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1657 else |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1658 result = std::pow (a, b); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1659 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1660 return result; |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1661 } |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1662 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1663 // -*- 8 -*- |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1664 octave_value |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1665 xpow (const FloatComplex& a, const FloatMatrix& b) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1666 { |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1667 octave_value retval; |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1668 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1669 octave_idx_type nr = b.rows (); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1670 octave_idx_type nc = b.cols (); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1671 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1672 if (nr == 0 || nc == 0 || nr != nc) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1673 error ("for x^A, A must be square"); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1674 else |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1675 { |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1676 FloatEIG b_eig (b); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1677 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1678 if (! error_state) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1679 { |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1680 FloatComplexColumnVector lambda (b_eig.eigenvalues ()); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1681 FloatComplexMatrix Q (b_eig.eigenvectors ()); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1682 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1683 for (octave_idx_type i = 0; i < nr; i++) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1684 { |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1685 FloatComplex elt = lambda(i); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1686 if (std::imag (elt) == 0.0) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1687 lambda(i) = std::pow (a, std::real (elt)); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1688 else |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1689 lambda(i) = std::pow (a, elt); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1690 } |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1691 FloatComplexDiagMatrix D (lambda); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1692 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1693 retval = FloatComplexMatrix (Q * D * Q.inverse ()); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1694 } |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1695 else |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1696 error ("xpow: matrix diagonalization failed"); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1697 } |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1698 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1699 return retval; |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1700 } |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1701 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1702 // -*- 9 -*- |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1703 octave_value |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1704 xpow (const FloatComplex& a, const FloatComplex& b) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1705 { |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1706 FloatComplex result; |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1707 result = std::pow (a, b); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1708 return result; |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1709 } |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1710 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1711 // -*- 10 -*- |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1712 octave_value |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1713 xpow (const FloatComplex& a, const FloatComplexMatrix& b) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1714 { |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1715 octave_value retval; |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1716 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1717 octave_idx_type nr = b.rows (); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1718 octave_idx_type nc = b.cols (); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1719 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1720 if (nr == 0 || nc == 0 || nr != nc) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1721 error ("for x^A, A must be square"); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1722 else |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1723 { |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1724 FloatEIG b_eig (b); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1725 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1726 if (! error_state) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1727 { |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1728 FloatComplexColumnVector lambda (b_eig.eigenvalues ()); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1729 FloatComplexMatrix Q (b_eig.eigenvectors ()); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1730 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1731 for (octave_idx_type i = 0; i < nr; i++) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1732 { |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1733 FloatComplex elt = lambda(i); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1734 if (std::imag (elt) == 0.0) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1735 lambda(i) = std::pow (a, std::real (elt)); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1736 else |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1737 lambda(i) = std::pow (a, elt); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1738 } |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1739 FloatComplexDiagMatrix D (lambda); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1740 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1741 retval = FloatComplexMatrix (Q * D * Q.inverse ()); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1742 } |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1743 else |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1744 error ("xpow: matrix diagonalization failed"); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1745 } |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1746 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1747 return retval; |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1748 } |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1749 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1750 // -*- 11 -*- |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1751 octave_value |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1752 xpow (const FloatComplexMatrix& a, float b) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1753 { |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1754 octave_value retval; |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1755 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1756 octave_idx_type nr = a.rows (); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1757 octave_idx_type nc = a.cols (); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1758 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1759 if (nr == 0 || nc == 0 || nr != nc) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1760 error ("for A^b, A must be square"); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1761 else |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1762 { |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1763 if (static_cast<int> (b) == b) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1764 { |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1765 int btmp = static_cast<int> (b); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1766 if (btmp == 0) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1767 { |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1768 retval = FloatDiagMatrix (nr, nr, 1.0); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1769 } |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1770 else |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1771 { |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1772 // Too much copying? |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1773 // FIXME -- we shouldn't do this if the exponent is |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1774 // large... |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1775 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1776 FloatComplexMatrix atmp; |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1777 if (btmp < 0) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1778 { |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1779 btmp = -btmp; |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1780 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1781 octave_idx_type info; |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1782 float rcond = 0.0; |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1783 MatrixType mattype (a); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1784 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1785 atmp = a.inverse (mattype, info, rcond, 1); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1786 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1787 if (info == -1) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1788 warning ("inverse: matrix singular to machine\ |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1789 precision, rcond = %g", rcond); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1790 } |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1791 else |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1792 atmp = a; |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1793 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1794 FloatComplexMatrix result (atmp); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1795 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1796 btmp--; |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1797 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1798 while (btmp > 0) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1799 { |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1800 if (btmp & 1) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1801 result = result * atmp; |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1802 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1803 btmp >>= 1; |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1804 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1805 if (btmp > 0) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1806 atmp = atmp * atmp; |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1807 } |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1808 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1809 retval = result; |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1810 } |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1811 } |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1812 else |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1813 { |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1814 FloatEIG a_eig (a); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1815 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1816 if (! error_state) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1817 { |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1818 FloatComplexColumnVector lambda (a_eig.eigenvalues ()); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1819 FloatComplexMatrix Q (a_eig.eigenvectors ()); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1820 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1821 for (octave_idx_type i = 0; i < nr; i++) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1822 lambda(i) = std::pow (lambda(i), b); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1823 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1824 FloatComplexDiagMatrix D (lambda); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1825 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1826 retval = FloatComplexMatrix (Q * D * Q.inverse ()); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1827 } |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1828 else |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1829 error ("xpow: matrix diagonalization failed"); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1830 } |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1831 } |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1832 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1833 return retval; |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1834 } |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1835 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1836 // -*- 12 -*- |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1837 octave_value |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1838 xpow (const FloatComplexMatrix& a, const FloatComplex& b) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1839 { |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1840 octave_value retval; |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1841 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1842 octave_idx_type nr = a.rows (); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1843 octave_idx_type nc = a.cols (); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1844 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1845 if (nr == 0 || nc == 0 || nr != nc) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1846 error ("for A^b, A must be square"); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1847 else |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1848 { |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1849 FloatEIG a_eig (a); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1850 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1851 if (! error_state) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1852 { |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1853 FloatComplexColumnVector lambda (a_eig.eigenvalues ()); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1854 FloatComplexMatrix Q (a_eig.eigenvectors ()); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1855 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1856 for (octave_idx_type i = 0; i < nr; i++) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1857 lambda(i) = std::pow (lambda(i), b); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1858 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1859 FloatComplexDiagMatrix D (lambda); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1860 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1861 retval = FloatComplexMatrix (Q * D * Q.inverse ()); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1862 } |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1863 else |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1864 error ("xpow: matrix diagonalization failed"); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1865 } |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1866 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1867 return retval; |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1868 } |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1869 |
8382
9b20a4847056
implement scalar powers of diag matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
1870 // -*- 12d -*- |
9b20a4847056
implement scalar powers of diag matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
1871 octave_value |
9b20a4847056
implement scalar powers of diag matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
1872 xpow (const FloatComplexDiagMatrix& a, const FloatComplex& b) |
9b20a4847056
implement scalar powers of diag matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
1873 { |
9b20a4847056
implement scalar powers of diag matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
1874 octave_value retval; |
9b20a4847056
implement scalar powers of diag matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
1875 |
9b20a4847056
implement scalar powers of diag matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
1876 octave_idx_type nr = a.rows (); |
9b20a4847056
implement scalar powers of diag matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
1877 octave_idx_type nc = a.cols (); |
9b20a4847056
implement scalar powers of diag matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
1878 |
9b20a4847056
implement scalar powers of diag matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
1879 if (nr == 0 || nc == 0 || nr != nc) |
9b20a4847056
implement scalar powers of diag matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
1880 error ("for A^b, A must be square"); |
9b20a4847056
implement scalar powers of diag matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
1881 else |
9b20a4847056
implement scalar powers of diag matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
1882 { |
9b20a4847056
implement scalar powers of diag matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
1883 FloatComplexDiagMatrix r (nr, nc); |
9b20a4847056
implement scalar powers of diag matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
1884 for (octave_idx_type i = 0; i < nc; i++) |
9b20a4847056
implement scalar powers of diag matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
1885 r(i, i) = std::pow (a(i, i), b); |
9b20a4847056
implement scalar powers of diag matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
1886 retval = r; |
9b20a4847056
implement scalar powers of diag matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
1887 } |
9b20a4847056
implement scalar powers of diag matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
1888 |
9b20a4847056
implement scalar powers of diag matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
1889 return retval; |
9b20a4847056
implement scalar powers of diag matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
1890 } |
9b20a4847056
implement scalar powers of diag matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
1891 |
9b20a4847056
implement scalar powers of diag matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
1892 // mixed |
9b20a4847056
implement scalar powers of diag matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
1893 octave_value |
9b20a4847056
implement scalar powers of diag matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
1894 xpow (const FloatComplexDiagMatrix& a, float b) |
9b20a4847056
implement scalar powers of diag matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
1895 { |
9b20a4847056
implement scalar powers of diag matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
1896 return xpow (a, static_cast<FloatComplex> (b)); |
9b20a4847056
implement scalar powers of diag matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
1897 } |
9b20a4847056
implement scalar powers of diag matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
1898 |
9b20a4847056
implement scalar powers of diag matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
1899 octave_value |
9b20a4847056
implement scalar powers of diag matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
1900 xpow (const FloatDiagMatrix& a, const FloatComplex& b) |
9b20a4847056
implement scalar powers of diag matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
1901 { |
9b20a4847056
implement scalar powers of diag matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
1902 return xpow (FloatComplexDiagMatrix (a), b); |
9b20a4847056
implement scalar powers of diag matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
1903 } |
9b20a4847056
implement scalar powers of diag matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
1904 |
7789
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1905 // Safer pow functions that work elementwise for matrices. |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1906 // |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1907 // op2 \ op1: s m cs cm |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1908 // +-- +---+---+----+----+ |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1909 // scalar | | * | 3 | * | 9 | |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1910 // +---+---+----+----+ |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1911 // matrix | 1 | 4 | 7 | 10 | |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1912 // +---+---+----+----+ |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1913 // complex_scalar | * | 5 | * | 11 | |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1914 // +---+---+----+----+ |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1915 // complex_matrix | 2 | 6 | 8 | 12 | |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1916 // +---+---+----+----+ |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1917 // |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1918 // * -> not needed. |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1919 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1920 // FIXME -- these functions need to be fixed so that things |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1921 // like |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1922 // |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1923 // a = -1; b = [ 0, 0.5, 1 ]; r = a .^ b |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1924 // |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1925 // and |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1926 // |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1927 // a = -1; b = [ 0, 0.5, 1 ]; for i = 1:3, r(i) = a .^ b(i), end |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1928 // |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1929 // produce identical results. Also, it would be nice if -1^0.5 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1930 // produced a pure imaginary result instead of a complex number with a |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1931 // small real part. But perhaps that's really a problem with the math |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1932 // library... |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1933 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1934 // -*- 1 -*- |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1935 octave_value |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1936 elem_xpow (float a, const FloatMatrix& b) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1937 { |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1938 octave_value retval; |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1939 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1940 octave_idx_type nr = b.rows (); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1941 octave_idx_type nc = b.cols (); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1942 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1943 float d1, d2; |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1944 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1945 if (a < 0.0 && ! b.all_integers (d1, d2)) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1946 { |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1947 FloatComplex atmp (a); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1948 FloatComplexMatrix result (nr, nc); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1949 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1950 for (octave_idx_type j = 0; j < nc; j++) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1951 for (octave_idx_type i = 0; i < nr; i++) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1952 { |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1953 OCTAVE_QUIT; |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1954 result (i, j) = std::pow (atmp, b (i, j)); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1955 } |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1956 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1957 retval = result; |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1958 } |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1959 else |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1960 { |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1961 FloatMatrix result (nr, nc); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1962 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1963 for (octave_idx_type j = 0; j < nc; j++) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1964 for (octave_idx_type i = 0; i < nr; i++) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1965 { |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1966 OCTAVE_QUIT; |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1967 result (i, j) = std::pow (a, b (i, j)); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1968 } |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1969 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1970 retval = result; |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1971 } |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1972 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1973 return retval; |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1974 } |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1975 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1976 // -*- 2 -*- |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1977 octave_value |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1978 elem_xpow (float a, const FloatComplexMatrix& b) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1979 { |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1980 octave_idx_type nr = b.rows (); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1981 octave_idx_type nc = b.cols (); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1982 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1983 FloatComplexMatrix result (nr, nc); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1984 FloatComplex atmp (a); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1985 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1986 for (octave_idx_type j = 0; j < nc; j++) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1987 for (octave_idx_type i = 0; i < nr; i++) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1988 { |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1989 OCTAVE_QUIT; |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1990 result (i, j) = std::pow (atmp, b (i, j)); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1991 } |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1992 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1993 return result; |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1994 } |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1995 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1996 // -*- 3 -*- |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1997 octave_value |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1998 elem_xpow (const FloatMatrix& a, float b) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
1999 { |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2000 octave_value retval; |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2001 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2002 octave_idx_type nr = a.rows (); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2003 octave_idx_type nc = a.cols (); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2004 |
8978
0a58c4cd1405
optimize element-wise power
Jaroslav Hajek <highegg@gmail.com>
parents:
8960
diff
changeset
|
2005 if (! xisint (b) && a.any_element_is_negative ()) |
7789
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2006 { |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2007 FloatComplexMatrix result (nr, nc); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2008 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2009 for (octave_idx_type j = 0; j < nc; j++) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2010 for (octave_idx_type i = 0; i < nr; i++) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2011 { |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2012 OCTAVE_QUIT; |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2013 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2014 FloatComplex atmp (a (i, j)); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2015 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2016 result (i, j) = std::pow (atmp, b); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2017 } |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2018 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2019 retval = result; |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2020 } |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2021 else |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2022 { |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2023 FloatMatrix result (nr, nc); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2024 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2025 for (octave_idx_type j = 0; j < nc; j++) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2026 for (octave_idx_type i = 0; i < nr; i++) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2027 { |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2028 OCTAVE_QUIT; |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2029 result (i, j) = std::pow (a (i, j), b); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2030 } |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2031 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2032 retval = result; |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2033 } |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2034 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2035 return retval; |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2036 } |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2037 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2038 // -*- 4 -*- |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2039 octave_value |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2040 elem_xpow (const FloatMatrix& a, const FloatMatrix& b) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2041 { |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2042 octave_value retval; |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2043 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2044 octave_idx_type nr = a.rows (); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2045 octave_idx_type nc = a.cols (); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2046 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2047 octave_idx_type b_nr = b.rows (); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2048 octave_idx_type b_nc = b.cols (); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2049 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2050 if (nr != b_nr || nc != b_nc) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2051 { |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2052 gripe_nonconformant ("operator .^", nr, nc, b_nr, b_nc); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2053 return octave_value (); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2054 } |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2055 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2056 int convert_to_complex = 0; |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2057 for (octave_idx_type j = 0; j < nc; j++) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2058 for (octave_idx_type i = 0; i < nr; i++) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2059 { |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2060 OCTAVE_QUIT; |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2061 float atmp = a (i, j); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2062 float btmp = b (i, j); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2063 if (atmp < 0.0 && static_cast<int> (btmp) != btmp) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2064 { |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2065 convert_to_complex = 1; |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2066 goto done; |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2067 } |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2068 } |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2069 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2070 done: |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2071 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2072 if (convert_to_complex) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2073 { |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2074 FloatComplexMatrix complex_result (nr, nc); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2075 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2076 for (octave_idx_type j = 0; j < nc; j++) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2077 for (octave_idx_type i = 0; i < nr; i++) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2078 { |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2079 OCTAVE_QUIT; |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2080 FloatComplex atmp (a (i, j)); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2081 FloatComplex btmp (b (i, j)); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2082 complex_result (i, j) = std::pow (atmp, btmp); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2083 } |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2084 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2085 retval = complex_result; |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2086 } |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2087 else |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2088 { |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2089 FloatMatrix result (nr, nc); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2090 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2091 for (octave_idx_type j = 0; j < nc; j++) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2092 for (octave_idx_type i = 0; i < nr; i++) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2093 { |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2094 OCTAVE_QUIT; |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2095 result (i, j) = std::pow (a (i, j), b (i, j)); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2096 } |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2097 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2098 retval = result; |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2099 } |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2100 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2101 return retval; |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2102 } |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2103 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2104 // -*- 5 -*- |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2105 octave_value |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2106 elem_xpow (const FloatMatrix& a, const FloatComplex& b) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2107 { |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2108 octave_idx_type nr = a.rows (); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2109 octave_idx_type nc = a.cols (); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2110 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2111 FloatComplexMatrix result (nr, nc); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2112 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2113 for (octave_idx_type j = 0; j < nc; j++) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2114 for (octave_idx_type i = 0; i < nr; i++) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2115 { |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2116 OCTAVE_QUIT; |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2117 result (i, j) = std::pow (FloatComplex (a (i, j)), b); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2118 } |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2119 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2120 return result; |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2121 } |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2122 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2123 // -*- 6 -*- |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2124 octave_value |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2125 elem_xpow (const FloatMatrix& a, const FloatComplexMatrix& b) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2126 { |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2127 octave_idx_type nr = a.rows (); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2128 octave_idx_type nc = a.cols (); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2129 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2130 octave_idx_type b_nr = b.rows (); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2131 octave_idx_type b_nc = b.cols (); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2132 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2133 if (nr != b_nr || nc != b_nc) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2134 { |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2135 gripe_nonconformant ("operator .^", nr, nc, b_nr, b_nc); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2136 return octave_value (); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2137 } |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2138 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2139 FloatComplexMatrix result (nr, nc); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2140 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2141 for (octave_idx_type j = 0; j < nc; j++) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2142 for (octave_idx_type i = 0; i < nr; i++) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2143 { |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2144 OCTAVE_QUIT; |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2145 result (i, j) = std::pow (FloatComplex (a (i, j)), b (i, j)); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2146 } |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2147 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2148 return result; |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2149 } |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2150 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2151 // -*- 7 -*- |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2152 octave_value |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2153 elem_xpow (const FloatComplex& a, const FloatMatrix& b) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2154 { |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2155 octave_idx_type nr = b.rows (); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2156 octave_idx_type nc = b.cols (); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2157 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2158 FloatComplexMatrix result (nr, nc); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2159 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2160 for (octave_idx_type j = 0; j < nc; j++) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2161 for (octave_idx_type i = 0; i < nr; i++) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2162 { |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2163 OCTAVE_QUIT; |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2164 float btmp = b (i, j); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2165 if (xisint (btmp)) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2166 result (i, j) = std::pow (a, static_cast<int> (btmp)); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2167 else |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2168 result (i, j) = std::pow (a, btmp); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2169 } |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2170 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2171 return result; |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2172 } |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2173 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2174 // -*- 8 -*- |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2175 octave_value |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2176 elem_xpow (const FloatComplex& a, const FloatComplexMatrix& b) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2177 { |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2178 octave_idx_type nr = b.rows (); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2179 octave_idx_type nc = b.cols (); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2180 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2181 FloatComplexMatrix result (nr, nc); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2182 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2183 for (octave_idx_type j = 0; j < nc; j++) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2184 for (octave_idx_type i = 0; i < nr; i++) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2185 { |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2186 OCTAVE_QUIT; |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2187 result (i, j) = std::pow (a, b (i, j)); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2188 } |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2189 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2190 return result; |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2191 } |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2192 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2193 // -*- 9 -*- |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2194 octave_value |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2195 elem_xpow (const FloatComplexMatrix& a, float b) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2196 { |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2197 octave_idx_type nr = a.rows (); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2198 octave_idx_type nc = a.cols (); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2199 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2200 FloatComplexMatrix result (nr, nc); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2201 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2202 if (xisint (b)) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2203 { |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2204 for (octave_idx_type j = 0; j < nc; j++) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2205 for (octave_idx_type i = 0; i < nr; i++) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2206 { |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2207 OCTAVE_QUIT; |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2208 result (i, j) = std::pow (a (i, j), static_cast<int> (b)); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2209 } |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2210 } |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2211 else |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2212 { |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2213 for (octave_idx_type j = 0; j < nc; j++) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2214 for (octave_idx_type i = 0; i < nr; i++) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2215 { |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2216 OCTAVE_QUIT; |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2217 result (i, j) = std::pow (a (i, j), b); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2218 } |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2219 } |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2220 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2221 return result; |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2222 } |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2223 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2224 // -*- 10 -*- |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2225 octave_value |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2226 elem_xpow (const FloatComplexMatrix& a, const FloatMatrix& b) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2227 { |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2228 octave_idx_type nr = a.rows (); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2229 octave_idx_type nc = a.cols (); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2230 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2231 octave_idx_type b_nr = b.rows (); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2232 octave_idx_type b_nc = b.cols (); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2233 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2234 if (nr != b_nr || nc != b_nc) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2235 { |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2236 gripe_nonconformant ("operator .^", nr, nc, b_nr, b_nc); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2237 return octave_value (); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2238 } |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2239 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2240 FloatComplexMatrix result (nr, nc); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2241 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2242 for (octave_idx_type j = 0; j < nc; j++) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2243 for (octave_idx_type i = 0; i < nr; i++) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2244 { |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2245 OCTAVE_QUIT; |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2246 float btmp = b (i, j); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2247 if (xisint (btmp)) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2248 result (i, j) = std::pow (a (i, j), static_cast<int> (btmp)); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2249 else |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2250 result (i, j) = std::pow (a (i, j), btmp); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2251 } |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2252 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2253 return result; |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2254 } |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2255 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2256 // -*- 11 -*- |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2257 octave_value |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2258 elem_xpow (const FloatComplexMatrix& a, const FloatComplex& b) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2259 { |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2260 octave_idx_type nr = a.rows (); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2261 octave_idx_type nc = a.cols (); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2262 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2263 FloatComplexMatrix result (nr, nc); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2264 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2265 for (octave_idx_type j = 0; j < nc; j++) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2266 for (octave_idx_type i = 0; i < nr; i++) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2267 { |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2268 OCTAVE_QUIT; |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2269 result (i, j) = std::pow (a (i, j), b); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2270 } |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2271 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2272 return result; |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2273 } |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2274 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2275 // -*- 12 -*- |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2276 octave_value |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2277 elem_xpow (const FloatComplexMatrix& a, const FloatComplexMatrix& b) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2278 { |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2279 octave_idx_type nr = a.rows (); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2280 octave_idx_type nc = a.cols (); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2281 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2282 octave_idx_type b_nr = b.rows (); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2283 octave_idx_type b_nc = b.cols (); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2284 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2285 if (nr != b_nr || nc != b_nc) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2286 { |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2287 gripe_nonconformant ("operator .^", nr, nc, b_nr, b_nc); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2288 return octave_value (); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2289 } |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2290 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2291 FloatComplexMatrix result (nr, nc); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2292 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2293 for (octave_idx_type j = 0; j < nc; j++) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2294 for (octave_idx_type i = 0; i < nr; i++) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2295 { |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2296 OCTAVE_QUIT; |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2297 result (i, j) = std::pow (a (i, j), b (i, j)); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2298 } |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2299 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2300 return result; |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2301 } |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2302 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2303 // Safer pow functions that work elementwise for N-d arrays. |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2304 // |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2305 // op2 \ op1: s nd cs cnd |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2306 // +-- +---+---+----+----+ |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2307 // scalar | | * | 3 | * | 9 | |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2308 // +---+---+----+----+ |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2309 // N_d | 1 | 4 | 7 | 10 | |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2310 // +---+---+----+----+ |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2311 // complex_scalar | * | 5 | * | 11 | |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2312 // +---+---+----+----+ |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2313 // complex_N_d | 2 | 6 | 8 | 12 | |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2314 // +---+---+----+----+ |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2315 // |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2316 // * -> not needed. |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2317 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2318 // FIXME -- these functions need to be fixed so that things |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2319 // like |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2320 // |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2321 // a = -1; b = [ 0, 0.5, 1 ]; r = a .^ b |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2322 // |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2323 // and |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2324 // |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2325 // a = -1; b = [ 0, 0.5, 1 ]; for i = 1:3, r(i) = a .^ b(i), end |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2326 // |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2327 // produce identical results. Also, it would be nice if -1^0.5 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2328 // produced a pure imaginary result instead of a complex number with a |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2329 // small real part. But perhaps that's really a problem with the math |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2330 // library... |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2331 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2332 // -*- 1 -*- |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2333 octave_value |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2334 elem_xpow (float a, const FloatNDArray& b) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2335 { |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2336 octave_value retval; |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2337 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2338 float d1, d2; |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2339 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2340 if (a < 0.0 && ! b.all_integers (d1, d2)) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2341 { |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2342 FloatComplex atmp (a); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2343 FloatComplexNDArray result (b.dims ()); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2344 for (octave_idx_type i = 0; i < b.length (); i++) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2345 { |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2346 OCTAVE_QUIT; |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2347 result(i) = std::pow (atmp, b(i)); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2348 } |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2349 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2350 retval = result; |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2351 } |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2352 else |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2353 { |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2354 FloatNDArray result (b.dims ()); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2355 for (octave_idx_type i = 0; i < b.length (); i++) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2356 { |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2357 OCTAVE_QUIT; |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2358 result (i) = std::pow (a, b(i)); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2359 } |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2360 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2361 retval = result; |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2362 } |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2363 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2364 return retval; |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2365 } |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2366 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2367 // -*- 2 -*- |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2368 octave_value |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2369 elem_xpow (float a, const FloatComplexNDArray& b) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2370 { |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2371 FloatComplexNDArray result (b.dims ()); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2372 FloatComplex atmp (a); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2373 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2374 for (octave_idx_type i = 0; i < b.length (); i++) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2375 { |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2376 OCTAVE_QUIT; |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2377 result(i) = std::pow (atmp, b(i)); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2378 } |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2379 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2380 return result; |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2381 } |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2382 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2383 // -*- 3 -*- |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2384 octave_value |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2385 elem_xpow (const FloatNDArray& a, float b) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2386 { |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2387 octave_value retval; |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2388 |
8978
0a58c4cd1405
optimize element-wise power
Jaroslav Hajek <highegg@gmail.com>
parents:
8960
diff
changeset
|
2389 if (! xisint (b)) |
7789
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2390 { |
8978
0a58c4cd1405
optimize element-wise power
Jaroslav Hajek <highegg@gmail.com>
parents:
8960
diff
changeset
|
2391 if (a.any_element_is_negative ()) |
0a58c4cd1405
optimize element-wise power
Jaroslav Hajek <highegg@gmail.com>
parents:
8960
diff
changeset
|
2392 { |
0a58c4cd1405
optimize element-wise power
Jaroslav Hajek <highegg@gmail.com>
parents:
8960
diff
changeset
|
2393 FloatComplexNDArray result (a.dims ()); |
0a58c4cd1405
optimize element-wise power
Jaroslav Hajek <highegg@gmail.com>
parents:
8960
diff
changeset
|
2394 |
0a58c4cd1405
optimize element-wise power
Jaroslav Hajek <highegg@gmail.com>
parents:
8960
diff
changeset
|
2395 for (octave_idx_type i = 0; i < a.length (); i++) |
0a58c4cd1405
optimize element-wise power
Jaroslav Hajek <highegg@gmail.com>
parents:
8960
diff
changeset
|
2396 { |
0a58c4cd1405
optimize element-wise power
Jaroslav Hajek <highegg@gmail.com>
parents:
8960
diff
changeset
|
2397 OCTAVE_QUIT; |
0a58c4cd1405
optimize element-wise power
Jaroslav Hajek <highegg@gmail.com>
parents:
8960
diff
changeset
|
2398 |
0a58c4cd1405
optimize element-wise power
Jaroslav Hajek <highegg@gmail.com>
parents:
8960
diff
changeset
|
2399 FloatComplex atmp (a (i)); |
0a58c4cd1405
optimize element-wise power
Jaroslav Hajek <highegg@gmail.com>
parents:
8960
diff
changeset
|
2400 |
0a58c4cd1405
optimize element-wise power
Jaroslav Hajek <highegg@gmail.com>
parents:
8960
diff
changeset
|
2401 result(i) = std::pow (atmp, b); |
0a58c4cd1405
optimize element-wise power
Jaroslav Hajek <highegg@gmail.com>
parents:
8960
diff
changeset
|
2402 } |
0a58c4cd1405
optimize element-wise power
Jaroslav Hajek <highegg@gmail.com>
parents:
8960
diff
changeset
|
2403 |
0a58c4cd1405
optimize element-wise power
Jaroslav Hajek <highegg@gmail.com>
parents:
8960
diff
changeset
|
2404 retval = result; |
0a58c4cd1405
optimize element-wise power
Jaroslav Hajek <highegg@gmail.com>
parents:
8960
diff
changeset
|
2405 } |
0a58c4cd1405
optimize element-wise power
Jaroslav Hajek <highegg@gmail.com>
parents:
8960
diff
changeset
|
2406 else |
0a58c4cd1405
optimize element-wise power
Jaroslav Hajek <highegg@gmail.com>
parents:
8960
diff
changeset
|
2407 { |
0a58c4cd1405
optimize element-wise power
Jaroslav Hajek <highegg@gmail.com>
parents:
8960
diff
changeset
|
2408 FloatNDArray result (a.dims ()); |
0a58c4cd1405
optimize element-wise power
Jaroslav Hajek <highegg@gmail.com>
parents:
8960
diff
changeset
|
2409 for (octave_idx_type i = 0; i < a.length (); i++) |
0a58c4cd1405
optimize element-wise power
Jaroslav Hajek <highegg@gmail.com>
parents:
8960
diff
changeset
|
2410 { |
0a58c4cd1405
optimize element-wise power
Jaroslav Hajek <highegg@gmail.com>
parents:
8960
diff
changeset
|
2411 OCTAVE_QUIT; |
0a58c4cd1405
optimize element-wise power
Jaroslav Hajek <highegg@gmail.com>
parents:
8960
diff
changeset
|
2412 result(i) = std::pow (a(i), b); |
0a58c4cd1405
optimize element-wise power
Jaroslav Hajek <highegg@gmail.com>
parents:
8960
diff
changeset
|
2413 } |
0a58c4cd1405
optimize element-wise power
Jaroslav Hajek <highegg@gmail.com>
parents:
8960
diff
changeset
|
2414 |
0a58c4cd1405
optimize element-wise power
Jaroslav Hajek <highegg@gmail.com>
parents:
8960
diff
changeset
|
2415 retval = result; |
0a58c4cd1405
optimize element-wise power
Jaroslav Hajek <highegg@gmail.com>
parents:
8960
diff
changeset
|
2416 } |
7789
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2417 } |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2418 else |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2419 { |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2420 FloatNDArray result (a.dims ()); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2421 |
8978
0a58c4cd1405
optimize element-wise power
Jaroslav Hajek <highegg@gmail.com>
parents:
8960
diff
changeset
|
2422 int ib = static_cast<int> (b); |
0a58c4cd1405
optimize element-wise power
Jaroslav Hajek <highegg@gmail.com>
parents:
8960
diff
changeset
|
2423 if (ib == 2) |
0a58c4cd1405
optimize element-wise power
Jaroslav Hajek <highegg@gmail.com>
parents:
8960
diff
changeset
|
2424 { |
0a58c4cd1405
optimize element-wise power
Jaroslav Hajek <highegg@gmail.com>
parents:
8960
diff
changeset
|
2425 for (octave_idx_type i = 0; i < a.length (); i++) |
0a58c4cd1405
optimize element-wise power
Jaroslav Hajek <highegg@gmail.com>
parents:
8960
diff
changeset
|
2426 result.xelem (i) = a(i) * a(i); |
0a58c4cd1405
optimize element-wise power
Jaroslav Hajek <highegg@gmail.com>
parents:
8960
diff
changeset
|
2427 } |
0a58c4cd1405
optimize element-wise power
Jaroslav Hajek <highegg@gmail.com>
parents:
8960
diff
changeset
|
2428 else |
0a58c4cd1405
optimize element-wise power
Jaroslav Hajek <highegg@gmail.com>
parents:
8960
diff
changeset
|
2429 { |
0a58c4cd1405
optimize element-wise power
Jaroslav Hajek <highegg@gmail.com>
parents:
8960
diff
changeset
|
2430 for (octave_idx_type i = 0; i < a.length (); i++) |
0a58c4cd1405
optimize element-wise power
Jaroslav Hajek <highegg@gmail.com>
parents:
8960
diff
changeset
|
2431 { |
0a58c4cd1405
optimize element-wise power
Jaroslav Hajek <highegg@gmail.com>
parents:
8960
diff
changeset
|
2432 OCTAVE_QUIT; |
0a58c4cd1405
optimize element-wise power
Jaroslav Hajek <highegg@gmail.com>
parents:
8960
diff
changeset
|
2433 result(i) = std::pow (a(i), ib); |
0a58c4cd1405
optimize element-wise power
Jaroslav Hajek <highegg@gmail.com>
parents:
8960
diff
changeset
|
2434 } |
0a58c4cd1405
optimize element-wise power
Jaroslav Hajek <highegg@gmail.com>
parents:
8960
diff
changeset
|
2435 } |
7789
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2436 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2437 retval = result; |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2438 } |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2439 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2440 return retval; |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2441 } |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2442 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2443 // -*- 4 -*- |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2444 octave_value |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2445 elem_xpow (const FloatNDArray& a, const FloatNDArray& b) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2446 { |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2447 octave_value retval; |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2448 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2449 dim_vector a_dims = a.dims (); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2450 dim_vector b_dims = b.dims (); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2451 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2452 if (a_dims != b_dims) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2453 { |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2454 gripe_nonconformant ("operator .^", a_dims, b_dims); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2455 return octave_value (); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2456 } |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2457 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2458 int len = a.length (); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2459 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2460 bool convert_to_complex = false; |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2461 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2462 for (octave_idx_type i = 0; i < len; i++) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2463 { |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2464 OCTAVE_QUIT; |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2465 float atmp = a(i); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2466 float btmp = b(i); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2467 if (atmp < 0.0 && static_cast<int> (btmp) != btmp) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2468 { |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2469 convert_to_complex = true; |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2470 goto done; |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2471 } |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2472 } |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2473 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2474 done: |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2475 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2476 if (convert_to_complex) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2477 { |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2478 FloatComplexNDArray complex_result (a_dims); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2479 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2480 for (octave_idx_type i = 0; i < len; i++) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2481 { |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2482 OCTAVE_QUIT; |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2483 FloatComplex atmp (a(i)); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2484 FloatComplex btmp (b(i)); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2485 complex_result(i) = std::pow (atmp, btmp); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2486 } |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2487 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2488 retval = complex_result; |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2489 } |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2490 else |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2491 { |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2492 FloatNDArray result (a_dims); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2493 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2494 for (octave_idx_type i = 0; i < len; i++) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2495 { |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2496 OCTAVE_QUIT; |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2497 result(i) = std::pow (a(i), b(i)); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2498 } |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2499 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2500 retval = result; |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2501 } |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2502 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2503 return retval; |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2504 } |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2505 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2506 // -*- 5 -*- |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2507 octave_value |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2508 elem_xpow (const FloatNDArray& a, const FloatComplex& b) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2509 { |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2510 FloatComplexNDArray result (a.dims ()); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2511 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2512 for (octave_idx_type i = 0; i < a.length (); i++) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2513 { |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2514 OCTAVE_QUIT; |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2515 result(i) = std::pow (FloatComplex (a(i)), b); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2516 } |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2517 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2518 return result; |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2519 } |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2520 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2521 // -*- 6 -*- |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2522 octave_value |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2523 elem_xpow (const FloatNDArray& a, const FloatComplexNDArray& b) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2524 { |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2525 dim_vector a_dims = a.dims (); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2526 dim_vector b_dims = b.dims (); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2527 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2528 if (a_dims != b_dims) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2529 { |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2530 gripe_nonconformant ("operator .^", a_dims, b_dims); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2531 return octave_value (); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2532 } |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2533 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2534 FloatComplexNDArray result (a_dims); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2535 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2536 for (octave_idx_type i = 0; i < a.length (); i++) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2537 { |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2538 OCTAVE_QUIT; |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2539 result(i) = std::pow (FloatComplex (a(i)), b(i)); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2540 } |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2541 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2542 return result; |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2543 } |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2544 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2545 // -*- 7 -*- |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2546 octave_value |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2547 elem_xpow (const FloatComplex& a, const FloatNDArray& b) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2548 { |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2549 FloatComplexNDArray result (b.dims ()); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2550 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2551 for (octave_idx_type i = 0; i < b.length (); i++) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2552 { |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2553 OCTAVE_QUIT; |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2554 float btmp = b(i); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2555 if (xisint (btmp)) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2556 result(i) = std::pow (a, static_cast<int> (btmp)); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2557 else |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2558 result(i) = std::pow (a, btmp); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2559 } |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2560 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2561 return result; |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2562 } |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2563 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2564 // -*- 8 -*- |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2565 octave_value |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2566 elem_xpow (const FloatComplex& a, const FloatComplexNDArray& b) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2567 { |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2568 FloatComplexNDArray result (b.dims ()); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2569 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2570 for (octave_idx_type i = 0; i < b.length (); i++) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2571 { |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2572 OCTAVE_QUIT; |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2573 result(i) = std::pow (a, b(i)); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2574 } |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2575 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2576 return result; |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2577 } |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2578 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2579 // -*- 9 -*- |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2580 octave_value |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2581 elem_xpow (const FloatComplexNDArray& a, float b) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2582 { |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2583 FloatComplexNDArray result (a.dims ()); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2584 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2585 if (xisint (b)) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2586 { |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2587 for (octave_idx_type i = 0; i < a.length (); i++) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2588 { |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2589 OCTAVE_QUIT; |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2590 result(i) = std::pow (a(i), static_cast<int> (b)); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2591 } |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2592 } |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2593 else |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2594 { |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2595 for (octave_idx_type i = 0; i < a.length (); i++) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2596 { |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2597 OCTAVE_QUIT; |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2598 result(i) = std::pow (a(i), b); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2599 } |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2600 } |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2601 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2602 return result; |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2603 } |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2604 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2605 // -*- 10 -*- |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2606 octave_value |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2607 elem_xpow (const FloatComplexNDArray& a, const FloatNDArray& b) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2608 { |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2609 dim_vector a_dims = a.dims (); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2610 dim_vector b_dims = b.dims (); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2611 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2612 if (a_dims != b_dims) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2613 { |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2614 gripe_nonconformant ("operator .^", a_dims, b_dims); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2615 return octave_value (); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2616 } |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2617 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2618 FloatComplexNDArray result (a_dims); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2619 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2620 for (octave_idx_type i = 0; i < a.length (); i++) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2621 { |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2622 OCTAVE_QUIT; |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2623 float btmp = b(i); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2624 if (xisint (btmp)) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2625 result(i) = std::pow (a(i), static_cast<int> (btmp)); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2626 else |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2627 result(i) = std::pow (a(i), btmp); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2628 } |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2629 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2630 return result; |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2631 } |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2632 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2633 // -*- 11 -*- |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2634 octave_value |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2635 elem_xpow (const FloatComplexNDArray& a, const FloatComplex& b) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2636 { |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2637 FloatComplexNDArray result (a.dims ()); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2638 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2639 for (octave_idx_type i = 0; i < a.length (); i++) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2640 { |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2641 OCTAVE_QUIT; |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2642 result(i) = std::pow (a(i), b); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2643 } |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2644 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2645 return result; |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2646 } |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2647 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2648 // -*- 12 -*- |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2649 octave_value |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2650 elem_xpow (const FloatComplexNDArray& a, const FloatComplexNDArray& b) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2651 { |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2652 dim_vector a_dims = a.dims (); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2653 dim_vector b_dims = b.dims (); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2654 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2655 if (a_dims != b_dims) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2656 { |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2657 gripe_nonconformant ("operator .^", a_dims, b_dims); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2658 return octave_value (); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2659 } |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2660 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2661 FloatComplexNDArray result (a_dims); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2662 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2663 for (octave_idx_type i = 0; i < a.length (); i++) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2664 { |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2665 OCTAVE_QUIT; |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2666 result(i) = std::pow (a(i), b(i)); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2667 } |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2668 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2669 return result; |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2670 } |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7578
diff
changeset
|
2671 |
1 | 2672 /* |
2673 ;;; Local Variables: *** | |
2674 ;;; mode: C++ *** | |
2675 ;;; End: *** | |
2676 */ |