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