|
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? |
|
|
185 // XXX FIXME XXX -- we shouldn't do this if the exponent is |
|
|
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; |
|
|
195 |
|
1656
|
196 atmp = a.inverse (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? |
|
|
381 // XXX FIXME XXX -- we shouldn't do this if the exponent is |
|
|
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; |
|
|
391 |
|
1656
|
392 atmp = a.inverse (info, rcond, 1); |
|
1655
|
393 |
|
|
394 if (info == -1) |
|
|
395 warning ("inverse: matrix singular to machine\ |
|
|
396 precision, rcond = %g", rcond); |
|
1567
|
397 } |
|
|
398 else |
|
|
399 atmp = a; |
|
|
400 |
|
|
401 ComplexMatrix result (atmp); |
|
3178
|
402 |
|
|
403 btmp--; |
|
|
404 |
|
|
405 while (btmp > 0) |
|
|
406 { |
|
|
407 if (btmp & 1) |
|
|
408 result = result * atmp; |
|
|
409 |
|
|
410 btmp >>= 1; |
|
|
411 |
|
|
412 if (btmp > 0) |
|
|
413 atmp = atmp * atmp; |
|
|
414 } |
|
1567
|
415 |
|
|
416 retval = result; |
|
|
417 } |
|
1
|
418 } |
|
|
419 else |
|
|
420 { |
|
1567
|
421 EIG a_eig (a); |
|
|
422 ComplexColumnVector lambda (a_eig.eigenvalues ()); |
|
|
423 ComplexMatrix Q (a_eig.eigenvectors ()); |
|
|
424 |
|
5275
|
425 for (octave_idx_type i = 0; i < nr; i++) |
|
5260
|
426 lambda (i) = std::pow (lambda (i), b); |
|
1567
|
427 |
|
|
428 ComplexDiagMatrix D (lambda); |
|
|
429 |
|
|
430 retval = ComplexMatrix (Q * D * Q.inverse ()); |
|
|
431 } |
|
|
432 } |
|
1358
|
433 |
|
1567
|
434 return retval; |
|
|
435 } |
|
1
|
436 |
|
1567
|
437 // -*- 12 -*- |
|
2086
|
438 octave_value |
|
1567
|
439 xpow (const ComplexMatrix& a, const Complex& b) |
|
|
440 { |
|
2086
|
441 octave_value retval; |
|
1
|
442 |
|
5275
|
443 octave_idx_type nr = a.rows (); |
|
|
444 octave_idx_type nc = a.cols (); |
|
1567
|
445 |
|
|
446 if (nr == 0 || nc == 0 || nr != nc) |
|
2365
|
447 error ("for A^b, A must be square"); |
|
1
|
448 else |
|
|
449 { |
|
|
450 EIG a_eig (a); |
|
|
451 ComplexColumnVector lambda (a_eig.eigenvalues ()); |
|
|
452 ComplexMatrix Q (a_eig.eigenvectors ()); |
|
|
453 |
|
5275
|
454 for (octave_idx_type i = 0; i < nr; i++) |
|
5260
|
455 lambda (i) = std::pow (lambda (i), b); |
|
1
|
456 |
|
|
457 ComplexDiagMatrix D (lambda); |
|
|
458 |
|
1567
|
459 retval = ComplexMatrix (Q * D * Q.inverse ()); |
|
1
|
460 } |
|
|
461 |
|
|
462 return retval; |
|
|
463 } |
|
|
464 |
|
767
|
465 // Safer pow functions that work elementwise for matrices. |
|
|
466 // |
|
|
467 // op2 \ op1: s m cs cm |
|
|
468 // +-- +---+---+----+----+ |
|
|
469 // scalar | | * | 3 | * | 9 | |
|
|
470 // +---+---+----+----+ |
|
|
471 // matrix | 1 | 4 | 7 | 10 | |
|
|
472 // +---+---+----+----+ |
|
|
473 // complex_scalar | * | 5 | * | 11 | |
|
|
474 // +---+---+----+----+ |
|
|
475 // complex_matrix | 2 | 6 | 8 | 12 | |
|
|
476 // +---+---+----+----+ |
|
|
477 // |
|
|
478 // * -> not needed. |
|
1
|
479 |
|
3162
|
480 // XXX FIXME XXX -- these functions need to be fixed so that things |
|
|
481 // like |
|
|
482 // |
|
|
483 // a = -1; b = [ 0, 0.5, 1 ]; r = a .^ b |
|
|
484 // |
|
|
485 // and |
|
|
486 // |
|
|
487 // a = -1; b = [ 0, 0.5, 1 ]; for i = 1:3, r(i) = a .^ b(i), end |
|
|
488 // |
|
|
489 // produce identical results. Also, it would be nice if -1^0.5 |
|
|
490 // produced a pure imaginary result instead of a complex number with a |
|
|
491 // small real part. But perhaps that's really a problem with the math |
|
|
492 // library... |
|
|
493 |
|
767
|
494 // -*- 1 -*- |
|
2086
|
495 octave_value |
|
164
|
496 elem_xpow (double a, const Matrix& b) |
|
1
|
497 { |
|
2086
|
498 octave_value retval; |
|
1
|
499 |
|
5275
|
500 octave_idx_type nr = b.rows (); |
|
|
501 octave_idx_type nc = b.cols (); |
|
1
|
502 |
|
3162
|
503 double d1, d2; |
|
1358
|
504 |
|
3162
|
505 if (a < 0.0 && ! b.all_integers (d1, d2)) |
|
1
|
506 { |
|
|
507 Complex atmp (a); |
|
|
508 ComplexMatrix result (nr, nc); |
|
5275
|
509 |
|
|
510 for (octave_idx_type j = 0; j < nc; j++) |
|
|
511 for (octave_idx_type i = 0; i < nr; i++) |
|
4153
|
512 { |
|
|
513 OCTAVE_QUIT; |
|
5260
|
514 result (i, j) = std::pow (atmp, b (i, j)); |
|
4153
|
515 } |
|
1
|
516 |
|
1567
|
517 retval = result; |
|
1
|
518 } |
|
|
519 else |
|
|
520 { |
|
|
521 Matrix result (nr, nc); |
|
5275
|
522 |
|
|
523 for (octave_idx_type j = 0; j < nc; j++) |
|
|
524 for (octave_idx_type i = 0; i < nr; i++) |
|
4153
|
525 { |
|
|
526 OCTAVE_QUIT; |
|
5260
|
527 result (i, j) = std::pow (a, b (i, j)); |
|
4153
|
528 } |
|
1
|
529 |
|
1567
|
530 retval = result; |
|
1
|
531 } |
|
|
532 |
|
|
533 return retval; |
|
|
534 } |
|
|
535 |
|
767
|
536 // -*- 2 -*- |
|
2086
|
537 octave_value |
|
164
|
538 elem_xpow (double a, const ComplexMatrix& b) |
|
1
|
539 { |
|
5275
|
540 octave_idx_type nr = b.rows (); |
|
|
541 octave_idx_type nc = b.cols (); |
|
1
|
542 |
|
|
543 ComplexMatrix result (nr, nc); |
|
3125
|
544 Complex atmp (a); |
|
5275
|
545 |
|
|
546 for (octave_idx_type j = 0; j < nc; j++) |
|
|
547 for (octave_idx_type i = 0; i < nr; i++) |
|
4153
|
548 { |
|
|
549 OCTAVE_QUIT; |
|
5260
|
550 result (i, j) = std::pow (atmp, b (i, j)); |
|
4153
|
551 } |
|
1
|
552 |
|
1567
|
553 return result; |
|
1
|
554 } |
|
|
555 |
|
767
|
556 // -*- 3 -*- |
|
2086
|
557 octave_value |
|
164
|
558 elem_xpow (const Matrix& a, double b) |
|
1
|
559 { |
|
2086
|
560 octave_value retval; |
|
1
|
561 |
|
5275
|
562 octave_idx_type nr = a.rows (); |
|
|
563 octave_idx_type nc = a.cols (); |
|
1
|
564 |
|
2800
|
565 if (static_cast<int> (b) != b && a.any_element_is_negative ()) |
|
1
|
566 { |
|
|
567 ComplexMatrix result (nr, nc); |
|
5275
|
568 |
|
|
569 for (octave_idx_type j = 0; j < nc; j++) |
|
|
570 for (octave_idx_type i = 0; i < nr; i++) |
|
1
|
571 { |
|
5665
|
572 OCTAVE_QUIT; |
|
|
573 |
|
|
574 Complex atmp (a (i, j)); |
|
|
575 |
|
|
576 result (i, j) = std::pow (atmp, b); |
|
1
|
577 } |
|
|
578 |
|
1567
|
579 retval = result; |
|
1
|
580 } |
|
|
581 else |
|
|
582 { |
|
|
583 Matrix result (nr, nc); |
|
5275
|
584 |
|
|
585 for (octave_idx_type j = 0; j < nc; j++) |
|
|
586 for (octave_idx_type i = 0; i < nr; i++) |
|
4153
|
587 { |
|
|
588 OCTAVE_QUIT; |
|
5260
|
589 result (i, j) = std::pow (a (i, j), b); |
|
4153
|
590 } |
|
1
|
591 |
|
1567
|
592 retval = result; |
|
1
|
593 } |
|
|
594 |
|
|
595 return retval; |
|
|
596 } |
|
|
597 |
|
767
|
598 // -*- 4 -*- |
|
2086
|
599 octave_value |
|
164
|
600 elem_xpow (const Matrix& a, const Matrix& b) |
|
1
|
601 { |
|
2086
|
602 octave_value retval; |
|
1567
|
603 |
|
5275
|
604 octave_idx_type nr = a.rows (); |
|
|
605 octave_idx_type nc = a.cols (); |
|
2365
|
606 |
|
5275
|
607 octave_idx_type b_nr = b.rows (); |
|
|
608 octave_idx_type b_nc = b.cols (); |
|
1
|
609 |
|
2365
|
610 if (nr != b_nr || nc != b_nc) |
|
|
611 { |
|
|
612 gripe_nonconformant ("operator .^", nr, nc, b_nr, b_nc); |
|
|
613 return octave_value (); |
|
|
614 } |
|
1
|
615 |
|
|
616 int convert_to_complex = 0; |
|
5275
|
617 for (octave_idx_type j = 0; j < nc; j++) |
|
|
618 for (octave_idx_type i = 0; i < nr; i++) |
|
1
|
619 { |
|
4153
|
620 OCTAVE_QUIT; |
|
2305
|
621 double atmp = a (i, j); |
|
|
622 double btmp = b (i, j); |
|
2800
|
623 if (atmp < 0.0 && static_cast<int> (btmp) != btmp) |
|
1
|
624 { |
|
|
625 convert_to_complex = 1; |
|
|
626 goto done; |
|
|
627 } |
|
|
628 } |
|
|
629 |
|
2365
|
630 done: |
|
1
|
631 |
|
|
632 if (convert_to_complex) |
|
|
633 { |
|
|
634 ComplexMatrix complex_result (nr, nc); |
|
|
635 |
|
5275
|
636 for (octave_idx_type j = 0; j < nc; j++) |
|
|
637 for (octave_idx_type i = 0; i < nr; i++) |
|
1
|
638 { |
|
4153
|
639 OCTAVE_QUIT; |
|
5665
|
640 Complex atmp (a (i, j)); |
|
|
641 Complex btmp (b (i, j)); |
|
|
642 complex_result (i, j) = std::pow (atmp, btmp); |
|
1
|
643 } |
|
1567
|
644 |
|
|
645 retval = complex_result; |
|
1
|
646 } |
|
|
647 else |
|
|
648 { |
|
|
649 Matrix result (nr, nc); |
|
|
650 |
|
5275
|
651 for (octave_idx_type j = 0; j < nc; j++) |
|
|
652 for (octave_idx_type i = 0; i < nr; i++) |
|
4153
|
653 { |
|
|
654 OCTAVE_QUIT; |
|
5260
|
655 result (i, j) = std::pow (a (i, j), b (i, j)); |
|
4153
|
656 } |
|
1
|
657 |
|
1567
|
658 retval = result; |
|
1
|
659 } |
|
1567
|
660 |
|
|
661 return retval; |
|
1
|
662 } |
|
|
663 |
|
767
|
664 // -*- 5 -*- |
|
2086
|
665 octave_value |
|
164
|
666 elem_xpow (const Matrix& a, const Complex& b) |
|
1
|
667 { |
|
5275
|
668 octave_idx_type nr = a.rows (); |
|
|
669 octave_idx_type nc = a.cols (); |
|
1
|
670 |
|
|
671 ComplexMatrix result (nr, nc); |
|
5275
|
672 |
|
|
673 for (octave_idx_type j = 0; j < nc; j++) |
|
|
674 for (octave_idx_type i = 0; i < nr; i++) |
|
4153
|
675 { |
|
|
676 OCTAVE_QUIT; |
|
5260
|
677 result (i, j) = std::pow (Complex (a (i, j)), b); |
|
4153
|
678 } |
|
1
|
679 |
|
1567
|
680 return result; |
|
1
|
681 } |
|
|
682 |
|
767
|
683 // -*- 6 -*- |
|
2086
|
684 octave_value |
|
164
|
685 elem_xpow (const Matrix& a, const ComplexMatrix& b) |
|
1
|
686 { |
|
5275
|
687 octave_idx_type nr = a.rows (); |
|
|
688 octave_idx_type nc = a.cols (); |
|
2365
|
689 |
|
5275
|
690 octave_idx_type b_nr = b.rows (); |
|
|
691 octave_idx_type b_nc = b.cols (); |
|
1
|
692 |
|
2365
|
693 if (nr != b_nr || nc != b_nc) |
|
|
694 { |
|
|
695 gripe_nonconformant ("operator .^", nr, nc, b_nr, b_nc); |
|
|
696 return octave_value (); |
|
|
697 } |
|
1
|
698 |
|
|
699 ComplexMatrix result (nr, nc); |
|
5275
|
700 |
|
|
701 for (octave_idx_type j = 0; j < nc; j++) |
|
|
702 for (octave_idx_type i = 0; i < nr; i++) |
|
4153
|
703 { |
|
|
704 OCTAVE_QUIT; |
|
5260
|
705 result (i, j) = std::pow (Complex (a (i, j)), b (i, j)); |
|
4153
|
706 } |
|
1
|
707 |
|
1567
|
708 return result; |
|
1
|
709 } |
|
|
710 |
|
767
|
711 // -*- 7 -*- |
|
2086
|
712 octave_value |
|
164
|
713 elem_xpow (const Complex& a, const Matrix& b) |
|
1
|
714 { |
|
5275
|
715 octave_idx_type nr = b.rows (); |
|
|
716 octave_idx_type nc = b.cols (); |
|
1
|
717 |
|
|
718 ComplexMatrix result (nr, nc); |
|
5275
|
719 |
|
|
720 for (octave_idx_type j = 0; j < nc; j++) |
|
|
721 for (octave_idx_type i = 0; i < nr; i++) |
|
1567
|
722 { |
|
4153
|
723 OCTAVE_QUIT; |
|
2305
|
724 double btmp = b (i, j); |
|
1567
|
725 if (xisint (btmp)) |
|
5260
|
726 result (i, j) = std::pow (a, static_cast<int> (btmp)); |
|
1567
|
727 else |
|
5260
|
728 result (i, j) = std::pow (a, btmp); |
|
1567
|
729 } |
|
1
|
730 |
|
1567
|
731 return result; |
|
1
|
732 } |
|
|
733 |
|
767
|
734 // -*- 8 -*- |
|
2086
|
735 octave_value |
|
164
|
736 elem_xpow (const Complex& a, const ComplexMatrix& b) |
|
1
|
737 { |
|
5275
|
738 octave_idx_type nr = b.rows (); |
|
|
739 octave_idx_type nc = b.cols (); |
|
1
|
740 |
|
|
741 ComplexMatrix result (nr, nc); |
|
5275
|
742 |
|
|
743 for (octave_idx_type j = 0; j < nc; j++) |
|
|
744 for (octave_idx_type i = 0; i < nr; i++) |
|
4153
|
745 { |
|
|
746 OCTAVE_QUIT; |
|
5260
|
747 result (i, j) = std::pow (a, b (i, j)); |
|
4153
|
748 } |
|
1
|
749 |
|
1567
|
750 return result; |
|
1
|
751 } |
|
|
752 |
|
767
|
753 // -*- 9 -*- |
|
2086
|
754 octave_value |
|
164
|
755 elem_xpow (const ComplexMatrix& a, double b) |
|
1
|
756 { |
|
5275
|
757 octave_idx_type nr = a.rows (); |
|
|
758 octave_idx_type nc = a.cols (); |
|
1
|
759 |
|
|
760 ComplexMatrix result (nr, nc); |
|
|
761 |
|
1567
|
762 if (xisint (b)) |
|
|
763 { |
|
5275
|
764 for (octave_idx_type j = 0; j < nc; j++) |
|
|
765 for (octave_idx_type i = 0; i < nr; i++) |
|
4153
|
766 { |
|
|
767 OCTAVE_QUIT; |
|
5260
|
768 result (i, j) = std::pow (a (i, j), static_cast<int> (b)); |
|
4153
|
769 } |
|
1567
|
770 } |
|
|
771 else |
|
|
772 { |
|
5275
|
773 for (octave_idx_type j = 0; j < nc; j++) |
|
|
774 for (octave_idx_type i = 0; i < nr; i++) |
|
4153
|
775 { |
|
|
776 OCTAVE_QUIT; |
|
5260
|
777 result (i, j) = std::pow (a (i, j), b); |
|
4153
|
778 } |
|
1567
|
779 } |
|
|
780 |
|
|
781 return result; |
|
1
|
782 } |
|
|
783 |
|
767
|
784 // -*- 10 -*- |
|
2086
|
785 octave_value |
|
164
|
786 elem_xpow (const ComplexMatrix& a, const Matrix& b) |
|
1
|
787 { |
|
5275
|
788 octave_idx_type nr = a.rows (); |
|
|
789 octave_idx_type nc = a.cols (); |
|
2365
|
790 |
|
5275
|
791 octave_idx_type b_nr = b.rows (); |
|
|
792 octave_idx_type b_nc = b.cols (); |
|
1
|
793 |
|
2365
|
794 if (nr != b_nr || nc != b_nc) |
|
|
795 { |
|
|
796 gripe_nonconformant ("operator .^", nr, nc, b_nr, b_nc); |
|
|
797 return octave_value (); |
|
|
798 } |
|
1
|
799 |
|
|
800 ComplexMatrix result (nr, nc); |
|
5275
|
801 |
|
|
802 for (octave_idx_type j = 0; j < nc; j++) |
|
|
803 for (octave_idx_type i = 0; i < nr; i++) |
|
1567
|
804 { |
|
4153
|
805 OCTAVE_QUIT; |
|
2305
|
806 double btmp = b (i, j); |
|
1567
|
807 if (xisint (btmp)) |
|
5260
|
808 result (i, j) = std::pow (a (i, j), static_cast<int> (btmp)); |
|
1567
|
809 else |
|
5260
|
810 result (i, j) = std::pow (a (i, j), btmp); |
|
1567
|
811 } |
|
1
|
812 |
|
1567
|
813 return result; |
|
1
|
814 } |
|
|
815 |
|
767
|
816 // -*- 11 -*- |
|
2086
|
817 octave_value |
|
164
|
818 elem_xpow (const ComplexMatrix& a, const Complex& b) |
|
1
|
819 { |
|
5275
|
820 octave_idx_type nr = a.rows (); |
|
|
821 octave_idx_type nc = a.cols (); |
|
1
|
822 |
|
|
823 ComplexMatrix result (nr, nc); |
|
5275
|
824 |
|
|
825 for (octave_idx_type j = 0; j < nc; j++) |
|
|
826 for (octave_idx_type i = 0; i < nr; i++) |
|
4153
|
827 { |
|
|
828 OCTAVE_QUIT; |
|
5260
|
829 result (i, j) = std::pow (a (i, j), b); |
|
4153
|
830 } |
|
1
|
831 |
|
1567
|
832 return result; |
|
1
|
833 } |
|
|
834 |
|
767
|
835 // -*- 12 -*- |
|
2086
|
836 octave_value |
|
164
|
837 elem_xpow (const ComplexMatrix& a, const ComplexMatrix& b) |
|
1
|
838 { |
|
5275
|
839 octave_idx_type nr = a.rows (); |
|
|
840 octave_idx_type nc = a.cols (); |
|
2365
|
841 |
|
5275
|
842 octave_idx_type b_nr = b.rows (); |
|
|
843 octave_idx_type b_nc = b.cols (); |
|
2365
|
844 |
|
|
845 if (nr != b_nr || nc != b_nc) |
|
|
846 { |
|
|
847 gripe_nonconformant ("operator .^", nr, nc, b_nr, b_nc); |
|
|
848 return octave_value (); |
|
|
849 } |
|
1
|
850 |
|
|
851 ComplexMatrix result (nr, nc); |
|
5275
|
852 |
|
|
853 for (octave_idx_type j = 0; j < nc; j++) |
|
|
854 for (octave_idx_type i = 0; i < nr; i++) |
|
4153
|
855 { |
|
|
856 OCTAVE_QUIT; |
|
5260
|
857 result (i, j) = std::pow (a (i, j), b (i, j)); |
|
4153
|
858 } |
|
1
|
859 |
|
1567
|
860 return result; |
|
1
|
861 } |
|
|
862 |
|
4543
|
863 // Safer pow functions that work elementwise for N-d arrays. |
|
|
864 // |
|
|
865 // op2 \ op1: s nd cs cnd |
|
|
866 // +-- +---+---+----+----+ |
|
|
867 // scalar | | * | 3 | * | 9 | |
|
|
868 // +---+---+----+----+ |
|
|
869 // N_d | 1 | 4 | 7 | 10 | |
|
|
870 // +---+---+----+----+ |
|
|
871 // complex_scalar | * | 5 | * | 11 | |
|
|
872 // +---+---+----+----+ |
|
|
873 // complex_N_d | 2 | 6 | 8 | 12 | |
|
|
874 // +---+---+----+----+ |
|
|
875 // |
|
|
876 // * -> not needed. |
|
|
877 |
|
|
878 // XXX FIXME XXX -- these functions need to be fixed so that things |
|
|
879 // like |
|
|
880 // |
|
|
881 // a = -1; b = [ 0, 0.5, 1 ]; r = a .^ b |
|
|
882 // |
|
|
883 // and |
|
|
884 // |
|
|
885 // a = -1; b = [ 0, 0.5, 1 ]; for i = 1:3, r(i) = a .^ b(i), end |
|
|
886 // |
|
|
887 // produce identical results. Also, it would be nice if -1^0.5 |
|
|
888 // produced a pure imaginary result instead of a complex number with a |
|
|
889 // small real part. But perhaps that's really a problem with the math |
|
|
890 // library... |
|
|
891 |
|
|
892 // -*- 1 -*- |
|
|
893 octave_value |
|
|
894 elem_xpow (double a, const NDArray& b) |
|
|
895 { |
|
|
896 octave_value retval; |
|
|
897 |
|
|
898 double d1, d2; |
|
|
899 |
|
|
900 if (a < 0.0 && ! b.all_integers (d1, d2)) |
|
|
901 { |
|
|
902 Complex atmp (a); |
|
|
903 ComplexNDArray result (b.dims ()); |
|
5275
|
904 for (octave_idx_type i = 0; i < b.length (); i++) |
|
4543
|
905 { |
|
|
906 OCTAVE_QUIT; |
|
5260
|
907 result(i) = std::pow (atmp, b(i)); |
|
4543
|
908 } |
|
|
909 |
|
|
910 retval = result; |
|
|
911 } |
|
|
912 else |
|
|
913 { |
|
|
914 NDArray result (b.dims ()); |
|
5275
|
915 for (octave_idx_type i = 0; i < b.length (); i++) |
|
4543
|
916 { |
|
|
917 OCTAVE_QUIT; |
|
5260
|
918 result (i) = std::pow (a, b(i)); |
|
4543
|
919 } |
|
|
920 |
|
|
921 retval = result; |
|
|
922 } |
|
|
923 |
|
|
924 return retval; |
|
|
925 } |
|
|
926 |
|
|
927 // -*- 2 -*- |
|
|
928 octave_value |
|
|
929 elem_xpow (double a, const ComplexNDArray& b) |
|
|
930 { |
|
|
931 ComplexNDArray result (b.dims ()); |
|
|
932 Complex atmp (a); |
|
5275
|
933 |
|
|
934 for (octave_idx_type i = 0; i < b.length (); i++) |
|
4543
|
935 { |
|
|
936 OCTAVE_QUIT; |
|
5260
|
937 result(i) = std::pow (atmp, b(i)); |
|
4543
|
938 } |
|
|
939 |
|
|
940 return result; |
|
|
941 } |
|
|
942 |
|
|
943 // -*- 3 -*- |
|
|
944 octave_value |
|
|
945 elem_xpow (const NDArray& a, double b) |
|
|
946 { |
|
|
947 octave_value retval; |
|
|
948 |
|
|
949 if (static_cast<int> (b) != b && a.any_element_is_negative ()) |
|
|
950 { |
|
|
951 ComplexNDArray result (a.dims ()); |
|
|
952 |
|
5275
|
953 for (octave_idx_type i = 0; i < a.length (); i++) |
|
4543
|
954 { |
|
|
955 OCTAVE_QUIT; |
|
5665
|
956 |
|
|
957 Complex atmp (a (i)); |
|
|
958 |
|
|
959 result(i) = std::pow (atmp, b); |
|
4543
|
960 } |
|
|
961 |
|
|
962 retval = result; |
|
|
963 } |
|
|
964 else |
|
|
965 { |
|
|
966 NDArray result (a.dims ()); |
|
|
967 |
|
5275
|
968 for (octave_idx_type i = 0; i < a.length (); i++) |
|
4543
|
969 { |
|
|
970 OCTAVE_QUIT; |
|
5260
|
971 result(i) = std::pow (a(i), b); |
|
4543
|
972 } |
|
|
973 |
|
|
974 retval = result; |
|
|
975 } |
|
|
976 |
|
|
977 return retval; |
|
|
978 } |
|
|
979 |
|
|
980 // -*- 4 -*- |
|
|
981 octave_value |
|
|
982 elem_xpow (const NDArray& a, const NDArray& b) |
|
|
983 { |
|
|
984 octave_value retval; |
|
|
985 |
|
|
986 dim_vector a_dims = a.dims (); |
|
|
987 dim_vector b_dims = b.dims (); |
|
|
988 |
|
|
989 if (a_dims != b_dims) |
|
|
990 { |
|
|
991 gripe_nonconformant ("operator .^", a_dims, b_dims); |
|
|
992 return octave_value (); |
|
|
993 } |
|
|
994 |
|
|
995 int len = a.length (); |
|
|
996 |
|
|
997 bool convert_to_complex = false; |
|
|
998 |
|
5275
|
999 for (octave_idx_type i = 0; i < len; i++) |
|
4543
|
1000 { |
|
|
1001 OCTAVE_QUIT; |
|
|
1002 double atmp = a(i); |
|
|
1003 double btmp = b(i); |
|
|
1004 if (atmp < 0.0 && static_cast<int> (btmp) != btmp) |
|
|
1005 { |
|
|
1006 convert_to_complex = true; |
|
|
1007 goto done; |
|
|
1008 } |
|
|
1009 } |
|
|
1010 |
|
|
1011 done: |
|
|
1012 |
|
|
1013 if (convert_to_complex) |
|
|
1014 { |
|
|
1015 ComplexNDArray complex_result (a_dims); |
|
|
1016 |
|
5275
|
1017 for (octave_idx_type i = 0; i < len; i++) |
|
4543
|
1018 { |
|
|
1019 OCTAVE_QUIT; |
|
5665
|
1020 Complex atmp (a(i)); |
|
|
1021 Complex btmp (b(i)); |
|
|
1022 complex_result(i) = std::pow (atmp, btmp); |
|
4543
|
1023 } |
|
|
1024 |
|
|
1025 retval = complex_result; |
|
|
1026 } |
|
|
1027 else |
|
|
1028 { |
|
|
1029 NDArray result (a_dims); |
|
|
1030 |
|
5275
|
1031 for (octave_idx_type i = 0; i < len; i++) |
|
4543
|
1032 { |
|
|
1033 OCTAVE_QUIT; |
|
5260
|
1034 result(i) = std::pow (a(i), b(i)); |
|
4543
|
1035 } |
|
|
1036 |
|
|
1037 retval = result; |
|
|
1038 } |
|
|
1039 |
|
|
1040 return retval; |
|
|
1041 } |
|
|
1042 |
|
|
1043 // -*- 5 -*- |
|
|
1044 octave_value |
|
|
1045 elem_xpow (const NDArray& a, const Complex& b) |
|
|
1046 { |
|
|
1047 ComplexNDArray result (a.dims ()); |
|
|
1048 |
|
5275
|
1049 for (octave_idx_type i = 0; i < a.length (); i++) |
|
4543
|
1050 { |
|
|
1051 OCTAVE_QUIT; |
|
5260
|
1052 result(i) = std::pow (Complex (a(i)), b); |
|
4543
|
1053 } |
|
|
1054 |
|
|
1055 return result; |
|
|
1056 } |
|
|
1057 |
|
|
1058 // -*- 6 -*- |
|
|
1059 octave_value |
|
|
1060 elem_xpow (const NDArray& a, const ComplexNDArray& b) |
|
|
1061 { |
|
|
1062 dim_vector a_dims = a.dims (); |
|
|
1063 dim_vector b_dims = b.dims (); |
|
|
1064 |
|
|
1065 if (a_dims != b_dims) |
|
|
1066 { |
|
|
1067 gripe_nonconformant ("operator .^", a_dims, b_dims); |
|
|
1068 return octave_value (); |
|
|
1069 } |
|
|
1070 |
|
|
1071 ComplexNDArray result (a_dims); |
|
5275
|
1072 |
|
|
1073 for (octave_idx_type i = 0; i < a.length (); i++) |
|
4543
|
1074 { |
|
|
1075 OCTAVE_QUIT; |
|
5260
|
1076 result(i) = std::pow (Complex (a(i)), b(i)); |
|
4543
|
1077 } |
|
|
1078 |
|
|
1079 return result; |
|
|
1080 } |
|
|
1081 |
|
|
1082 // -*- 7 -*- |
|
|
1083 octave_value |
|
|
1084 elem_xpow (const Complex& a, const NDArray& b) |
|
|
1085 { |
|
|
1086 ComplexNDArray result (b.dims ()); |
|
5275
|
1087 |
|
|
1088 for (octave_idx_type i = 0; i < b.length (); i++) |
|
4543
|
1089 { |
|
|
1090 OCTAVE_QUIT; |
|
|
1091 double btmp = b(i); |
|
|
1092 if (xisint (btmp)) |
|
5260
|
1093 result(i) = std::pow (a, static_cast<int> (btmp)); |
|
4543
|
1094 else |
|
5260
|
1095 result(i) = std::pow (a, btmp); |
|
4543
|
1096 } |
|
|
1097 |
|
|
1098 return result; |
|
|
1099 } |
|
|
1100 |
|
|
1101 // -*- 8 -*- |
|
|
1102 octave_value |
|
|
1103 elem_xpow (const Complex& a, const ComplexNDArray& b) |
|
|
1104 { |
|
|
1105 ComplexNDArray result (b.dims ()); |
|
5275
|
1106 |
|
|
1107 for (octave_idx_type i = 0; i < b.length (); i++) |
|
4543
|
1108 { |
|
|
1109 OCTAVE_QUIT; |
|
5260
|
1110 result(i) = std::pow (a, b(i)); |
|
4543
|
1111 } |
|
|
1112 |
|
|
1113 return result; |
|
|
1114 } |
|
|
1115 |
|
|
1116 // -*- 9 -*- |
|
|
1117 octave_value |
|
|
1118 elem_xpow (const ComplexNDArray& a, double b) |
|
|
1119 { |
|
|
1120 ComplexNDArray result (a.dims ()); |
|
|
1121 |
|
|
1122 if (xisint (b)) |
|
|
1123 { |
|
5275
|
1124 for (octave_idx_type i = 0; i < a.length (); i++) |
|
4543
|
1125 { |
|
|
1126 OCTAVE_QUIT; |
|
5260
|
1127 result(i) = std::pow (a(i), static_cast<int> (b)); |
|
4543
|
1128 } |
|
|
1129 } |
|
|
1130 else |
|
|
1131 { |
|
5275
|
1132 for (octave_idx_type i = 0; i < a.length (); i++) |
|
4543
|
1133 { |
|
|
1134 OCTAVE_QUIT; |
|
5260
|
1135 result(i) = std::pow (a(i), b); |
|
4543
|
1136 } |
|
|
1137 } |
|
|
1138 |
|
|
1139 return result; |
|
|
1140 } |
|
|
1141 |
|
|
1142 // -*- 10 -*- |
|
|
1143 octave_value |
|
|
1144 elem_xpow (const ComplexNDArray& a, const NDArray& b) |
|
|
1145 { |
|
|
1146 dim_vector a_dims = a.dims (); |
|
|
1147 dim_vector b_dims = b.dims (); |
|
|
1148 |
|
|
1149 if (a_dims != b_dims) |
|
|
1150 { |
|
|
1151 gripe_nonconformant ("operator .^", a_dims, b_dims); |
|
|
1152 return octave_value (); |
|
|
1153 } |
|
|
1154 |
|
|
1155 ComplexNDArray result (a_dims); |
|
5275
|
1156 |
|
|
1157 for (octave_idx_type i = 0; i < a.length (); i++) |
|
4543
|
1158 { |
|
|
1159 OCTAVE_QUIT; |
|
|
1160 double btmp = b(i); |
|
|
1161 if (xisint (btmp)) |
|
5260
|
1162 result(i) = std::pow (a(i), static_cast<int> (btmp)); |
|
4543
|
1163 else |
|
5260
|
1164 result(i) = std::pow (a(i), btmp); |
|
4543
|
1165 } |
|
|
1166 |
|
|
1167 return result; |
|
|
1168 } |
|
|
1169 |
|
|
1170 // -*- 11 -*- |
|
|
1171 octave_value |
|
|
1172 elem_xpow (const ComplexNDArray& a, const Complex& b) |
|
|
1173 { |
|
|
1174 ComplexNDArray result (a.dims ()); |
|
5275
|
1175 |
|
|
1176 for (octave_idx_type i = 0; i < a.length (); i++) |
|
4543
|
1177 { |
|
|
1178 OCTAVE_QUIT; |
|
5260
|
1179 result(i) = std::pow (a(i), b); |
|
4543
|
1180 } |
|
|
1181 |
|
|
1182 return result; |
|
|
1183 } |
|
|
1184 |
|
|
1185 // -*- 12 -*- |
|
|
1186 octave_value |
|
|
1187 elem_xpow (const ComplexNDArray& a, const ComplexNDArray& b) |
|
|
1188 { |
|
|
1189 dim_vector a_dims = a.dims (); |
|
|
1190 dim_vector b_dims = b.dims (); |
|
|
1191 |
|
|
1192 if (a_dims != b_dims) |
|
|
1193 { |
|
|
1194 gripe_nonconformant ("operator .^", a_dims, b_dims); |
|
|
1195 return octave_value (); |
|
|
1196 } |
|
|
1197 |
|
|
1198 ComplexNDArray result (a_dims); |
|
5275
|
1199 |
|
|
1200 for (octave_idx_type i = 0; i < a.length (); i++) |
|
4543
|
1201 { |
|
|
1202 OCTAVE_QUIT; |
|
5260
|
1203 result(i) = std::pow (a(i), b(i)); |
|
4543
|
1204 } |
|
|
1205 |
|
|
1206 return result; |
|
|
1207 } |
|
|
1208 |
|
1
|
1209 /* |
|
|
1210 ;;; Local Variables: *** |
|
|
1211 ;;; mode: C++ *** |
|
|
1212 ;;; End: *** |
|
|
1213 */ |
