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