Mercurial > octave-nkf

comparison src/xpow.cc @ 5275:23b37da9fd5b

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 2005-04-08 16:07:35 by jwe]
author jwe
date Fri, 08 Apr 2005 16:07:37 +0000
parents deed800e7bef
children 4c8a2e4e0717
comparison
equal deleted inserted replaced
5274:eae7b40388e9 5275:23b37da9fd5b
40 40
41 #include "error.h" 41 #include "error.h"
42 #include "oct-obj.h" 42 #include "oct-obj.h"
43 #include "utils.h" 43 #include "utils.h"
44 #include "xpow.h" 44 #include "xpow.h"
45
46 #ifdef _OPENMP
47 #include <omp.h>
48 #endif
45 49
46 static inline int 50 static inline int
47 xisint (double x) 51 xisint (double x)
48 { 52 {
49 return (D_NINT (x) == x 53 return (D_NINT (x) == x
82 octave_value 86 octave_value
83 xpow (double a, const Matrix& b) 87 xpow (double a, const Matrix& b)
84 { 88 {
85 octave_value retval; 89 octave_value retval;
86 90
87 int nr = b.rows (); 91 octave_idx_type nr = b.rows ();
88 int nc = b.cols (); 92 octave_idx_type nc = b.cols ();
89 93
90 if (nr == 0 || nc == 0 || nr != nc) 94 if (nr == 0 || nc == 0 || nr != nc)
91 error ("for x^A, A must be square"); 95 error ("for x^A, A must be square");
92 else 96 else
93 { 97 {
94 EIG b_eig (b); 98 EIG b_eig (b);
95 ComplexColumnVector lambda (b_eig.eigenvalues ()); 99 ComplexColumnVector lambda (b_eig.eigenvalues ());
96 ComplexMatrix Q (b_eig.eigenvectors ()); 100 ComplexMatrix Q (b_eig.eigenvectors ());
97 101
98 for (int i = 0; i < nr; i++) 102 for (octave_idx_type i = 0; i < nr; i++)
99 { 103 {
100 Complex elt = lambda (i); 104 Complex elt = lambda (i);
101 if (std::imag (elt) == 0.0) 105 if (std::imag (elt) == 0.0)
102 lambda (i) = std::pow (a, std::real (elt)); 106 lambda (i) = std::pow (a, std::real (elt));
103 else 107 else
125 octave_value 129 octave_value
126 xpow (double a, const ComplexMatrix& b) 130 xpow (double a, const ComplexMatrix& b)
127 { 131 {
128 octave_value retval; 132 octave_value retval;
129 133
130 int nr = b.rows (); 134 octave_idx_type nr = b.rows ();
131 int nc = b.cols (); 135 octave_idx_type nc = b.cols ();
132 136
133 if (nr == 0 || nc == 0 || nr != nc) 137 if (nr == 0 || nc == 0 || nr != nc)
134 error ("for x^A, A must be square"); 138 error ("for x^A, A must be square");
135 else 139 else
136 { 140 {
137 EIG b_eig (b); 141 EIG b_eig (b);
138 ComplexColumnVector lambda (b_eig.eigenvalues ()); 142 ComplexColumnVector lambda (b_eig.eigenvalues ());
139 ComplexMatrix Q (b_eig.eigenvectors ()); 143 ComplexMatrix Q (b_eig.eigenvectors ());
140 144
141 for (int i = 0; i < nr; i++) 145 for (octave_idx_type i = 0; i < nr; i++)
142 { 146 {
143 Complex elt = lambda (i); 147 Complex elt = lambda (i);
144 if (std::imag (elt) == 0.0) 148 if (std::imag (elt) == 0.0)
145 lambda (i) = std::pow (a, std::real (elt)); 149 lambda (i) = std::pow (a, std::real (elt));
146 else 150 else
158 octave_value 162 octave_value
159 xpow (const Matrix& a, double b) 163 xpow (const Matrix& a, double b)
160 { 164 {
161 octave_value retval; 165 octave_value retval;
162 166
163 int nr = a.rows (); 167 octave_idx_type nr = a.rows ();
164 int nc = a.cols (); 168 octave_idx_type nc = a.cols ();
165 169
166 if (nr == 0 || nc == 0 || nr != nc) 170 if (nr == 0 || nc == 0 || nr != nc)
167 error ("for A^b, A must be square"); 171 error ("for A^b, A must be square");
168 else 172 else
169 { 173 {
183 Matrix atmp; 187 Matrix atmp;
184 if (btmp < 0) 188 if (btmp < 0)
185 { 189 {
186 btmp = -btmp; 190 btmp = -btmp;
187 191
188 int info; 192 octave_idx_type info;
189 double rcond = 0.0; 193 double rcond = 0.0;
190 194
191 atmp = a.inverse (info, rcond, 1); 195 atmp = a.inverse (info, rcond, 1);
192 196
193 if (info == -1) 197 if (info == -1)
219 { 223 {
220 EIG a_eig (a); 224 EIG a_eig (a);
221 ComplexColumnVector lambda (a_eig.eigenvalues ()); 225 ComplexColumnVector lambda (a_eig.eigenvalues ());
222 ComplexMatrix Q (a_eig.eigenvectors ()); 226 ComplexMatrix Q (a_eig.eigenvectors ());
223 227
224 for (int i = 0; i < nr; i++) 228 for (octave_idx_type i = 0; i < nr; i++)
225 lambda (i) = std::pow (lambda (i), b); 229 lambda (i) = std::pow (lambda (i), b);
226 230
227 ComplexDiagMatrix D (lambda); 231 ComplexDiagMatrix D (lambda);
228 232
229 retval = ComplexMatrix (Q * D * Q.inverse ()); 233 retval = ComplexMatrix (Q * D * Q.inverse ());
237 octave_value 241 octave_value
238 xpow (const Matrix& a, const Complex& b) 242 xpow (const Matrix& a, const Complex& b)
239 { 243 {
240 octave_value retval; 244 octave_value retval;
241 245
242 int nr = a.rows (); 246 octave_idx_type nr = a.rows ();
243 int nc = a.cols (); 247 octave_idx_type nc = a.cols ();
244 248
245 if (nr == 0 || nc == 0 || nr != nc) 249 if (nr == 0 || nc == 0 || nr != nc)
246 error ("for A^b, A must be square"); 250 error ("for A^b, A must be square");
247 else 251 else
248 { 252 {
249 EIG a_eig (a); 253 EIG a_eig (a);
250 ComplexColumnVector lambda (a_eig.eigenvalues ()); 254 ComplexColumnVector lambda (a_eig.eigenvalues ());
251 ComplexMatrix Q (a_eig.eigenvectors ()); 255 ComplexMatrix Q (a_eig.eigenvectors ());
252 256
253 for (int i = 0; i < nr; i++) 257 for (octave_idx_type i = 0; i < nr; i++)
254 lambda (i) = std::pow (lambda (i), b); 258 lambda (i) = std::pow (lambda (i), b);
255 259
256 ComplexDiagMatrix D (lambda); 260 ComplexDiagMatrix D (lambda);
257 261
258 retval = ComplexMatrix (Q * D * Q.inverse ()); 262 retval = ComplexMatrix (Q * D * Q.inverse ());
279 octave_value 283 octave_value
280 xpow (const Complex& a, const Matrix& b) 284 xpow (const Complex& a, const Matrix& b)
281 { 285 {
282 octave_value retval; 286 octave_value retval;
283 287
284 int nr = b.rows (); 288 octave_idx_type nr = b.rows ();
285 int nc = b.cols (); 289 octave_idx_type nc = b.cols ();
286 290
287 if (nr == 0 || nc == 0 || nr != nc) 291 if (nr == 0 || nc == 0 || nr != nc)
288 error ("for x^A, A must be square"); 292 error ("for x^A, A must be square");
289 else 293 else
290 { 294 {
291 EIG b_eig (b); 295 EIG b_eig (b);
292 ComplexColumnVector lambda (b_eig.eigenvalues ()); 296 ComplexColumnVector lambda (b_eig.eigenvalues ());
293 ComplexMatrix Q (b_eig.eigenvectors ()); 297 ComplexMatrix Q (b_eig.eigenvectors ());
294 298
295 for (int i = 0; i < nr; i++) 299 for (octave_idx_type i = 0; i < nr; i++)
296 { 300 {
297 Complex elt = lambda (i); 301 Complex elt = lambda (i);
298 if (std::imag (elt) == 0.0) 302 if (std::imag (elt) == 0.0)
299 lambda (i) = std::pow (a, std::real (elt)); 303 lambda (i) = std::pow (a, std::real (elt));
300 else 304 else
321 octave_value 325 octave_value
322 xpow (const Complex& a, const ComplexMatrix& b) 326 xpow (const Complex& a, const ComplexMatrix& b)
323 { 327 {
324 octave_value retval; 328 octave_value retval;
325 329
326 int nr = b.rows (); 330 octave_idx_type nr = b.rows ();
327 int nc = b.cols (); 331 octave_idx_type nc = b.cols ();
328 332
329 if (nr == 0 || nc == 0 || nr != nc) 333 if (nr == 0 || nc == 0 || nr != nc)
330 error ("for x^A, A must be square"); 334 error ("for x^A, A must be square");
331 else 335 else
332 { 336 {
333 EIG b_eig (b); 337 EIG b_eig (b);
334 ComplexColumnVector lambda (b_eig.eigenvalues ()); 338 ComplexColumnVector lambda (b_eig.eigenvalues ());
335 ComplexMatrix Q (b_eig.eigenvectors ()); 339 ComplexMatrix Q (b_eig.eigenvectors ());
336 340
337 for (int i = 0; i < nr; i++) 341 for (octave_idx_type i = 0; i < nr; i++)
338 { 342 {
339 Complex elt = lambda (i); 343 Complex elt = lambda (i);
340 if (std::imag (elt) == 0.0) 344 if (std::imag (elt) == 0.0)
341 lambda (i) = std::pow (a, std::real (elt)); 345 lambda (i) = std::pow (a, std::real (elt));
342 else 346 else
354 octave_value 358 octave_value
355 xpow (const ComplexMatrix& a, double b) 359 xpow (const ComplexMatrix& a, double b)
356 { 360 {
357 octave_value retval; 361 octave_value retval;
358 362
359 int nr = a.rows (); 363 octave_idx_type nr = a.rows ();
360 int nc = a.cols (); 364 octave_idx_type nc = a.cols ();
361 365
362 if (nr == 0 || nc == 0 || nr != nc) 366 if (nr == 0 || nc == 0 || nr != nc)
363 error ("for A^b, A must be square"); 367 error ("for A^b, A must be square");
364 else 368 else
365 { 369 {
379 ComplexMatrix atmp; 383 ComplexMatrix atmp;
380 if (btmp < 0) 384 if (btmp < 0)
381 { 385 {
382 btmp = -btmp; 386 btmp = -btmp;
383 387
384 int info; 388 octave_idx_type info;
385 double rcond = 0.0; 389 double rcond = 0.0;
386 390
387 atmp = a.inverse (info, rcond, 1); 391 atmp = a.inverse (info, rcond, 1);
388 392
389 if (info == -1) 393 if (info == -1)
415 { 419 {
416 EIG a_eig (a); 420 EIG a_eig (a);
417 ComplexColumnVector lambda (a_eig.eigenvalues ()); 421 ComplexColumnVector lambda (a_eig.eigenvalues ());
418 ComplexMatrix Q (a_eig.eigenvectors ()); 422 ComplexMatrix Q (a_eig.eigenvectors ());
419 423
420 for (int i = 0; i < nr; i++) 424 for (octave_idx_type i = 0; i < nr; i++)
421 lambda (i) = std::pow (lambda (i), b); 425 lambda (i) = std::pow (lambda (i), b);
422 426
423 ComplexDiagMatrix D (lambda); 427 ComplexDiagMatrix D (lambda);
424 428
425 retval = ComplexMatrix (Q * D * Q.inverse ()); 429 retval = ComplexMatrix (Q * D * Q.inverse ());
433 octave_value 437 octave_value
434 xpow (const ComplexMatrix& a, const Complex& b) 438 xpow (const ComplexMatrix& a, const Complex& b)
435 { 439 {
436 octave_value retval; 440 octave_value retval;
437 441
438 int nr = a.rows (); 442 octave_idx_type nr = a.rows ();
439 int nc = a.cols (); 443 octave_idx_type nc = a.cols ();
440 444
441 if (nr == 0 || nc == 0 || nr != nc) 445 if (nr == 0 || nc == 0 || nr != nc)
442 error ("for A^b, A must be square"); 446 error ("for A^b, A must be square");
443 else 447 else
444 { 448 {
445 EIG a_eig (a); 449 EIG a_eig (a);
446 ComplexColumnVector lambda (a_eig.eigenvalues ()); 450 ComplexColumnVector lambda (a_eig.eigenvalues ());
447 ComplexMatrix Q (a_eig.eigenvectors ()); 451 ComplexMatrix Q (a_eig.eigenvectors ());
448 452
449 for (int i = 0; i < nr; i++) 453 for (octave_idx_type i = 0; i < nr; i++)
450 lambda (i) = std::pow (lambda (i), b); 454 lambda (i) = std::pow (lambda (i), b);
451 455
452 ComplexDiagMatrix D (lambda); 456 ComplexDiagMatrix D (lambda);
453 457
454 retval = ComplexMatrix (Q * D * Q.inverse ()); 458 retval = ComplexMatrix (Q * D * Q.inverse ());
490 octave_value 494 octave_value
491 elem_xpow (double a, const Matrix& b) 495 elem_xpow (double a, const Matrix& b)
492 { 496 {
493 octave_value retval; 497 octave_value retval;
494 498
495 int nr = b.rows (); 499 octave_idx_type nr = b.rows ();
496 int nc = b.cols (); 500 octave_idx_type nc = b.cols ();
497 501
498 double d1, d2; 502 double d1, d2;
499 503
500 if (a < 0.0 && ! b.all_integers (d1, d2)) 504 if (a < 0.0 && ! b.all_integers (d1, d2))
501 { 505 {
502 Complex atmp (a); 506 Complex atmp (a);
503 ComplexMatrix result (nr, nc); 507 ComplexMatrix result (nr, nc);
504 for (int j = 0; j < nc; j++) 508
505 for (int i = 0; i < nr; i++) 509 for (octave_idx_type j = 0; j < nc; j++)
510 for (octave_idx_type i = 0; i < nr; i++)
506 { 511 {
507 OCTAVE_QUIT; 512 OCTAVE_QUIT;
508 result (i, j) = std::pow (atmp, b (i, j)); 513 result (i, j) = std::pow (atmp, b (i, j));
509 } 514 }
510 515
511 retval = result; 516 retval = result;
512 } 517 }
513 else 518 else
514 { 519 {
515 Matrix result (nr, nc); 520 Matrix result (nr, nc);
516 for (int j = 0; j < nc; j++) 521
517 for (int i = 0; i < nr; i++) 522 for (octave_idx_type j = 0; j < nc; j++)
523 for (octave_idx_type i = 0; i < nr; i++)
518 { 524 {
519 OCTAVE_QUIT; 525 OCTAVE_QUIT;
520 result (i, j) = std::pow (a, b (i, j)); 526 result (i, j) = std::pow (a, b (i, j));
521 } 527 }
522 528
528 534
529 // -*- 2 -*- 535 // -*- 2 -*-
530 octave_value 536 octave_value
531 elem_xpow (double a, const ComplexMatrix& b) 537 elem_xpow (double a, const ComplexMatrix& b)
532 { 538 {
533 int nr = b.rows (); 539 octave_idx_type nr = b.rows ();
534 int nc = b.cols (); 540 octave_idx_type nc = b.cols ();
535 541
536 ComplexMatrix result (nr, nc); 542 ComplexMatrix result (nr, nc);
537 Complex atmp (a); 543 Complex atmp (a);
538 for (int j = 0; j < nc; j++) 544
539 for (int i = 0; i < nr; i++) 545 for (octave_idx_type j = 0; j < nc; j++)
546 for (octave_idx_type i = 0; i < nr; i++)
540 { 547 {
541 OCTAVE_QUIT; 548 OCTAVE_QUIT;
542 result (i, j) = std::pow (atmp, b (i, j)); 549 result (i, j) = std::pow (atmp, b (i, j));
543 } 550 }
544 551
549 octave_value 556 octave_value
550 elem_xpow (const Matrix& a, double b) 557 elem_xpow (const Matrix& a, double b)
551 { 558 {
552 octave_value retval; 559 octave_value retval;
553 560
554 int nr = a.rows (); 561 octave_idx_type nr = a.rows ();
555 int nc = a.cols (); 562 octave_idx_type nc = a.cols ();
556 563
557 if (static_cast<int> (b) != b && a.any_element_is_negative ()) 564 if (static_cast<int> (b) != b && a.any_element_is_negative ())
558 { 565 {
559 ComplexMatrix result (nr, nc); 566 ComplexMatrix result (nr, nc);
560 for (int j = 0; j < nc; j++) 567
561 for (int i = 0; i < nr; i++) 568 for (octave_idx_type j = 0; j < nc; j++)
569 for (octave_idx_type i = 0; i < nr; i++)
562 { 570 {
563 OCTAVE_QUIT; 571 OCTAVE_QUIT;
564 572
565 Complex atmp (a (i, j)); 573 Complex atmp (a (i, j));
566 574
567 result (i, j) = std::pow (atmp, b); 575 result (i, j) = std::pow (atmp, b);
568 } 576 }
570 retval = result; 578 retval = result;
571 } 579 }
572 else 580 else
573 { 581 {
574 Matrix result (nr, nc); 582 Matrix result (nr, nc);
575 for (int j = 0; j < nc; j++) 583
576 for (int i = 0; i < nr; i++) 584 for (octave_idx_type j = 0; j < nc; j++)
585 for (octave_idx_type i = 0; i < nr; i++)
577 { 586 {
578 OCTAVE_QUIT; 587 OCTAVE_QUIT;
579 result (i, j) = std::pow (a (i, j), b); 588 result (i, j) = std::pow (a (i, j), b);
580 } 589 }
581 590
589 octave_value 598 octave_value
590 elem_xpow (const Matrix& a, const Matrix& b) 599 elem_xpow (const Matrix& a, const Matrix& b)
591 { 600 {
592 octave_value retval; 601 octave_value retval;
593 602
594 int nr = a.rows (); 603 octave_idx_type nr = a.rows ();
595 int nc = a.cols (); 604 octave_idx_type nc = a.cols ();
596 605
597 int b_nr = b.rows (); 606 octave_idx_type b_nr = b.rows ();
598 int b_nc = b.cols (); 607 octave_idx_type b_nc = b.cols ();
599 608
600 if (nr != b_nr || nc != b_nc) 609 if (nr != b_nr || nc != b_nc)
601 { 610 {
602 gripe_nonconformant ("operator .^", nr, nc, b_nr, b_nc); 611 gripe_nonconformant ("operator .^", nr, nc, b_nr, b_nc);
603 return octave_value (); 612 return octave_value ();
604 } 613 }
605 614
606 int convert_to_complex = 0; 615 int convert_to_complex = 0;
607 for (int j = 0; j < nc; j++) 616 for (octave_idx_type j = 0; j < nc; j++)
608 for (int i = 0; i < nr; i++) 617 for (octave_idx_type i = 0; i < nr; i++)
609 { 618 {
610 OCTAVE_QUIT; 619 OCTAVE_QUIT;
611 double atmp = a (i, j); 620 double atmp = a (i, j);
612 double btmp = b (i, j); 621 double btmp = b (i, j);
613 if (atmp < 0.0 && static_cast<int> (btmp) != btmp) 622 if (atmp < 0.0 && static_cast<int> (btmp) != btmp)
621 630
622 if (convert_to_complex) 631 if (convert_to_complex)
623 { 632 {
624 ComplexMatrix complex_result (nr, nc); 633 ComplexMatrix complex_result (nr, nc);
625 634
626 for (int j = 0; j < nc; j++) 635 for (octave_idx_type j = 0; j < nc; j++)
627 for (int i = 0; i < nr; i++) 636 for (octave_idx_type i = 0; i < nr; i++)
628 { 637 {
629 OCTAVE_QUIT; 638 OCTAVE_QUIT;
630 Complex atmp (a (i, j)); 639 Complex atmp (a (i, j));
631 Complex btmp (b (i, j)); 640 Complex btmp (b (i, j));
632 complex_result (i, j) = std::pow (atmp, btmp); 641 complex_result (i, j) = std::pow (atmp, btmp);
636 } 645 }
637 else 646 else
638 { 647 {
639 Matrix result (nr, nc); 648 Matrix result (nr, nc);
640 649
641 for (int j = 0; j < nc; j++) 650 for (octave_idx_type j = 0; j < nc; j++)
642 for (int i = 0; i < nr; i++) 651 for (octave_idx_type i = 0; i < nr; i++)
643 { 652 {
644 OCTAVE_QUIT; 653 OCTAVE_QUIT;
645 result (i, j) = std::pow (a (i, j), b (i, j)); 654 result (i, j) = std::pow (a (i, j), b (i, j));
646 } 655 }
647 656
653 662
654 // -*- 5 -*- 663 // -*- 5 -*-
655 octave_value 664 octave_value
656 elem_xpow (const Matrix& a, const Complex& b) 665 elem_xpow (const Matrix& a, const Complex& b)
657 { 666 {
658 int nr = a.rows (); 667 octave_idx_type nr = a.rows ();
659 int nc = a.cols (); 668 octave_idx_type nc = a.cols ();
660 669
661 ComplexMatrix result (nr, nc); 670 ComplexMatrix result (nr, nc);
662 for (int j = 0; j < nc; j++) 671
663 for (int i = 0; i < nr; i++) 672 for (octave_idx_type j = 0; j < nc; j++)
673 for (octave_idx_type i = 0; i < nr; i++)
664 { 674 {
665 OCTAVE_QUIT; 675 OCTAVE_QUIT;
666 result (i, j) = std::pow (Complex (a (i, j)), b); 676 result (i, j) = std::pow (Complex (a (i, j)), b);
667 } 677 }
668 678
671 681
672 // -*- 6 -*- 682 // -*- 6 -*-
673 octave_value 683 octave_value
674 elem_xpow (const Matrix& a, const ComplexMatrix& b) 684 elem_xpow (const Matrix& a, const ComplexMatrix& b)
675 { 685 {
676 int nr = a.rows (); 686 octave_idx_type nr = a.rows ();
677 int nc = a.cols (); 687 octave_idx_type nc = a.cols ();
678 688
679 int b_nr = b.rows (); 689 octave_idx_type b_nr = b.rows ();
680 int b_nc = b.cols (); 690 octave_idx_type b_nc = b.cols ();
681 691
682 if (nr != b_nr || nc != b_nc) 692 if (nr != b_nr || nc != b_nc)
683 { 693 {
684 gripe_nonconformant ("operator .^", nr, nc, b_nr, b_nc); 694 gripe_nonconformant ("operator .^", nr, nc, b_nr, b_nc);
685 return octave_value (); 695 return octave_value ();
686 } 696 }
687 697
688 ComplexMatrix result (nr, nc); 698 ComplexMatrix result (nr, nc);
689 for (int j = 0; j < nc; j++) 699
690 for (int i = 0; i < nr; i++) 700 for (octave_idx_type j = 0; j < nc; j++)
701 for (octave_idx_type i = 0; i < nr; i++)
691 { 702 {
692 OCTAVE_QUIT; 703 OCTAVE_QUIT;
693 result (i, j) = std::pow (Complex (a (i, j)), b (i, j)); 704 result (i, j) = std::pow (Complex (a (i, j)), b (i, j));
694 } 705 }
695 706
698 709
699 // -*- 7 -*- 710 // -*- 7 -*-
700 octave_value 711 octave_value
701 elem_xpow (const Complex& a, const Matrix& b) 712 elem_xpow (const Complex& a, const Matrix& b)
702 { 713 {
703 int nr = b.rows (); 714 octave_idx_type nr = b.rows ();
704 int nc = b.cols (); 715 octave_idx_type nc = b.cols ();
705 716
706 ComplexMatrix result (nr, nc); 717 ComplexMatrix result (nr, nc);
707 for (int j = 0; j < nc; j++) 718
708 for (int i = 0; i < nr; i++) 719 for (octave_idx_type j = 0; j < nc; j++)
720 for (octave_idx_type i = 0; i < nr; i++)
709 { 721 {
710 OCTAVE_QUIT; 722 OCTAVE_QUIT;
711 double btmp = b (i, j); 723 double btmp = b (i, j);
712 if (xisint (btmp)) 724 if (xisint (btmp))
713 result (i, j) = std::pow (a, static_cast<int> (btmp)); 725 result (i, j) = std::pow (a, static_cast<int> (btmp));
720 732
721 // -*- 8 -*- 733 // -*- 8 -*-
722 octave_value 734 octave_value
723 elem_xpow (const Complex& a, const ComplexMatrix& b) 735 elem_xpow (const Complex& a, const ComplexMatrix& b)
724 { 736 {
725 int nr = b.rows (); 737 octave_idx_type nr = b.rows ();
726 int nc = b.cols (); 738 octave_idx_type nc = b.cols ();
727 739
728 ComplexMatrix result (nr, nc); 740 ComplexMatrix result (nr, nc);
729 for (int j = 0; j < nc; j++) 741
730 for (int i = 0; i < nr; i++) 742 for (octave_idx_type j = 0; j < nc; j++)
743 for (octave_idx_type i = 0; i < nr; i++)
731 { 744 {
732 OCTAVE_QUIT; 745 OCTAVE_QUIT;
733 result (i, j) = std::pow (a, b (i, j)); 746 result (i, j) = std::pow (a, b (i, j));
734 } 747 }
735 748
738 751
739 // -*- 9 -*- 752 // -*- 9 -*-
740 octave_value 753 octave_value
741 elem_xpow (const ComplexMatrix& a, double b) 754 elem_xpow (const ComplexMatrix& a, double b)
742 { 755 {
743 int nr = a.rows (); 756 octave_idx_type nr = a.rows ();
744 int nc = a.cols (); 757 octave_idx_type nc = a.cols ();
745 758
746 ComplexMatrix result (nr, nc); 759 ComplexMatrix result (nr, nc);
747 760
748 if (xisint (b)) 761 if (xisint (b))
749 { 762 {
750 for (int j = 0; j < nc; j++) 763 for (octave_idx_type j = 0; j < nc; j++)
751 for (int i = 0; i < nr; i++) 764 for (octave_idx_type i = 0; i < nr; i++)
752 { 765 {
753 OCTAVE_QUIT; 766 OCTAVE_QUIT;
754 result (i, j) = std::pow (a (i, j), static_cast<int> (b)); 767 result (i, j) = std::pow (a (i, j), static_cast<int> (b));
755 } 768 }
756 } 769 }
757 else 770 else
758 { 771 {
759 for (int j = 0; j < nc; j++) 772 for (octave_idx_type j = 0; j < nc; j++)
760 for (int i = 0; i < nr; i++) 773 for (octave_idx_type i = 0; i < nr; i++)
761 { 774 {
762 OCTAVE_QUIT; 775 OCTAVE_QUIT;
763 result (i, j) = std::pow (a (i, j), b); 776 result (i, j) = std::pow (a (i, j), b);
764 } 777 }
765 } 778 }
769 782
770 // -*- 10 -*- 783 // -*- 10 -*-
771 octave_value 784 octave_value
772 elem_xpow (const ComplexMatrix& a, const Matrix& b) 785 elem_xpow (const ComplexMatrix& a, const Matrix& b)
773 { 786 {
774 int nr = a.rows (); 787 octave_idx_type nr = a.rows ();
775 int nc = a.cols (); 788 octave_idx_type nc = a.cols ();
776 789
777 int b_nr = b.rows (); 790 octave_idx_type b_nr = b.rows ();
778 int b_nc = b.cols (); 791 octave_idx_type b_nc = b.cols ();
779 792
780 if (nr != b_nr || nc != b_nc) 793 if (nr != b_nr || nc != b_nc)
781 { 794 {
782 gripe_nonconformant ("operator .^", nr, nc, b_nr, b_nc); 795 gripe_nonconformant ("operator .^", nr, nc, b_nr, b_nc);
783 return octave_value (); 796 return octave_value ();
784 } 797 }
785 798
786 ComplexMatrix result (nr, nc); 799 ComplexMatrix result (nr, nc);
787 for (int j = 0; j < nc; j++) 800
788 for (int i = 0; i < nr; i++) 801 for (octave_idx_type j = 0; j < nc; j++)
802 for (octave_idx_type i = 0; i < nr; i++)
789 { 803 {
790 OCTAVE_QUIT; 804 OCTAVE_QUIT;
791 double btmp = b (i, j); 805 double btmp = b (i, j);
792 if (xisint (btmp)) 806 if (xisint (btmp))
793 result (i, j) = std::pow (a (i, j), static_cast<int> (btmp)); 807 result (i, j) = std::pow (a (i, j), static_cast<int> (btmp));
800 814
801 // -*- 11 -*- 815 // -*- 11 -*-
802 octave_value 816 octave_value
803 elem_xpow (const ComplexMatrix& a, const Complex& b) 817 elem_xpow (const ComplexMatrix& a, const Complex& b)
804 { 818 {
805 int nr = a.rows (); 819 octave_idx_type nr = a.rows ();
806 int nc = a.cols (); 820 octave_idx_type nc = a.cols ();
807 821
808 ComplexMatrix result (nr, nc); 822 ComplexMatrix result (nr, nc);
809 for (int j = 0; j < nc; j++) 823
810 for (int i = 0; i < nr; i++) 824 for (octave_idx_type j = 0; j < nc; j++)
825 for (octave_idx_type i = 0; i < nr; i++)
811 { 826 {
812 OCTAVE_QUIT; 827 OCTAVE_QUIT;
813 result (i, j) = std::pow (a (i, j), b); 828 result (i, j) = std::pow (a (i, j), b);
814 } 829 }
815 830
818 833
819 // -*- 12 -*- 834 // -*- 12 -*-
820 octave_value 835 octave_value
821 elem_xpow (const ComplexMatrix& a, const ComplexMatrix& b) 836 elem_xpow (const ComplexMatrix& a, const ComplexMatrix& b)
822 { 837 {
823 int nr = a.rows (); 838 octave_idx_type nr = a.rows ();
824 int nc = a.cols (); 839 octave_idx_type nc = a.cols ();
825 840
826 int b_nr = b.rows (); 841 octave_idx_type b_nr = b.rows ();
827 int b_nc = b.cols (); 842 octave_idx_type b_nc = b.cols ();
828 843
829 if (nr != b_nr || nc != b_nc) 844 if (nr != b_nr || nc != b_nc)
830 { 845 {
831 gripe_nonconformant ("operator .^", nr, nc, b_nr, b_nc); 846 gripe_nonconformant ("operator .^", nr, nc, b_nr, b_nc);
832 return octave_value (); 847 return octave_value ();
833 } 848 }
834 849
835 ComplexMatrix result (nr, nc); 850 ComplexMatrix result (nr, nc);
836 for (int j = 0; j < nc; j++) 851
837 for (int i = 0; i < nr; i++) 852 for (octave_idx_type j = 0; j < nc; j++)
853 for (octave_idx_type i = 0; i < nr; i++)
838 { 854 {
839 OCTAVE_QUIT; 855 OCTAVE_QUIT;
840 result (i, j) = std::pow (a (i, j), b (i, j)); 856 result (i, j) = std::pow (a (i, j), b (i, j));
841 } 857 }
842 858
882 898
883 if (a < 0.0 && ! b.all_integers (d1, d2)) 899 if (a < 0.0 && ! b.all_integers (d1, d2))
884 { 900 {
885 Complex atmp (a); 901 Complex atmp (a);
886 ComplexNDArray result (b.dims ()); 902 ComplexNDArray result (b.dims ());
887 for (int i = 0; i < b.length (); i++) 903 for (octave_idx_type i = 0; i < b.length (); i++)
888 { 904 {
889 OCTAVE_QUIT; 905 OCTAVE_QUIT;
890 result(i) = std::pow (atmp, b(i)); 906 result(i) = std::pow (atmp, b(i));
891 } 907 }
892 908
893 retval = result; 909 retval = result;
894 } 910 }
895 else 911 else
896 { 912 {
897 NDArray result (b.dims ()); 913 NDArray result (b.dims ());
898 for (int i = 0; i < b.length (); i++) 914 for (octave_idx_type i = 0; i < b.length (); i++)
899 { 915 {
900 OCTAVE_QUIT; 916 OCTAVE_QUIT;
901 result (i) = std::pow (a, b(i)); 917 result (i) = std::pow (a, b(i));
902 } 918 }
903 919
911 octave_value 927 octave_value
912 elem_xpow (double a, const ComplexNDArray& b) 928 elem_xpow (double a, const ComplexNDArray& b)
913 { 929 {
914 ComplexNDArray result (b.dims ()); 930 ComplexNDArray result (b.dims ());
915 Complex atmp (a); 931 Complex atmp (a);
916 for (int i = 0; i < b.length (); i++) 932
933 for (octave_idx_type i = 0; i < b.length (); i++)
917 { 934 {
918 OCTAVE_QUIT; 935 OCTAVE_QUIT;
919 result(i) = std::pow (atmp, b(i)); 936 result(i) = std::pow (atmp, b(i));
920 } 937 }
921 938
930 947
931 if (static_cast<int> (b) != b && a.any_element_is_negative ()) 948 if (static_cast<int> (b) != b && a.any_element_is_negative ())
932 { 949 {
933 ComplexNDArray result (a.dims ()); 950 ComplexNDArray result (a.dims ());
934 951
935 for (int i = 0; i < a.length (); i++) 952 for (octave_idx_type i = 0; i < a.length (); i++)
936 { 953 {
937 OCTAVE_QUIT; 954 OCTAVE_QUIT;
938 955
939 Complex atmp (a (i)); 956 Complex atmp (a (i));
940 957
945 } 962 }
946 else 963 else
947 { 964 {
948 NDArray result (a.dims ()); 965 NDArray result (a.dims ());
949 966
950 for (int i = 0; i < a.length (); i++) 967 for (octave_idx_type i = 0; i < a.length (); i++)
951 { 968 {
952 OCTAVE_QUIT; 969 OCTAVE_QUIT;
953 result(i) = std::pow (a(i), b); 970 result(i) = std::pow (a(i), b);
954 } 971 }
955 972
976 993
977 int len = a.length (); 994 int len = a.length ();
978 995
979 bool convert_to_complex = false; 996 bool convert_to_complex = false;
980 997
981 for (int i = 0; i < len; i++) 998 for (octave_idx_type i = 0; i < len; i++)
982 { 999 {
983 OCTAVE_QUIT; 1000 OCTAVE_QUIT;
984 double atmp = a(i); 1001 double atmp = a(i);
985 double btmp = b(i); 1002 double btmp = b(i);
986 if (atmp < 0.0 && static_cast<int> (btmp) != btmp) 1003 if (atmp < 0.0 && static_cast<int> (btmp) != btmp)
994 1011
995 if (convert_to_complex) 1012 if (convert_to_complex)
996 { 1013 {
997 ComplexNDArray complex_result (a_dims); 1014 ComplexNDArray complex_result (a_dims);
998 1015
999 for (int i = 0; i < len; i++) 1016 for (octave_idx_type i = 0; i < len; i++)
1000 { 1017 {
1001 OCTAVE_QUIT; 1018 OCTAVE_QUIT;
1002 Complex atmp (a(i)); 1019 Complex atmp (a(i));
1003 Complex btmp (b(i)); 1020 Complex btmp (b(i));
1004 complex_result(i) = std::pow (atmp, btmp); 1021 complex_result(i) = std::pow (atmp, btmp);
1008 } 1025 }
1009 else 1026 else
1010 { 1027 {
1011 NDArray result (a_dims); 1028 NDArray result (a_dims);
1012 1029
1013 for (int i = 0; i < len; i++) 1030 for (octave_idx_type i = 0; i < len; i++)
1014 { 1031 {
1015 OCTAVE_QUIT; 1032 OCTAVE_QUIT;
1016 result(i) = std::pow (a(i), b(i)); 1033 result(i) = std::pow (a(i), b(i));
1017 } 1034 }
1018 1035
1026 octave_value 1043 octave_value
1027 elem_xpow (const NDArray& a, const Complex& b) 1044 elem_xpow (const NDArray& a, const Complex& b)
1028 { 1045 {
1029 ComplexNDArray result (a.dims ()); 1046 ComplexNDArray result (a.dims ());
1030 1047
1031 for (int i = 0; i < a.length (); i++) 1048 for (octave_idx_type i = 0; i < a.length (); i++)
1032 { 1049 {
1033 OCTAVE_QUIT; 1050 OCTAVE_QUIT;
1034 result(i) = std::pow (Complex (a(i)), b); 1051 result(i) = std::pow (Complex (a(i)), b);
1035 } 1052 }
1036 1053
1049 gripe_nonconformant ("operator .^", a_dims, b_dims); 1066 gripe_nonconformant ("operator .^", a_dims, b_dims);
1050 return octave_value (); 1067 return octave_value ();
1051 } 1068 }
1052 1069
1053 ComplexNDArray result (a_dims); 1070 ComplexNDArray result (a_dims);
1054 for (int i = 0; i < a.length (); i++) 1071
1072 for (octave_idx_type i = 0; i < a.length (); i++)
1055 { 1073 {
1056 OCTAVE_QUIT; 1074 OCTAVE_QUIT;
1057 result(i) = std::pow (Complex (a(i)), b(i)); 1075 result(i) = std::pow (Complex (a(i)), b(i));
1058 } 1076 }
1059 1077
1063 // -*- 7 -*- 1081 // -*- 7 -*-
1064 octave_value 1082 octave_value
1065 elem_xpow (const Complex& a, const NDArray& b) 1083 elem_xpow (const Complex& a, const NDArray& b)
1066 { 1084 {
1067 ComplexNDArray result (b.dims ()); 1085 ComplexNDArray result (b.dims ());
1068 for (int i = 0; i < b.length (); i++) 1086
1087 for (octave_idx_type i = 0; i < b.length (); i++)
1069 { 1088 {
1070 OCTAVE_QUIT; 1089 OCTAVE_QUIT;
1071 double btmp = b(i); 1090 double btmp = b(i);
1072 if (xisint (btmp)) 1091 if (xisint (btmp))
1073 result(i) = std::pow (a, static_cast<int> (btmp)); 1092 result(i) = std::pow (a, static_cast<int> (btmp));
1081 // -*- 8 -*- 1100 // -*- 8 -*-
1082 octave_value 1101 octave_value
1083 elem_xpow (const Complex& a, const ComplexNDArray& b) 1102 elem_xpow (const Complex& a, const ComplexNDArray& b)
1084 { 1103 {
1085 ComplexNDArray result (b.dims ()); 1104 ComplexNDArray result (b.dims ());
1086 for (int i = 0; i < b.length (); i++) 1105
1106 for (octave_idx_type i = 0; i < b.length (); i++)
1087 { 1107 {
1088 OCTAVE_QUIT; 1108 OCTAVE_QUIT;
1089 result(i) = std::pow (a, b(i)); 1109 result(i) = std::pow (a, b(i));
1090 } 1110 }
1091 1111
1098 { 1118 {
1099 ComplexNDArray result (a.dims ()); 1119 ComplexNDArray result (a.dims ());
1100 1120
1101 if (xisint (b)) 1121 if (xisint (b))
1102 { 1122 {
1103 for (int i = 0; i < a.length (); i++) 1123 for (octave_idx_type i = 0; i < a.length (); i++)
1104 { 1124 {
1105 OCTAVE_QUIT; 1125 OCTAVE_QUIT;
1106 result(i) = std::pow (a(i), static_cast<int> (b)); 1126 result(i) = std::pow (a(i), static_cast<int> (b));
1107 } 1127 }
1108 } 1128 }
1109 else 1129 else
1110 { 1130 {
1111 for (int i = 0; i < a.length (); i++) 1131 for (octave_idx_type i = 0; i < a.length (); i++)
1112 { 1132 {
1113 OCTAVE_QUIT; 1133 OCTAVE_QUIT;
1114 result(i) = std::pow (a(i), b); 1134 result(i) = std::pow (a(i), b);
1115 } 1135 }
1116 } 1136 }
1130 gripe_nonconformant ("operator .^", a_dims, b_dims); 1150 gripe_nonconformant ("operator .^", a_dims, b_dims);
1131 return octave_value (); 1151 return octave_value ();
1132 } 1152 }
1133 1153
1134 ComplexNDArray result (a_dims); 1154 ComplexNDArray result (a_dims);
1135 for (int i = 0; i < a.length (); i++) 1155
1156 for (octave_idx_type i = 0; i < a.length (); i++)
1136 { 1157 {
1137 OCTAVE_QUIT; 1158 OCTAVE_QUIT;
1138 double btmp = b(i); 1159 double btmp = b(i);
1139 if (xisint (btmp)) 1160 if (xisint (btmp))
1140 result(i) = std::pow (a(i), static_cast<int> (btmp)); 1161 result(i) = std::pow (a(i), static_cast<int> (btmp));
1148 // -*- 11 -*- 1169 // -*- 11 -*-
1149 octave_value 1170 octave_value
1150 elem_xpow (const ComplexNDArray& a, const Complex& b) 1171 elem_xpow (const ComplexNDArray& a, const Complex& b)
1151 { 1172 {
1152 ComplexNDArray result (a.dims ()); 1173 ComplexNDArray result (a.dims ());
1153 for (int i = 0; i < a.length (); i++) 1174
1175 for (octave_idx_type i = 0; i < a.length (); i++)
1154 { 1176 {
1155 OCTAVE_QUIT; 1177 OCTAVE_QUIT;
1156 result(i) = std::pow (a(i), b); 1178 result(i) = std::pow (a(i), b);
1157 } 1179 }
1158 1180
1171 gripe_nonconformant ("operator .^", a_dims, b_dims); 1193 gripe_nonconformant ("operator .^", a_dims, b_dims);
1172 return octave_value (); 1194 return octave_value ();
1173 } 1195 }
1174 1196
1175 ComplexNDArray result (a_dims); 1197 ComplexNDArray result (a_dims);
1176 for (int i = 0; i < a.length (); i++) 1198
1199 for (octave_idx_type i = 0; i < a.length (); i++)
1177 { 1200 {
1178 OCTAVE_QUIT; 1201 OCTAVE_QUIT;
1179 result(i) = std::pow (a(i), b(i)); 1202 result(i) = std::pow (a(i), b(i));
1180 } 1203 }
1181 1204