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