Mercurial > hg > octave-nkf

comparison src/xpow.cc @ 2365:7c60f8a6e6a0

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 1996-10-11 23:03:52 by jwe]
author jwe
date Fri, 11 Oct 1996 23:05:38 +0000
parents 5a3f1d00a474
children 9aeba8e006a4
comparison
equal deleted inserted replaced
2364:5eb0af0730d6 2365:7c60f8a6e6a0
34 #include "dDiagMatrix.h" 34 #include "dDiagMatrix.h"
35 #include "dMatrix.h" 35 #include "dMatrix.h"
36 #include "oct-cmplx.h" 36 #include "oct-cmplx.h"
37 37
38 #include "error.h" 38 #include "error.h"
39 #include "pt-const.h" 39 #include "ov.h"
40 #include "utils.h" 40 #include "utils.h"
41 #include "xpow.h" 41 #include "xpow.h"
42
43 // This function also appears in tree-const.cc. Maybe it should be a
44 // member function of the Matrix class.
45
46 static int
47 any_element_is_negative (const Matrix& a)
48 {
49 int nr = a.rows ();
50 int nc = a.columns ();
51 for (int j = 0; j < nc; j++)
52 for (int i = 0; i < nr; i++)
53 if (a (i, j) < 0.0)
54 return 1;
55 return 0;
56 }
57 42
58 static inline int 43 static inline int
59 xisint (double x) 44 xisint (double x)
60 { 45 {
61 return (D_NINT (x) == x 46 return (D_NINT (x) == x
67 // 52 //
68 // op2 \ op1: s m cs cm 53 // op2 \ op1: s m cs cm
69 // +-- +---+---+----+----+ 54 // +-- +---+---+----+----+
70 // scalar | | 1 | 5 | 7 | 11 | 55 // scalar | | 1 | 5 | 7 | 11 |
71 // +---+---+----+----+ 56 // +---+---+----+----+
72 // matrix | 2 | E | 8 | E | 57 // matrix | 2 | * | 8 | * |
73 // +---+---+----+----+ 58 // +---+---+----+----+
74 // complex_scalar | 3 | 6 | 9 | 12 | 59 // complex_scalar | 3 | 6 | 9 | 12 |
75 // +---+---+----+----+ 60 // +---+---+----+----+
76 // complex_matrix | 4 | E | 10 | E | 61 // complex_matrix | 4 | * | 10 | * |
77 // +---+---+----+----+ 62 // +---+---+----+----+
78 //
79 // E -> error, trapped in arith-ops.cc.
80 63
81 // -*- 1 -*- 64 // -*- 1 -*-
82 octave_value 65 octave_value
83 xpow (double a, double b) 66 xpow (double a, double b)
84 { 67 {
96 xpow (double a, const Matrix& b) 79 xpow (double a, const Matrix& b)
97 { 80 {
98 octave_value retval; 81 octave_value retval;
99 82
100 int nr = b.rows (); 83 int nr = b.rows ();
101 int nc = b.columns (); 84 int nc = b.cols ();
102 85
103 if (nr == 0 || nc == 0 || nr != nc) 86 if (nr == 0 || nc == 0 || nr != nc)
104 error ("for x^A, A must be square"); 87 error ("for x^A, A must be square");
105 else 88 else
106 { 89 {
139 xpow (double a, const ComplexMatrix& b) 122 xpow (double a, const ComplexMatrix& b)
140 { 123 {
141 octave_value retval; 124 octave_value retval;
142 125
143 int nr = b.rows (); 126 int nr = b.rows ();
144 int nc = b.columns (); 127 int nc = b.cols ();
145 128
146 if (nr == 0 || nc == 0 || nr != nc) 129 if (nr == 0 || nc == 0 || nr != nc)
147 { 130 error ("for x^A, A must be square");
148 error ("for x^A, A must be square");
149 }
150 else 131 else
151 { 132 {
152 EIG b_eig (b); 133 EIG b_eig (b);
153 ComplexColumnVector lambda (b_eig.eigenvalues ()); 134 ComplexColumnVector lambda (b_eig.eigenvalues ());
154 ComplexMatrix Q (b_eig.eigenvectors ()); 135 ComplexMatrix Q (b_eig.eigenvectors ());
174 xpow (const Matrix& a, double b) 155 xpow (const Matrix& a, double b)
175 { 156 {
176 octave_value retval; 157 octave_value retval;
177 158
178 int nr = a.rows (); 159 int nr = a.rows ();
179 int nc = a.columns (); 160 int nc = a.cols ();
180 161
181 if (nr == 0 || nc == 0 || nr != nc) 162 if (nr == 0 || nc == 0 || nr != nc)
182 { 163 error ("for A^b, A must be square");
183 error ("for A^b, A must be square");
184 }
185 else 164 else
186 { 165 {
187 if ((int) b == b) 166 if ((int) b == b)
188 { 167 {
189 int btmp = (int) b; 168 int btmp = (int) b;
244 xpow (const Matrix& a, const Complex& b) 223 xpow (const Matrix& a, const Complex& b)
245 { 224 {
246 octave_value retval; 225 octave_value retval;
247 226
248 int nr = a.rows (); 227 int nr = a.rows ();
249 int nc = a.columns (); 228 int nc = a.cols ();
250 229
251 if (nr == 0 || nc == 0 || nr != nc) 230 if (nr == 0 || nc == 0 || nr != nc)
252 { 231 error ("for A^b, A must be square");
253 error ("for A^b, A must be square");
254 }
255 else 232 else
256 { 233 {
257 EIG a_eig (a); 234 EIG a_eig (a);
258 ComplexColumnVector lambda (a_eig.eigenvalues ()); 235 ComplexColumnVector lambda (a_eig.eigenvalues ());
259 ComplexMatrix Q (a_eig.eigenvectors ()); 236 ComplexMatrix Q (a_eig.eigenvectors ());
288 xpow (const Complex& a, const Matrix& b) 265 xpow (const Complex& a, const Matrix& b)
289 { 266 {
290 octave_value retval; 267 octave_value retval;
291 268
292 int nr = b.rows (); 269 int nr = b.rows ();
293 int nc = b.columns (); 270 int nc = b.cols ();
294 271
295 if (nr == 0 || nc == 0 || nr != nc) 272 if (nr == 0 || nc == 0 || nr != nc)
296 { 273 error ("for x^A, A must be square");
297 error ("for x^A, A must be square");
298 }
299 else 274 else
300 { 275 {
301 EIG b_eig (b); 276 EIG b_eig (b);
302 ComplexColumnVector lambda (b_eig.eigenvalues ()); 277 ComplexColumnVector lambda (b_eig.eigenvalues ());
303 ComplexMatrix Q (b_eig.eigenvectors ()); 278 ComplexMatrix Q (b_eig.eigenvectors ());
332 xpow (const Complex& a, const ComplexMatrix& b) 307 xpow (const Complex& a, const ComplexMatrix& b)
333 { 308 {
334 octave_value retval; 309 octave_value retval;
335 310
336 int nr = b.rows (); 311 int nr = b.rows ();
337 int nc = b.columns (); 312 int nc = b.cols ();
338 313
339 if (nr == 0 || nc == 0 || nr != nc) 314 if (nr == 0 || nc == 0 || nr != nc)
340 { 315 error ("for x^A, A must be square");
341 error ("for x^A, A must be square");
342 }
343 else 316 else
344 { 317 {
345 EIG b_eig (b); 318 EIG b_eig (b);
346 ComplexColumnVector lambda (b_eig.eigenvalues ()); 319 ComplexColumnVector lambda (b_eig.eigenvalues ());
347 ComplexMatrix Q (b_eig.eigenvectors ()); 320 ComplexMatrix Q (b_eig.eigenvectors ());
367 xpow (const ComplexMatrix& a, double b) 340 xpow (const ComplexMatrix& a, double b)
368 { 341 {
369 octave_value retval; 342 octave_value retval;
370 343
371 int nr = a.rows (); 344 int nr = a.rows ();
372 int nc = a.columns (); 345 int nc = a.cols ();
373 346
374 if (nr == 0 || nc == 0 || nr != nc) 347 if (nr == 0 || nc == 0 || nr != nc)
375 { 348 error ("for A^b, A must be square");
376 error ("for A^b, A must be square");
377 }
378 else 349 else
379 { 350 {
380 if ((int) b == b) 351 if ((int) b == b)
381 { 352 {
382 int btmp = (int) b; 353 int btmp = (int) b;
437 xpow (const ComplexMatrix& a, const Complex& b) 408 xpow (const ComplexMatrix& a, const Complex& b)
438 { 409 {
439 octave_value retval; 410 octave_value retval;
440 411
441 int nr = a.rows (); 412 int nr = a.rows ();
442 int nc = a.columns (); 413 int nc = a.cols ();
443 414
444 if (nr == 0 || nc == 0 || nr != nc) 415 if (nr == 0 || nc == 0 || nr != nc)
445 { 416 error ("for A^b, A must be square");
446 error ("for A^b, A must be square");
447 }
448 else 417 else
449 { 418 {
450 EIG a_eig (a); 419 EIG a_eig (a);
451 ComplexColumnVector lambda (a_eig.eigenvalues ()); 420 ComplexColumnVector lambda (a_eig.eigenvalues ());
452 ComplexMatrix Q (a_eig.eigenvectors ()); 421 ComplexMatrix Q (a_eig.eigenvectors ());
482 elem_xpow (double a, const Matrix& b) 451 elem_xpow (double a, const Matrix& b)
483 { 452 {
484 octave_value retval; 453 octave_value retval;
485 454
486 int nr = b.rows (); 455 int nr = b.rows ();
487 int nc = b.columns (); 456 int nc = b.cols ();
488 457
489 // For now, assume the worst. 458 // For now, assume the worst.
490 459
491 if (a < 0.0) 460 if (a < 0.0)
492 { 461 {
514 // -*- 2 -*- 483 // -*- 2 -*-
515 octave_value 484 octave_value
516 elem_xpow (double a, const ComplexMatrix& b) 485 elem_xpow (double a, const ComplexMatrix& b)
517 { 486 {
518 int nr = b.rows (); 487 int nr = b.rows ();
519 int nc = b.columns (); 488 int nc = b.cols ();
520 489
521 ComplexMatrix result (nr, nc); 490 ComplexMatrix result (nr, nc);
522 for (int j = 0; j < nc; j++) 491 for (int j = 0; j < nc; j++)
523 for (int i = 0; i < nr; i++) 492 for (int i = 0; i < nr; i++)
524 result (i, j) = pow (a, b (i, j)); 493 result (i, j) = pow (a, b (i, j));
531 elem_xpow (const Matrix& a, double b) 500 elem_xpow (const Matrix& a, double b)
532 { 501 {
533 octave_value retval; 502 octave_value retval;
534 503
535 int nr = a.rows (); 504 int nr = a.rows ();
536 int nc = a.columns (); 505 int nc = a.cols ();
537 506
538 if ((int) b != b && any_element_is_negative (a)) 507 if ((int) b != b && a.any_element_is_negative ())
539 { 508 {
540 ComplexMatrix result (nr, nc); 509 ComplexMatrix result (nr, nc);
541 for (int j = 0; j < nc; j++) 510 for (int j = 0; j < nc; j++)
542 for (int i = 0; i < nr; i++) 511 for (int i = 0; i < nr; i++)
543 { 512 {
565 elem_xpow (const Matrix& a, const Matrix& b) 534 elem_xpow (const Matrix& a, const Matrix& b)
566 { 535 {
567 octave_value retval; 536 octave_value retval;
568 537
569 int nr = a.rows (); 538 int nr = a.rows ();
570 int nc = a.columns (); 539 int nc = a.cols ();
571 540
572 assert (nr == b.rows () && nc == b.columns ()); 541 int b_nr = b.rows ();
542 int b_nc = b.cols ();
543
544 if (nr != b_nr || nc != b_nc)
545 {
546 gripe_nonconformant ("operator .^", nr, nc, b_nr, b_nc);
547 return octave_value ();
548 }
573 549
574 int convert_to_complex = 0; 550 int convert_to_complex = 0;
575 for (int j = 0; j < nc; j++) 551 for (int j = 0; j < nc; j++)
576 for (int i = 0; i < nr; i++) 552 for (int i = 0; i < nr; i++)
577 { 553 {
582 convert_to_complex = 1; 558 convert_to_complex = 1;
583 goto done; 559 goto done;
584 } 560 }
585 } 561 }
586 562
587 done: 563 done:
588 564
589 if (convert_to_complex) 565 if (convert_to_complex)
590 { 566 {
591 ComplexMatrix complex_result (nr, nc); 567 ComplexMatrix complex_result (nr, nc);
592 568
617 // -*- 5 -*- 593 // -*- 5 -*-
618 octave_value 594 octave_value
619 elem_xpow (const Matrix& a, const Complex& b) 595 elem_xpow (const Matrix& a, const Complex& b)
620 { 596 {
621 int nr = a.rows (); 597 int nr = a.rows ();
622 int nc = a.columns (); 598 int nc = a.cols ();
623 599
624 ComplexMatrix result (nr, nc); 600 ComplexMatrix result (nr, nc);
625 for (int j = 0; j < nc; j++) 601 for (int j = 0; j < nc; j++)
626 for (int i = 0; i < nr; i++) 602 for (int i = 0; i < nr; i++)
627 result (i, j) = pow (a (i, j), b); 603 result (i, j) = pow (a (i, j), b);
632 // -*- 6 -*- 608 // -*- 6 -*-
633 octave_value 609 octave_value
634 elem_xpow (const Matrix& a, const ComplexMatrix& b) 610 elem_xpow (const Matrix& a, const ComplexMatrix& b)
635 { 611 {
636 int nr = a.rows (); 612 int nr = a.rows ();
637 int nc = a.columns (); 613 int nc = a.cols ();
638 614
639 assert (nr == b.rows () && nc == b.columns ()); 615 int b_nr = b.rows ();
616 int b_nc = b.cols ();
617
618 if (nr != b_nr || nc != b_nc)
619 {
620 gripe_nonconformant ("operator .^", nr, nc, b_nr, b_nc);
621 return octave_value ();
622 }
640 623
641 ComplexMatrix result (nr, nc); 624 ComplexMatrix result (nr, nc);
642 for (int j = 0; j < nc; j++) 625 for (int j = 0; j < nc; j++)
643 for (int i = 0; i < nr; i++) 626 for (int i = 0; i < nr; i++)
644 result (i, j) = pow (a (i, j), b (i, j)); 627 result (i, j) = pow (a (i, j), b (i, j));
649 // -*- 7 -*- 632 // -*- 7 -*-
650 octave_value 633 octave_value
651 elem_xpow (const Complex& a, const Matrix& b) 634 elem_xpow (const Complex& a, const Matrix& b)
652 { 635 {
653 int nr = b.rows (); 636 int nr = b.rows ();
654 int nc = b.columns (); 637 int nc = b.cols ();
655 638
656 ComplexMatrix result (nr, nc); 639 ComplexMatrix result (nr, nc);
657 for (int j = 0; j < nc; j++) 640 for (int j = 0; j < nc; j++)
658 for (int i = 0; i < nr; i++) 641 for (int i = 0; i < nr; i++)
659 { 642 {
670 // -*- 8 -*- 653 // -*- 8 -*-
671 octave_value 654 octave_value
672 elem_xpow (const Complex& a, const ComplexMatrix& b) 655 elem_xpow (const Complex& a, const ComplexMatrix& b)
673 { 656 {
674 int nr = b.rows (); 657 int nr = b.rows ();
675 int nc = b.columns (); 658 int nc = b.cols ();
676 659
677 ComplexMatrix result (nr, nc); 660 ComplexMatrix result (nr, nc);
678 for (int j = 0; j < nc; j++) 661 for (int j = 0; j < nc; j++)
679 for (int i = 0; i < nr; i++) 662 for (int i = 0; i < nr; i++)
680 result (i, j) = pow (a, b (i, j)); 663 result (i, j) = pow (a, b (i, j));
685 // -*- 9 -*- 668 // -*- 9 -*-
686 octave_value 669 octave_value
687 elem_xpow (const ComplexMatrix& a, double b) 670 elem_xpow (const ComplexMatrix& a, double b)
688 { 671 {
689 int nr = a.rows (); 672 int nr = a.rows ();
690 int nc = a.columns (); 673 int nc = a.cols ();
691 674
692 ComplexMatrix result (nr, nc); 675 ComplexMatrix result (nr, nc);
693 676
694 if (xisint (b)) 677 if (xisint (b))
695 { 678 {
710 // -*- 10 -*- 693 // -*- 10 -*-
711 octave_value 694 octave_value
712 elem_xpow (const ComplexMatrix& a, const Matrix& b) 695 elem_xpow (const ComplexMatrix& a, const Matrix& b)
713 { 696 {
714 int nr = a.rows (); 697 int nr = a.rows ();
715 int nc = a.columns (); 698 int nc = a.cols ();
716 699
717 assert (nr == b.rows () && nc == b.columns ()); 700 int b_nr = b.rows ();
701 int b_nc = b.cols ();
702
703 if (nr != b_nr || nc != b_nc)
704 {
705 gripe_nonconformant ("operator .^", nr, nc, b_nr, b_nc);
706 return octave_value ();
707 }
718 708
719 ComplexMatrix result (nr, nc); 709 ComplexMatrix result (nr, nc);
720 for (int j = 0; j < nc; j++) 710 for (int j = 0; j < nc; j++)
721 for (int i = 0; i < nr; i++) 711 for (int i = 0; i < nr; i++)
722 { 712 {
733 // -*- 11 -*- 723 // -*- 11 -*-
734 octave_value 724 octave_value
735 elem_xpow (const ComplexMatrix& a, const Complex& b) 725 elem_xpow (const ComplexMatrix& a, const Complex& b)
736 { 726 {
737 int nr = a.rows (); 727 int nr = a.rows ();
738 int nc = a.columns (); 728 int nc = a.cols ();
739 729
740 ComplexMatrix result (nr, nc); 730 ComplexMatrix result (nr, nc);
741 for (int j = 0; j < nc; j++) 731 for (int j = 0; j < nc; j++)
742 for (int i = 0; i < nr; i++) 732 for (int i = 0; i < nr; i++)
743 result (i, j) = pow (a (i, j), b); 733 result (i, j) = pow (a (i, j), b);
748 // -*- 12 -*- 738 // -*- 12 -*-
749 octave_value 739 octave_value
750 elem_xpow (const ComplexMatrix& a, const ComplexMatrix& b) 740 elem_xpow (const ComplexMatrix& a, const ComplexMatrix& b)
751 { 741 {
752 int nr = a.rows (); 742 int nr = a.rows ();
753 int nc = a.columns (); 743 int nc = a.cols ();
754 744
755 ComplexMatrix result (nr, nc); 745 int b_nr = b.rows ();
756 746 int b_nc = b.cols ();
747
748 if (nr != b_nr || nc != b_nc)
749 {
750 gripe_nonconformant ("operator .^", nr, nc, b_nr, b_nc);
751 return octave_value ();
752 }
753
754 ComplexMatrix result (nr, nc);
757 for (int j = 0; j < nc; j++) 755 for (int j = 0; j < nc; j++)
758 for (int i = 0; i < nr; i++) 756 for (int i = 0; i < nr; i++)
759 result (i, j) = pow (a (i, j), b (i, j)); 757 result (i, j) = pow (a (i, j), b (i, j));
760 758
761 return result; 759 return result;