|
1
|
1 /* |
|
|
2 |
|
2847
|
3 Copyright (C) 1996, 1997 John W. Eaton |
|
1
|
4 |
|
|
5 This file is part of Octave. |
|
|
6 |
|
|
7 Octave is free software; you can redistribute it and/or modify it |
|
|
8 under the terms of the GNU General Public License as published by the |
|
7016
|
9 Free Software Foundation; either version 3 of the License, or (at your |
|
|
10 option) any later version. |
|
1
|
11 |
|
|
12 Octave is distributed in the hope that it will be useful, but WITHOUT |
|
|
13 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or |
|
|
14 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License |
|
|
15 for more details. |
|
|
16 |
|
|
17 You should have received a copy of the GNU General Public License |
|
7016
|
18 along with Octave; see the file COPYING. If not, see |
|
|
19 <http://www.gnu.org/licenses/>. |
|
1
|
20 |
|
|
21 */ |
|
|
22 |
|
240
|
23 #ifdef HAVE_CONFIG_H |
|
1192
|
24 #include <config.h> |
|
1
|
25 #endif |
|
|
26 |
|
1343
|
27 #include <cassert> |
|
1580
|
28 #include <climits> |
|
1343
|
29 |
|
4669
|
30 #include "Array-util.h" |
|
1352
|
31 #include "CColVector.h" |
|
453
|
32 #include "CDiagMatrix.h" |
|
1352
|
33 #include "CMatrix.h" |
|
453
|
34 #include "EIG.h" |
|
1352
|
35 #include "dDiagMatrix.h" |
|
|
36 #include "dMatrix.h" |
|
3585
|
37 #include "mx-cm-cdm.h" |
|
1651
|
38 #include "oct-cmplx.h" |
|
4153
|
39 #include "quit.h" |
|
1352
|
40 |
|
|
41 #include "error.h" |
|
4055
|
42 #include "oct-obj.h" |
|
1567
|
43 #include "utils.h" |
|
1352
|
44 #include "xpow.h" |
|
1
|
45 |
|
5275
|
46 #ifdef _OPENMP |
|
|
47 #include <omp.h> |
|
|
48 #endif |
|
|
49 |
|
1567
|
50 static inline int |
|
|
51 xisint (double x) |
|
|
52 { |
|
|
53 return (D_NINT (x) == x |
|
|
54 && ((x >= 0 && x < INT_MAX) |
|
|
55 || (x <= 0 && x > INT_MIN))); |
|
|
56 } |
|
|
57 |
|
767
|
58 // Safer pow functions. |
|
|
59 // |
|
|
60 // op2 \ op1: s m cs cm |
|
|
61 // +-- +---+---+----+----+ |
|
|
62 // scalar | | 1 | 5 | 7 | 11 | |
|
|
63 // +---+---+----+----+ |
|
2365
|
64 // matrix | 2 | * | 8 | * | |
|
767
|
65 // +---+---+----+----+ |
|
|
66 // complex_scalar | 3 | 6 | 9 | 12 | |
|
|
67 // +---+---+----+----+ |
|
2365
|
68 // complex_matrix | 4 | * | 10 | * | |
|
767
|
69 // +---+---+----+----+ |
|
1
|
70 |
|
767
|
71 // -*- 1 -*- |
|
2086
|
72 octave_value |
|
1
|
73 xpow (double a, double b) |
|
|
74 { |
|
2800
|
75 if (a < 0.0 && static_cast<int> (b) != b) |
|
1
|
76 { |
|
5260
|
77 Complex atmp (a); |
|
4682
|
78 |
|
5260
|
79 return std::pow (atmp, b); |
|
1
|
80 } |
|
|
81 else |
|
5260
|
82 return std::pow (a, b); |
|
1
|
83 } |
|
|
84 |
|
767
|
85 // -*- 2 -*- |
|
2086
|
86 octave_value |
|
164
|
87 xpow (double a, const Matrix& b) |
|
1
|
88 { |
|
2086
|
89 octave_value retval; |
|
1
|
90 |
|
5275
|
91 octave_idx_type nr = b.rows (); |
|
|
92 octave_idx_type nc = b.cols (); |
|
1
|
93 |
|
|
94 if (nr == 0 || nc == 0 || nr != nc) |
|
|
95 error ("for x^A, A must be square"); |
|
|
96 else |
|
|
97 { |
|
|
98 EIG b_eig (b); |
|
|
99 ComplexColumnVector lambda (b_eig.eigenvalues ()); |
|
|
100 ComplexMatrix Q (b_eig.eigenvectors ()); |
|
|
101 |
|
5275
|
102 for (octave_idx_type i = 0; i < nr; i++) |
|
1
|
103 { |
|
2305
|
104 Complex elt = lambda (i); |
|
5260
|
105 if (std::imag (elt) == 0.0) |
|
|
106 lambda (i) = std::pow (a, std::real (elt)); |
|
1
|
107 else |
|
5260
|
108 lambda (i) = std::pow (a, elt); |
|
1
|
109 } |
|
|
110 ComplexDiagMatrix D (lambda); |
|
|
111 |
|
1567
|
112 retval = ComplexMatrix (Q * D * Q.inverse ()); |
|
1
|
113 } |
|
|
114 |
|
|
115 return retval; |
|
|
116 } |
|
|
117 |
|
767
|
118 // -*- 3 -*- |
|
2086
|
119 octave_value |
|
164
|
120 xpow (double a, const Complex& b) |
|
1
|
121 { |
|
|
122 Complex result; |
|
|
123 Complex atmp (a); |
|
5260
|
124 result = std::pow (atmp, b); |
|
1567
|
125 return result; |
|
1
|
126 } |
|
|
127 |
|
767
|
128 // -*- 4 -*- |
|
2086
|
129 octave_value |
|
164
|
130 xpow (double a, const ComplexMatrix& b) |
|
1
|
131 { |
|
2086
|
132 octave_value retval; |
|
1
|
133 |
|
5275
|
134 octave_idx_type nr = b.rows (); |
|
|
135 octave_idx_type nc = b.cols (); |
|
1
|
136 |
|
|
137 if (nr == 0 || nc == 0 || nr != nc) |
|
2365
|
138 error ("for x^A, A must be square"); |
|
1
|
139 else |
|
|
140 { |
|
|
141 EIG b_eig (b); |
|
|
142 ComplexColumnVector lambda (b_eig.eigenvalues ()); |
|
|
143 ComplexMatrix Q (b_eig.eigenvectors ()); |
|
|
144 |
|
5275
|
145 for (octave_idx_type i = 0; i < nr; i++) |
|
1
|
146 { |
|
2305
|
147 Complex elt = lambda (i); |
|
5260
|
148 if (std::imag (elt) == 0.0) |
|
|
149 lambda (i) = std::pow (a, std::real (elt)); |
|
1
|
150 else |
|
5260
|
151 lambda (i) = std::pow (a, elt); |
|
1
|
152 } |
|
|
153 ComplexDiagMatrix D (lambda); |
|
|
154 |
|
1567
|
155 retval = ComplexMatrix (Q * D * Q.inverse ()); |
|
1
|
156 } |
|
|
157 |
|
|
158 return retval; |
|
|
159 } |
|
|
160 |
|
767
|
161 // -*- 5 -*- |
|
2086
|
162 octave_value |
|
164
|
163 xpow (const Matrix& a, double b) |
|
1
|
164 { |
|
2086
|
165 octave_value retval; |
|
1
|
166 |
|
5275
|
167 octave_idx_type nr = a.rows (); |
|
|
168 octave_idx_type nc = a.cols (); |
|
1
|
169 |
|
|
170 if (nr == 0 || nc == 0 || nr != nc) |
|
2365
|
171 error ("for A^b, A must be square"); |
|
1567
|
172 else |
|
1
|
173 { |
|
2800
|
174 if (static_cast<int> (b) == b) |
|
1
|
175 { |
|
2804
|
176 int btmp = static_cast<int> (b); |
|
1567
|
177 if (btmp == 0) |
|
|
178 { |
|
|
179 retval = DiagMatrix (nr, nr, 1.0); |
|
|
180 } |
|
|
181 else |
|
|
182 { |
|
|
183 // Too much copying? |
|
5775
|
184 // FIXME -- we shouldn't do this if the exponent is |
|
1567
|
185 // large... |
|
|
186 |
|
|
187 Matrix atmp; |
|
|
188 if (btmp < 0) |
|
|
189 { |
|
|
190 btmp = -btmp; |
|
1655
|
191 |
|
5275
|
192 octave_idx_type info; |
|
1655
|
193 double rcond = 0.0; |
|
6207
|
194 MatrixType mattype (a); |
|
1655
|
195 |
|
6207
|
196 atmp = a.inverse (mattype, info, rcond, 1); |
|
1655
|
197 |
|
|
198 if (info == -1) |
|
|
199 warning ("inverse: matrix singular to machine\ |
|
|
200 precision, rcond = %g", rcond); |
|
1567
|
201 } |
|
|
202 else |
|
|
203 atmp = a; |
|
|
204 |
|
|
205 Matrix result (atmp); |
|
3178
|
206 |
|
|
207 btmp--; |
|
|
208 |
|
|
209 while (btmp > 0) |
|
|
210 { |
|
|
211 if (btmp & 1) |
|
|
212 result = result * atmp; |
|
|
213 |
|
|
214 btmp >>= 1; |
|
|
215 |
|
|
216 if (btmp > 0) |
|
|
217 atmp = atmp * atmp; |
|
|
218 } |
|
1567
|
219 |
|
|
220 retval = result; |
|
|
221 } |
|
1
|
222 } |
|
|
223 else |
|
|
224 { |
|
1567
|
225 EIG a_eig (a); |
|
|
226 ComplexColumnVector lambda (a_eig.eigenvalues ()); |
|
|
227 ComplexMatrix Q (a_eig.eigenvectors ()); |
|
|
228 |
|
5275
|
229 for (octave_idx_type i = 0; i < nr; i++) |
|
5260
|
230 lambda (i) = std::pow (lambda (i), b); |
|
1567
|
231 |
|
|
232 ComplexDiagMatrix D (lambda); |
|
|
233 |
|
|
234 retval = ComplexMatrix (Q * D * Q.inverse ()); |
|
|
235 } |
|
|
236 } |
|
1358
|
237 |
|
1567
|
238 return retval; |
|
|
239 } |
|
1
|
240 |
|
1567
|
241 // -*- 6 -*- |
|
2086
|
242 octave_value |
|
1567
|
243 xpow (const Matrix& a, const Complex& b) |
|
|
244 { |
|
2086
|
245 octave_value retval; |
|
1
|
246 |
|
5275
|
247 octave_idx_type nr = a.rows (); |
|
|
248 octave_idx_type nc = a.cols (); |
|
1567
|
249 |
|
|
250 if (nr == 0 || nc == 0 || nr != nc) |
|
2365
|
251 error ("for A^b, A must be square"); |
|
1
|
252 else |
|
|
253 { |
|
|
254 EIG a_eig (a); |
|
|
255 ComplexColumnVector lambda (a_eig.eigenvalues ()); |
|
|
256 ComplexMatrix Q (a_eig.eigenvectors ()); |
|
|
257 |
|
5275
|
258 for (octave_idx_type i = 0; i < nr; i++) |
|
5260
|
259 lambda (i) = std::pow (lambda (i), b); |
|
1
|
260 |
|
|
261 ComplexDiagMatrix D (lambda); |
|
|
262 |
|
1567
|
263 retval = ComplexMatrix (Q * D * Q.inverse ()); |
|
1
|
264 } |
|
|
265 |
|
|
266 return retval; |
|
|
267 } |
|
|
268 |
|
767
|
269 // -*- 7 -*- |
|
2086
|
270 octave_value |
|
164
|
271 xpow (const Complex& a, double b) |
|
1
|
272 { |
|
|
273 Complex result; |
|
1567
|
274 |
|
|
275 if (xisint (b)) |
|
5260
|
276 result = std::pow (a, static_cast<int> (b)); |
|
1567
|
277 else |
|
5260
|
278 result = std::pow (a, b); |
|
1567
|
279 |
|
|
280 return result; |
|
1
|
281 } |
|
|
282 |
|
767
|
283 // -*- 8 -*- |
|
2086
|
284 octave_value |
|
164
|
285 xpow (const Complex& a, const Matrix& b) |
|
1
|
286 { |
|
2086
|
287 octave_value retval; |
|
1
|
288 |
|
5275
|
289 octave_idx_type nr = b.rows (); |
|
|
290 octave_idx_type nc = b.cols (); |
|
1
|
291 |
|
|
292 if (nr == 0 || nc == 0 || nr != nc) |
|
2365
|
293 error ("for x^A, A must be square"); |
|
1
|
294 else |
|
|
295 { |
|
|
296 EIG b_eig (b); |
|
|
297 ComplexColumnVector lambda (b_eig.eigenvalues ()); |
|
|
298 ComplexMatrix Q (b_eig.eigenvectors ()); |
|
|
299 |
|
5275
|
300 for (octave_idx_type i = 0; i < nr; i++) |
|
1
|
301 { |
|
2305
|
302 Complex elt = lambda (i); |
|
5260
|
303 if (std::imag (elt) == 0.0) |
|
|
304 lambda (i) = std::pow (a, std::real (elt)); |
|
1
|
305 else |
|
5260
|
306 lambda (i) = std::pow (a, elt); |
|
1
|
307 } |
|
|
308 ComplexDiagMatrix D (lambda); |
|
|
309 |
|
1567
|
310 retval = ComplexMatrix (Q * D * Q.inverse ()); |
|
1
|
311 } |
|
|
312 |
|
|
313 return retval; |
|
|
314 } |
|
|
315 |
|
767
|
316 // -*- 9 -*- |
|
2086
|
317 octave_value |
|
164
|
318 xpow (const Complex& a, const Complex& b) |
|
1
|
319 { |
|
|
320 Complex result; |
|
5260
|
321 result = std::pow (a, b); |
|
1567
|
322 return result; |
|
1
|
323 } |
|
|
324 |
|
767
|
325 // -*- 10 -*- |
|
2086
|
326 octave_value |
|
164
|
327 xpow (const Complex& a, const ComplexMatrix& b) |
|
1
|
328 { |
|
2086
|
329 octave_value retval; |
|
1
|
330 |
|
5275
|
331 octave_idx_type nr = b.rows (); |
|
|
332 octave_idx_type nc = b.cols (); |
|
1
|
333 |
|
|
334 if (nr == 0 || nc == 0 || nr != nc) |
|
2365
|
335 error ("for x^A, A must be square"); |
|
1
|
336 else |
|
|
337 { |
|
|
338 EIG b_eig (b); |
|
|
339 ComplexColumnVector lambda (b_eig.eigenvalues ()); |
|
|
340 ComplexMatrix Q (b_eig.eigenvectors ()); |
|
|
341 |
|
5275
|
342 for (octave_idx_type i = 0; i < nr; i++) |
|
1
|
343 { |
|
2305
|
344 Complex elt = lambda (i); |
|
5260
|
345 if (std::imag (elt) == 0.0) |
|
|
346 lambda (i) = std::pow (a, std::real (elt)); |
|
1
|
347 else |
|
5260
|
348 lambda (i) = std::pow (a, elt); |
|
1
|
349 } |
|
|
350 ComplexDiagMatrix D (lambda); |
|
|
351 |
|
1567
|
352 retval = ComplexMatrix (Q * D * Q.inverse ()); |
|
1
|
353 } |
|
|
354 |
|
|
355 return retval; |
|
|
356 } |
|
|
357 |
|
767
|
358 // -*- 11 -*- |
|
2086
|
359 octave_value |
|
164
|
360 xpow (const ComplexMatrix& a, double b) |
|
1
|
361 { |
|
2086
|
362 octave_value retval; |
|
1
|
363 |
|
5275
|
364 octave_idx_type nr = a.rows (); |
|
|
365 octave_idx_type nc = a.cols (); |
|
1
|
366 |
|
|
367 if (nr == 0 || nc == 0 || nr != nc) |
|
2365
|
368 error ("for A^b, A must be square"); |
|
1567
|
369 else |
|
1
|
370 { |
|
2800
|
371 if (static_cast<int> (b) == b) |
|
1
|
372 { |
|
2804
|
373 int btmp = static_cast<int> (b); |
|
1567
|
374 if (btmp == 0) |
|
|
375 { |
|
|
376 retval = DiagMatrix (nr, nr, 1.0); |
|
|
377 } |
|
|
378 else |
|
|
379 { |
|
|
380 // Too much copying? |
|
5775
|
381 // FIXME -- we shouldn't do this if the exponent is |
|
1567
|
382 // large... |
|
|
383 |
|
|
384 ComplexMatrix atmp; |
|
|
385 if (btmp < 0) |
|
|
386 { |
|
|
387 btmp = -btmp; |
|
1655
|
388 |
|
5275
|
389 octave_idx_type info; |
|
1655
|
390 double rcond = 0.0; |
|
6207
|
391 MatrixType mattype (a); |
|
1655
|
392 |
|
6207
|
393 atmp = a.inverse (mattype, info, rcond, 1); |
|
1655
|
394 |
|
|
395 if (info == -1) |
|
|
396 warning ("inverse: matrix singular to machine\ |
|
|
397 precision, rcond = %g", rcond); |
|
1567
|
398 } |
|
|
399 else |
|
|
400 atmp = a; |
|
|
401 |
|
|
402 ComplexMatrix result (atmp); |
|
3178
|
403 |
|
|
404 btmp--; |
|
|
405 |
|
|
406 while (btmp > 0) |
|
|
407 { |
|
|
408 if (btmp & 1) |
|
|
409 result = result * atmp; |
|
|
410 |
|
|
411 btmp >>= 1; |
|
|
412 |
|
|
413 if (btmp > 0) |
|
|
414 atmp = atmp * atmp; |
|
|
415 } |
|
1567
|
416 |
|
|
417 retval = result; |
|
|
418 } |
|
1
|
419 } |
|
|
420 else |
|
|
421 { |
|
1567
|
422 EIG a_eig (a); |
|
|
423 ComplexColumnVector lambda (a_eig.eigenvalues ()); |
|
|
424 ComplexMatrix Q (a_eig.eigenvectors ()); |
|
|
425 |
|
5275
|
426 for (octave_idx_type i = 0; i < nr; i++) |
|
5260
|
427 lambda (i) = std::pow (lambda (i), b); |
|
1567
|
428 |
|
|
429 ComplexDiagMatrix D (lambda); |
|
|
430 |
|
|
431 retval = ComplexMatrix (Q * D * Q.inverse ()); |
|
|
432 } |
|
|
433 } |
|
1358
|
434 |
|
1567
|
435 return retval; |
|
|
436 } |
|
1
|
437 |
|
1567
|
438 // -*- 12 -*- |
|
2086
|
439 octave_value |
|
1567
|
440 xpow (const ComplexMatrix& a, const Complex& b) |
|
|
441 { |
|
2086
|
442 octave_value retval; |
|
1
|
443 |
|
5275
|
444 octave_idx_type nr = a.rows (); |
|
|
445 octave_idx_type nc = a.cols (); |
|
1567
|
446 |
|
|
447 if (nr == 0 || nc == 0 || nr != nc) |
|
2365
|
448 error ("for A^b, A must be square"); |
|
1
|
449 else |
|
|
450 { |
|
|
451 EIG a_eig (a); |
|
|
452 ComplexColumnVector lambda (a_eig.eigenvalues ()); |
|
|
453 ComplexMatrix Q (a_eig.eigenvectors ()); |
|
|
454 |
|
5275
|
455 for (octave_idx_type i = 0; i < nr; i++) |
|
5260
|
456 lambda (i) = std::pow (lambda (i), b); |
|
1
|
457 |
|
|
458 ComplexDiagMatrix D (lambda); |
|
|
459 |
|
1567
|
460 retval = ComplexMatrix (Q * D * Q.inverse ()); |
|
1
|
461 } |
|
|
462 |
|
|
463 return retval; |
|
|
464 } |
|
|
465 |
|
767
|
466 // Safer pow functions that work elementwise for matrices. |
|
|
467 // |
|
|
468 // op2 \ op1: s m cs cm |
|
|
469 // +-- +---+---+----+----+ |
|
|
470 // scalar | | * | 3 | * | 9 | |
|
|
471 // +---+---+----+----+ |
|
|
472 // matrix | 1 | 4 | 7 | 10 | |
|
|
473 // +---+---+----+----+ |
|
|
474 // complex_scalar | * | 5 | * | 11 | |
|
|
475 // +---+---+----+----+ |
|
|
476 // complex_matrix | 2 | 6 | 8 | 12 | |
|
|
477 // +---+---+----+----+ |
|
|
478 // |
|
|
479 // * -> not needed. |
|
1
|
480 |
|
5775
|
481 // FIXME -- these functions need to be fixed so that things |
|
3162
|
482 // like |
|
|
483 // |
|
|
484 // a = -1; b = [ 0, 0.5, 1 ]; r = a .^ b |
|
|
485 // |
|
|
486 // and |
|
|
487 // |
|
|
488 // a = -1; b = [ 0, 0.5, 1 ]; for i = 1:3, r(i) = a .^ b(i), end |
|
|
489 // |
|
|
490 // produce identical results. Also, it would be nice if -1^0.5 |
|
|
491 // produced a pure imaginary result instead of a complex number with a |
|
|
492 // small real part. But perhaps that's really a problem with the math |
|
|
493 // library... |
|
|
494 |
|
767
|
495 // -*- 1 -*- |
|
2086
|
496 octave_value |
|
164
|
497 elem_xpow (double a, const Matrix& b) |
|
1
|
498 { |
|
2086
|
499 octave_value retval; |
|
1
|
500 |
|
5275
|
501 octave_idx_type nr = b.rows (); |
|
|
502 octave_idx_type nc = b.cols (); |
|
1
|
503 |
|
3162
|
504 double d1, d2; |
|
1358
|
505 |
|
3162
|
506 if (a < 0.0 && ! b.all_integers (d1, d2)) |
|
1
|
507 { |
|
|
508 Complex atmp (a); |
|
|
509 ComplexMatrix result (nr, nc); |
|
5275
|
510 |
|
|
511 for (octave_idx_type j = 0; j < nc; j++) |
|
|
512 for (octave_idx_type i = 0; i < nr; i++) |
|
4153
|
513 { |
|
|
514 OCTAVE_QUIT; |
|
5260
|
515 result (i, j) = std::pow (atmp, b (i, j)); |
|
4153
|
516 } |
|
1
|
517 |
|
1567
|
518 retval = result; |
|
1
|
519 } |
|
|
520 else |
|
|
521 { |
|
|
522 Matrix result (nr, nc); |
|
5275
|
523 |
|
|
524 for (octave_idx_type j = 0; j < nc; j++) |
|
|
525 for (octave_idx_type i = 0; i < nr; i++) |
|
4153
|
526 { |
|
|
527 OCTAVE_QUIT; |
|
5260
|
528 result (i, j) = std::pow (a, b (i, j)); |
|
4153
|
529 } |
|
1
|
530 |
|
1567
|
531 retval = result; |
|
1
|
532 } |
|
|
533 |
|
|
534 return retval; |
|
|
535 } |
|
|
536 |
|
767
|
537 // -*- 2 -*- |
|
2086
|
538 octave_value |
|
164
|
539 elem_xpow (double a, const ComplexMatrix& b) |
|
1
|
540 { |
|
5275
|
541 octave_idx_type nr = b.rows (); |
|
|
542 octave_idx_type nc = b.cols (); |
|
1
|
543 |
|
|
544 ComplexMatrix result (nr, nc); |
|
3125
|
545 Complex atmp (a); |
|
5275
|
546 |
|
|
547 for (octave_idx_type j = 0; j < nc; j++) |
|
|
548 for (octave_idx_type i = 0; i < nr; i++) |
|
4153
|
549 { |
|
|
550 OCTAVE_QUIT; |
|
5260
|
551 result (i, j) = std::pow (atmp, b (i, j)); |
|
4153
|
552 } |
|
1
|
553 |
|
1567
|
554 return result; |
|
1
|
555 } |
|
|
556 |
|
767
|
557 // -*- 3 -*- |
|
2086
|
558 octave_value |
|
164
|
559 elem_xpow (const Matrix& a, double b) |
|
1
|
560 { |
|
2086
|
561 octave_value retval; |
|
1
|
562 |
|
5275
|
563 octave_idx_type nr = a.rows (); |
|
|
564 octave_idx_type nc = a.cols (); |
|
1
|
565 |
|
2800
|
566 if (static_cast<int> (b) != b && a.any_element_is_negative ()) |
|
1
|
567 { |
|
|
568 ComplexMatrix result (nr, nc); |
|
5275
|
569 |
|
|
570 for (octave_idx_type j = 0; j < nc; j++) |
|
|
571 for (octave_idx_type i = 0; i < nr; i++) |
|
1
|
572 { |
|
5665
|
573 OCTAVE_QUIT; |
|
|
574 |
|
|
575 Complex atmp (a (i, j)); |
|
|
576 |
|
|
577 result (i, j) = std::pow (atmp, b); |
|
1
|
578 } |
|
|
579 |
|
1567
|
580 retval = result; |
|
1
|
581 } |
|
|
582 else |
|
|
583 { |
|
|
584 Matrix result (nr, nc); |
|
5275
|
585 |
|
|
586 for (octave_idx_type j = 0; j < nc; j++) |
|
|
587 for (octave_idx_type i = 0; i < nr; i++) |
|
4153
|
588 { |
|
|
589 OCTAVE_QUIT; |
|
5260
|
590 result (i, j) = std::pow (a (i, j), b); |
|
4153
|
591 } |
|
1
|
592 |
|
1567
|
593 retval = result; |
|
1
|
594 } |
|
|
595 |
|
|
596 return retval; |
|
|
597 } |
|
|
598 |
|
767
|
599 // -*- 4 -*- |
|
2086
|
600 octave_value |
|
164
|
601 elem_xpow (const Matrix& a, const Matrix& b) |
|
1
|
602 { |
|
2086
|
603 octave_value retval; |
|
1567
|
604 |
|
5275
|
605 octave_idx_type nr = a.rows (); |
|
|
606 octave_idx_type nc = a.cols (); |
|
2365
|
607 |
|
5275
|
608 octave_idx_type b_nr = b.rows (); |
|
|
609 octave_idx_type b_nc = b.cols (); |
|
1
|
610 |
|
2365
|
611 if (nr != b_nr || nc != b_nc) |
|
|
612 { |
|
|
613 gripe_nonconformant ("operator .^", nr, nc, b_nr, b_nc); |
|
|
614 return octave_value (); |
|
|
615 } |
|
1
|
616 |
|
|
617 int convert_to_complex = 0; |
|
5275
|
618 for (octave_idx_type j = 0; j < nc; j++) |
|
|
619 for (octave_idx_type i = 0; i < nr; i++) |
|
1
|
620 { |
|
4153
|
621 OCTAVE_QUIT; |
|
2305
|
622 double atmp = a (i, j); |
|
|
623 double btmp = b (i, j); |
|
2800
|
624 if (atmp < 0.0 && static_cast<int> (btmp) != btmp) |
|
1
|
625 { |
|
|
626 convert_to_complex = 1; |
|
|
627 goto done; |
|
|
628 } |
|
|
629 } |
|
|
630 |
|
2365
|
631 done: |
|
1
|
632 |
|
|
633 if (convert_to_complex) |
|
|
634 { |
|
|
635 ComplexMatrix complex_result (nr, nc); |
|
|
636 |
|
5275
|
637 for (octave_idx_type j = 0; j < nc; j++) |
|
|
638 for (octave_idx_type i = 0; i < nr; i++) |
|
1
|
639 { |
|
4153
|
640 OCTAVE_QUIT; |
|
5665
|
641 Complex atmp (a (i, j)); |
|
|
642 Complex btmp (b (i, j)); |
|
|
643 complex_result (i, j) = std::pow (atmp, btmp); |
|
1
|
644 } |
|
1567
|
645 |
|
|
646 retval = complex_result; |
|
1
|
647 } |
|
|
648 else |
|
|
649 { |
|
|
650 Matrix result (nr, nc); |
|
|
651 |
|
5275
|
652 for (octave_idx_type j = 0; j < nc; j++) |
|
|
653 for (octave_idx_type i = 0; i < nr; i++) |
|
4153
|
654 { |
|
|
655 OCTAVE_QUIT; |
|
5260
|
656 result (i, j) = std::pow (a (i, j), b (i, j)); |
|
4153
|
657 } |
|
1
|
658 |
|
1567
|
659 retval = result; |
|
1
|
660 } |
|
1567
|
661 |
|
|
662 return retval; |
|
1
|
663 } |
|
|
664 |
|
767
|
665 // -*- 5 -*- |
|
2086
|
666 octave_value |
|
164
|
667 elem_xpow (const Matrix& a, const Complex& b) |
|
1
|
668 { |
|
5275
|
669 octave_idx_type nr = a.rows (); |
|
|
670 octave_idx_type nc = a.cols (); |
|
1
|
671 |
|
|
672 ComplexMatrix result (nr, nc); |
|
5275
|
673 |
|
|
674 for (octave_idx_type j = 0; j < nc; j++) |
|
|
675 for (octave_idx_type i = 0; i < nr; i++) |
|
4153
|
676 { |
|
|
677 OCTAVE_QUIT; |
|
5260
|
678 result (i, j) = std::pow (Complex (a (i, j)), b); |
|
4153
|
679 } |
|
1
|
680 |
|
1567
|
681 return result; |
|
1
|
682 } |
|
|
683 |
|
767
|
684 // -*- 6 -*- |
|
2086
|
685 octave_value |
|
164
|
686 elem_xpow (const Matrix& a, const ComplexMatrix& b) |
|
1
|
687 { |
|
5275
|
688 octave_idx_type nr = a.rows (); |
|
|
689 octave_idx_type nc = a.cols (); |
|
2365
|
690 |
|
5275
|
691 octave_idx_type b_nr = b.rows (); |
|
|
692 octave_idx_type b_nc = b.cols (); |
|
1
|
693 |
|
2365
|
694 if (nr != b_nr || nc != b_nc) |
|
|
695 { |
|
|
696 gripe_nonconformant ("operator .^", nr, nc, b_nr, b_nc); |
|
|
697 return octave_value (); |
|
|
698 } |
|
1
|
699 |
|
|
700 ComplexMatrix result (nr, nc); |
|
5275
|
701 |
|
|
702 for (octave_idx_type j = 0; j < nc; j++) |
|
|
703 for (octave_idx_type i = 0; i < nr; i++) |
|
4153
|
704 { |
|
|
705 OCTAVE_QUIT; |
|
5260
|
706 result (i, j) = std::pow (Complex (a (i, j)), b (i, j)); |
|
4153
|
707 } |
|
1
|
708 |
|
1567
|
709 return result; |
|
1
|
710 } |
|
|
711 |
|
767
|
712 // -*- 7 -*- |
|
2086
|
713 octave_value |
|
164
|
714 elem_xpow (const Complex& a, const Matrix& b) |
|
1
|
715 { |
|
5275
|
716 octave_idx_type nr = b.rows (); |
|
|
717 octave_idx_type nc = b.cols (); |
|
1
|
718 |
|
|
719 ComplexMatrix result (nr, nc); |
|
5275
|
720 |
|
|
721 for (octave_idx_type j = 0; j < nc; j++) |
|
|
722 for (octave_idx_type i = 0; i < nr; i++) |
|
1567
|
723 { |
|
4153
|
724 OCTAVE_QUIT; |
|
2305
|
725 double btmp = b (i, j); |
|
1567
|
726 if (xisint (btmp)) |
|
5260
|
727 result (i, j) = std::pow (a, static_cast<int> (btmp)); |
|
1567
|
728 else |
|
5260
|
729 result (i, j) = std::pow (a, btmp); |
|
1567
|
730 } |
|
1
|
731 |
|
1567
|
732 return result; |
|
1
|
733 } |
|
|
734 |
|
767
|
735 // -*- 8 -*- |
|
2086
|
736 octave_value |
|
164
|
737 elem_xpow (const Complex& a, const ComplexMatrix& b) |
|
1
|
738 { |
|
5275
|
739 octave_idx_type nr = b.rows (); |
|
|
740 octave_idx_type nc = b.cols (); |
|
1
|
741 |
|
|
742 ComplexMatrix result (nr, nc); |
|
5275
|
743 |
|
|
744 for (octave_idx_type j = 0; j < nc; j++) |
|
|
745 for (octave_idx_type i = 0; i < nr; i++) |
|
4153
|
746 { |
|
|
747 OCTAVE_QUIT; |
|
5260
|
748 result (i, j) = std::pow (a, b (i, j)); |
|
4153
|
749 } |
|
1
|
750 |
|
1567
|
751 return result; |
|
1
|
752 } |
|
|
753 |
|
767
|
754 // -*- 9 -*- |
|
2086
|
755 octave_value |
|
164
|
756 elem_xpow (const ComplexMatrix& a, double b) |
|
1
|
757 { |
|
5275
|
758 octave_idx_type nr = a.rows (); |
|
|
759 octave_idx_type nc = a.cols (); |
|
1
|
760 |
|
|
761 ComplexMatrix result (nr, nc); |
|
|
762 |
|
1567
|
763 if (xisint (b)) |
|
|
764 { |
|
5275
|
765 for (octave_idx_type j = 0; j < nc; j++) |
|
|
766 for (octave_idx_type i = 0; i < nr; i++) |
|
4153
|
767 { |
|
|
768 OCTAVE_QUIT; |
|
5260
|
769 result (i, j) = std::pow (a (i, j), static_cast<int> (b)); |
|
4153
|
770 } |
|
1567
|
771 } |
|
|
772 else |
|
|
773 { |
|
5275
|
774 for (octave_idx_type j = 0; j < nc; j++) |
|
|
775 for (octave_idx_type i = 0; i < nr; i++) |
|
4153
|
776 { |
|
|
777 OCTAVE_QUIT; |
|
5260
|
778 result (i, j) = std::pow (a (i, j), b); |
|
4153
|
779 } |
|
1567
|
780 } |
|
|
781 |
|
|
782 return result; |
|
1
|
783 } |
|
|
784 |
|
767
|
785 // -*- 10 -*- |
|
2086
|
786 octave_value |
|
164
|
787 elem_xpow (const ComplexMatrix& a, const Matrix& b) |
|
1
|
788 { |
|
5275
|
789 octave_idx_type nr = a.rows (); |
|
|
790 octave_idx_type nc = a.cols (); |
|
2365
|
791 |
|
5275
|
792 octave_idx_type b_nr = b.rows (); |
|
|
793 octave_idx_type b_nc = b.cols (); |
|
1
|
794 |
|
2365
|
795 if (nr != b_nr || nc != b_nc) |
|
|
796 { |
|
|
797 gripe_nonconformant ("operator .^", nr, nc, b_nr, b_nc); |
|
|
798 return octave_value (); |
|
|
799 } |
|
1
|
800 |
|
|
801 ComplexMatrix result (nr, nc); |
|
5275
|
802 |
|
|
803 for (octave_idx_type j = 0; j < nc; j++) |
|
|
804 for (octave_idx_type i = 0; i < nr; i++) |
|
1567
|
805 { |
|
4153
|
806 OCTAVE_QUIT; |
|
2305
|
807 double btmp = b (i, j); |
|
1567
|
808 if (xisint (btmp)) |
|
5260
|
809 result (i, j) = std::pow (a (i, j), static_cast<int> (btmp)); |
|
1567
|
810 else |
|
5260
|
811 result (i, j) = std::pow (a (i, j), btmp); |
|
1567
|
812 } |
|
1
|
813 |
|
1567
|
814 return result; |
|
1
|
815 } |
|
|
816 |
|
767
|
817 // -*- 11 -*- |
|
2086
|
818 octave_value |
|
164
|
819 elem_xpow (const ComplexMatrix& a, const Complex& b) |
|
1
|
820 { |
|
5275
|
821 octave_idx_type nr = a.rows (); |
|
|
822 octave_idx_type nc = a.cols (); |
|
1
|
823 |
|
|
824 ComplexMatrix result (nr, nc); |
|
5275
|
825 |
|
|
826 for (octave_idx_type j = 0; j < nc; j++) |
|
|
827 for (octave_idx_type i = 0; i < nr; i++) |
|
4153
|
828 { |
|
|
829 OCTAVE_QUIT; |
|
5260
|
830 result (i, j) = std::pow (a (i, j), b); |
|
4153
|
831 } |
|
1
|
832 |
|
1567
|
833 return result; |
|
1
|
834 } |
|
|
835 |
|
767
|
836 // -*- 12 -*- |
|
2086
|
837 octave_value |
|
164
|
838 elem_xpow (const ComplexMatrix& a, const ComplexMatrix& b) |
|
1
|
839 { |
|
5275
|
840 octave_idx_type nr = a.rows (); |
|
|
841 octave_idx_type nc = a.cols (); |
|
2365
|
842 |
|
5275
|
843 octave_idx_type b_nr = b.rows (); |
|
|
844 octave_idx_type b_nc = b.cols (); |
|
2365
|
845 |
|
|
846 if (nr != b_nr || nc != b_nc) |
|
|
847 { |
|
|
848 gripe_nonconformant ("operator .^", nr, nc, b_nr, b_nc); |
|
|
849 return octave_value (); |
|
|
850 } |
|
1
|
851 |
|
|
852 ComplexMatrix result (nr, nc); |
|
5275
|
853 |
|
|
854 for (octave_idx_type j = 0; j < nc; j++) |
|
|
855 for (octave_idx_type i = 0; i < nr; i++) |
|
4153
|
856 { |
|
|
857 OCTAVE_QUIT; |
|
5260
|
858 result (i, j) = std::pow (a (i, j), b (i, j)); |
|
4153
|
859 } |
|
1
|
860 |
|
1567
|
861 return result; |
|
1
|
862 } |
|
|
863 |
|
4543
|
864 // Safer pow functions that work elementwise for N-d arrays. |
|
|
865 // |
|
|
866 // op2 \ op1: s nd cs cnd |
|
|
867 // +-- +---+---+----+----+ |
|
|
868 // scalar | | * | 3 | * | 9 | |
|
|
869 // +---+---+----+----+ |
|
|
870 // N_d | 1 | 4 | 7 | 10 | |
|
|
871 // +---+---+----+----+ |
|
|
872 // complex_scalar | * | 5 | * | 11 | |
|
|
873 // +---+---+----+----+ |
|
|
874 // complex_N_d | 2 | 6 | 8 | 12 | |
|
|
875 // +---+---+----+----+ |
|
|
876 // |
|
|
877 // * -> not needed. |
|
|
878 |
|
5775
|
879 // FIXME -- these functions need to be fixed so that things |
|
4543
|
880 // like |
|
|
881 // |
|
|
882 // a = -1; b = [ 0, 0.5, 1 ]; r = a .^ b |
|
|
883 // |
|
|
884 // and |
|
|
885 // |
|
|
886 // a = -1; b = [ 0, 0.5, 1 ]; for i = 1:3, r(i) = a .^ b(i), end |
|
|
887 // |
|
|
888 // produce identical results. Also, it would be nice if -1^0.5 |
|
|
889 // produced a pure imaginary result instead of a complex number with a |
|
|
890 // small real part. But perhaps that's really a problem with the math |
|
|
891 // library... |
|
|
892 |
|
|
893 // -*- 1 -*- |
|
|
894 octave_value |
|
|
895 elem_xpow (double a, const NDArray& b) |
|
|
896 { |
|
|
897 octave_value retval; |
|
|
898 |
|
|
899 double d1, d2; |
|
|
900 |
|
|
901 if (a < 0.0 && ! b.all_integers (d1, d2)) |
|
|
902 { |
|
|
903 Complex atmp (a); |
|
|
904 ComplexNDArray result (b.dims ()); |
|
5275
|
905 for (octave_idx_type i = 0; i < b.length (); i++) |
|
4543
|
906 { |
|
|
907 OCTAVE_QUIT; |
|
5260
|
908 result(i) = std::pow (atmp, b(i)); |
|
4543
|
909 } |
|
|
910 |
|
|
911 retval = result; |
|
|
912 } |
|
|
913 else |
|
|
914 { |
|
|
915 NDArray result (b.dims ()); |
|
5275
|
916 for (octave_idx_type i = 0; i < b.length (); i++) |
|
4543
|
917 { |
|
|
918 OCTAVE_QUIT; |
|
5260
|
919 result (i) = std::pow (a, b(i)); |
|
4543
|
920 } |
|
|
921 |
|
|
922 retval = result; |
|
|
923 } |
|
|
924 |
|
|
925 return retval; |
|
|
926 } |
|
|
927 |
|
|
928 // -*- 2 -*- |
|
|
929 octave_value |
|
|
930 elem_xpow (double a, const ComplexNDArray& b) |
|
|
931 { |
|
|
932 ComplexNDArray result (b.dims ()); |
|
|
933 Complex atmp (a); |
|
5275
|
934 |
|
|
935 for (octave_idx_type i = 0; i < b.length (); i++) |
|
4543
|
936 { |
|
|
937 OCTAVE_QUIT; |
|
5260
|
938 result(i) = std::pow (atmp, b(i)); |
|
4543
|
939 } |
|
|
940 |
|
|
941 return result; |
|
|
942 } |
|
|
943 |
|
|
944 // -*- 3 -*- |
|
|
945 octave_value |
|
|
946 elem_xpow (const NDArray& a, double b) |
|
|
947 { |
|
|
948 octave_value retval; |
|
|
949 |
|
|
950 if (static_cast<int> (b) != b && a.any_element_is_negative ()) |
|
|
951 { |
|
|
952 ComplexNDArray result (a.dims ()); |
|
|
953 |
|
5275
|
954 for (octave_idx_type i = 0; i < a.length (); i++) |
|
4543
|
955 { |
|
|
956 OCTAVE_QUIT; |
|
5665
|
957 |
|
|
958 Complex atmp (a (i)); |
|
|
959 |
|
|
960 result(i) = std::pow (atmp, b); |
|
4543
|
961 } |
|
|
962 |
|
|
963 retval = result; |
|
|
964 } |
|
|
965 else |
|
|
966 { |
|
|
967 NDArray result (a.dims ()); |
|
|
968 |
|
5275
|
969 for (octave_idx_type i = 0; i < a.length (); i++) |
|
4543
|
970 { |
|
|
971 OCTAVE_QUIT; |
|
5260
|
972 result(i) = std::pow (a(i), b); |
|
4543
|
973 } |
|
|
974 |
|
|
975 retval = result; |
|
|
976 } |
|
|
977 |
|
|
978 return retval; |
|
|
979 } |
|
|
980 |
|
|
981 // -*- 4 -*- |
|
|
982 octave_value |
|
|
983 elem_xpow (const NDArray& a, const NDArray& b) |
|
|
984 { |
|
|
985 octave_value retval; |
|
|
986 |
|
|
987 dim_vector a_dims = a.dims (); |
|
|
988 dim_vector b_dims = b.dims (); |
|
|
989 |
|
|
990 if (a_dims != b_dims) |
|
|
991 { |
|
|
992 gripe_nonconformant ("operator .^", a_dims, b_dims); |
|
|
993 return octave_value (); |
|
|
994 } |
|
|
995 |
|
|
996 int len = a.length (); |
|
|
997 |
|
|
998 bool convert_to_complex = false; |
|
|
999 |
|
5275
|
1000 for (octave_idx_type i = 0; i < len; i++) |
|
4543
|
1001 { |
|
|
1002 OCTAVE_QUIT; |
|
|
1003 double atmp = a(i); |
|
|
1004 double btmp = b(i); |
|
|
1005 if (atmp < 0.0 && static_cast<int> (btmp) != btmp) |
|
|
1006 { |
|
|
1007 convert_to_complex = true; |
|
|
1008 goto done; |
|
|
1009 } |
|
|
1010 } |
|
|
1011 |
|
|
1012 done: |
|
|
1013 |
|
|
1014 if (convert_to_complex) |
|
|
1015 { |
|
|
1016 ComplexNDArray complex_result (a_dims); |
|
|
1017 |
|
5275
|
1018 for (octave_idx_type i = 0; i < len; i++) |
|
4543
|
1019 { |
|
|
1020 OCTAVE_QUIT; |
|
5665
|
1021 Complex atmp (a(i)); |
|
|
1022 Complex btmp (b(i)); |
|
|
1023 complex_result(i) = std::pow (atmp, btmp); |
|
4543
|
1024 } |
|
|
1025 |
|
|
1026 retval = complex_result; |
|
|
1027 } |
|
|
1028 else |
|
|
1029 { |
|
|
1030 NDArray result (a_dims); |
|
|
1031 |
|
5275
|
1032 for (octave_idx_type i = 0; i < len; i++) |
|
4543
|
1033 { |
|
|
1034 OCTAVE_QUIT; |
|
5260
|
1035 result(i) = std::pow (a(i), b(i)); |
|
4543
|
1036 } |
|
|
1037 |
|
|
1038 retval = result; |
|
|
1039 } |
|
|
1040 |
|
|
1041 return retval; |
|
|
1042 } |
|
|
1043 |
|
|
1044 // -*- 5 -*- |
|
|
1045 octave_value |
|
|
1046 elem_xpow (const NDArray& a, const Complex& b) |
|
|
1047 { |
|
|
1048 ComplexNDArray result (a.dims ()); |
|
|
1049 |
|
5275
|
1050 for (octave_idx_type i = 0; i < a.length (); i++) |
|
4543
|
1051 { |
|
|
1052 OCTAVE_QUIT; |
|
5260
|
1053 result(i) = std::pow (Complex (a(i)), b); |
|
4543
|
1054 } |
|
|
1055 |
|
|
1056 return result; |
|
|
1057 } |
|
|
1058 |
|
|
1059 // -*- 6 -*- |
|
|
1060 octave_value |
|
|
1061 elem_xpow (const NDArray& a, const ComplexNDArray& b) |
|
|
1062 { |
|
|
1063 dim_vector a_dims = a.dims (); |
|
|
1064 dim_vector b_dims = b.dims (); |
|
|
1065 |
|
|
1066 if (a_dims != b_dims) |
|
|
1067 { |
|
|
1068 gripe_nonconformant ("operator .^", a_dims, b_dims); |
|
|
1069 return octave_value (); |
|
|
1070 } |
|
|
1071 |
|
|
1072 ComplexNDArray result (a_dims); |
|
5275
|
1073 |
|
|
1074 for (octave_idx_type i = 0; i < a.length (); i++) |
|
4543
|
1075 { |
|
|
1076 OCTAVE_QUIT; |
|
5260
|
1077 result(i) = std::pow (Complex (a(i)), b(i)); |
|
4543
|
1078 } |
|
|
1079 |
|
|
1080 return result; |
|
|
1081 } |
|
|
1082 |
|
|
1083 // -*- 7 -*- |
|
|
1084 octave_value |
|
|
1085 elem_xpow (const Complex& a, const NDArray& b) |
|
|
1086 { |
|
|
1087 ComplexNDArray result (b.dims ()); |
|
5275
|
1088 |
|
|
1089 for (octave_idx_type i = 0; i < b.length (); i++) |
|
4543
|
1090 { |
|
|
1091 OCTAVE_QUIT; |
|
|
1092 double btmp = b(i); |
|
|
1093 if (xisint (btmp)) |
|
5260
|
1094 result(i) = std::pow (a, static_cast<int> (btmp)); |
|
4543
|
1095 else |
|
5260
|
1096 result(i) = std::pow (a, btmp); |
|
4543
|
1097 } |
|
|
1098 |
|
|
1099 return result; |
|
|
1100 } |
|
|
1101 |
|
|
1102 // -*- 8 -*- |
|
|
1103 octave_value |
|
|
1104 elem_xpow (const Complex& a, const ComplexNDArray& b) |
|
|
1105 { |
|
|
1106 ComplexNDArray result (b.dims ()); |
|
5275
|
1107 |
|
|
1108 for (octave_idx_type i = 0; i < b.length (); i++) |
|
4543
|
1109 { |
|
|
1110 OCTAVE_QUIT; |
|
5260
|
1111 result(i) = std::pow (a, b(i)); |
|
4543
|
1112 } |
|
|
1113 |
|
|
1114 return result; |
|
|
1115 } |
|
|
1116 |
|
|
1117 // -*- 9 -*- |
|
|
1118 octave_value |
|
|
1119 elem_xpow (const ComplexNDArray& a, double b) |
|
|
1120 { |
|
|
1121 ComplexNDArray result (a.dims ()); |
|
|
1122 |
|
|
1123 if (xisint (b)) |
|
|
1124 { |
|
5275
|
1125 for (octave_idx_type i = 0; i < a.length (); i++) |
|
4543
|
1126 { |
|
|
1127 OCTAVE_QUIT; |
|
5260
|
1128 result(i) = std::pow (a(i), static_cast<int> (b)); |
|
4543
|
1129 } |
|
|
1130 } |
|
|
1131 else |
|
|
1132 { |
|
5275
|
1133 for (octave_idx_type i = 0; i < a.length (); i++) |
|
4543
|
1134 { |
|
|
1135 OCTAVE_QUIT; |
|
5260
|
1136 result(i) = std::pow (a(i), b); |
|
4543
|
1137 } |
|
|
1138 } |
|
|
1139 |
|
|
1140 return result; |
|
|
1141 } |
|
|
1142 |
|
|
1143 // -*- 10 -*- |
|
|
1144 octave_value |
|
|
1145 elem_xpow (const ComplexNDArray& a, const NDArray& b) |
|
|
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 ComplexNDArray result (a_dims); |
|
5275
|
1157 |
|
|
1158 for (octave_idx_type i = 0; i < a.length (); i++) |
|
4543
|
1159 { |
|
|
1160 OCTAVE_QUIT; |
|
|
1161 double btmp = b(i); |
|
|
1162 if (xisint (btmp)) |
|
5260
|
1163 result(i) = std::pow (a(i), static_cast<int> (btmp)); |
|
4543
|
1164 else |
|
5260
|
1165 result(i) = std::pow (a(i), btmp); |
|
4543
|
1166 } |
|
|
1167 |
|
|
1168 return result; |
|
|
1169 } |
|
|
1170 |
|
|
1171 // -*- 11 -*- |
|
|
1172 octave_value |
|
|
1173 elem_xpow (const ComplexNDArray& a, const Complex& b) |
|
|
1174 { |
|
|
1175 ComplexNDArray result (a.dims ()); |
|
5275
|
1176 |
|
|
1177 for (octave_idx_type i = 0; i < a.length (); i++) |
|
4543
|
1178 { |
|
|
1179 OCTAVE_QUIT; |
|
5260
|
1180 result(i) = std::pow (a(i), b); |
|
4543
|
1181 } |
|
|
1182 |
|
|
1183 return result; |
|
|
1184 } |
|
|
1185 |
|
|
1186 // -*- 12 -*- |
|
|
1187 octave_value |
|
|
1188 elem_xpow (const ComplexNDArray& a, const ComplexNDArray& b) |
|
|
1189 { |
|
|
1190 dim_vector a_dims = a.dims (); |
|
|
1191 dim_vector b_dims = b.dims (); |
|
|
1192 |
|
|
1193 if (a_dims != b_dims) |
|
|
1194 { |
|
|
1195 gripe_nonconformant ("operator .^", a_dims, b_dims); |
|
|
1196 return octave_value (); |
|
|
1197 } |
|
|
1198 |
|
|
1199 ComplexNDArray result (a_dims); |
|
5275
|
1200 |
|
|
1201 for (octave_idx_type i = 0; i < a.length (); i++) |
|
4543
|
1202 { |
|
|
1203 OCTAVE_QUIT; |
|
5260
|
1204 result(i) = std::pow (a(i), b(i)); |
|
4543
|
1205 } |
|
|
1206 |
|
|
1207 return result; |
|
|
1208 } |
|
|
1209 |
|
1
|
1210 /* |
|
|
1211 ;;; Local Variables: *** |
|
|
1212 ;;; mode: C++ *** |
|
|
1213 ;;; End: *** |
|
|
1214 */ |
