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