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