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