Mercurial > hg > octave-nkf

comparison src/arith-ops.cc @ 1:78fd87e624cb

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 1993-08-08 01:13:40 by jwe] Initial revision
author jwe
date Sun, 08 Aug 1993 01:13:40 +0000
parents
children 7849db4b6dbc
comparison
equal deleted inserted replaced
0:22412e3a4641 1:78fd87e624cb
1 // Helper functions for arithmetic operations. -*- C++ -*-
2 // Used by the tree class.
3 /*
4
5 Copyright (C) 1992, 1993 John W. Eaton
6
7 This file is part of Octave.
8
9 Octave is free software; you can redistribute it and/or modify it
10 under the terms of the GNU General Public License as published by the
11 Free Software Foundation; either version 2, or (at your option) any
12 later version.
13
14 Octave is distributed in the hope that it will be useful, but WITHOUT
15 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
16 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
17 for more details.
18
19 You should have received a copy of the GNU General Public License
20 along with Octave; see the file COPYING. If not, write to the Free
21 Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
22
23 */
24
25 #ifdef __GNUG__
26 #pragma implementation
27 #endif
28
29 #include <ctype.h>
30 #include <setjmp.h>
31 #include <math.h>
32
33 #include "error.h"
34 #include "gripes.h"
35 #include "utils.h"
36 #include "mappers.h"
37 #include "user-prefs.h"
38 #include "tree-const.h"
39 #include "arith-ops.h"
40 #include "unwind-prot.h"
41 #include "xpow.h"
42 #include "xdiv.h"
43
44 #if defined (HAVE_ISINF) || (defined (HAVE_FINITE) && defined (HAVE_ISNAN))
45 #define DIVIDE_BY_ZERO_ERROR \
46 do \
47 { \
48 if (user_pref.warn_divide_by_zero) \
49 warning ("division by zero"); \
50 } \
51 while (0)
52 #else
53 #define DIVIDE_BY_ZERO_ERROR \
54 do \
55 { \
56 error ("division by zero attempted"); \
57 return tree_constant (); \
58 } \
59 while (0)
60 #endif
61
62 // But first, some stupid functions that don\'t deserve to be in the
63 // Matrix class...
64
65 enum
66 Matrix_bool_op
67 {
68 Matrix_LT,
69 Matrix_LE,
70 Matrix_EQ,
71 Matrix_GE,
72 Matrix_GT,
73 Matrix_NE,
74 Matrix_AND,
75 Matrix_OR,
76 };
77
78 /*
79 * Stupid binary comparison operations like the ones Matlab provides.
80 * One for each type combination, in the order given here:
81 *
82 * op2 \ op1: s m cs cm
83 * +-- +---+---+----+----+
84 * scalar | | * | 3 | * | 9 |
85 * +---+---+----+----+
86 * matrix | 1 | 4 | 7 | 10 |
87 * +---+---+----+----+
88 * complex_scalar | * | 5 | * | 11 |
89 * +---+---+----+----+
90 * complex_matrix | 2 | 6 | 8 | 12 |
91 * +---+---+----+----+
92 */
93
94 /* 1 */
95 static Matrix
96 mx_stupid_bool_op (Matrix_bool_op op, double s, Matrix& a)
97 {
98 int ar = a.rows ();
99 int ac = a.columns ();
100 Matrix t (ar, ac);
101 for (int j = 0; j < ac; j++)
102 for (int i = 0; i < ar; i++)
103 {
104 switch (op)
105 {
106 case Matrix_LT:
107 t.elem (i,j) = s < a.elem (i,j);
108 break;
109 case Matrix_LE:
110 t.elem (i,j) = s <= a.elem (i,j);
111 break;
112 case Matrix_EQ:
113 t.elem (i,j) = s == a.elem (i,j);
114 break;
115 case Matrix_GE:
116 t.elem (i,j) = s >= a.elem (i,j);
117 break;
118 case Matrix_GT:
119 t.elem (i,j) = s > a.elem (i,j);
120 break;
121 case Matrix_NE:
122 t.elem (i,j) = s != a.elem (i,j);
123 break;
124 case Matrix_AND:
125 t.elem (i,j) = s && a.elem (i,j);
126 break;
127 case Matrix_OR:
128 t.elem (i,j) = s || a.elem (i,j);
129 break;
130 default:
131 panic_impossible ();
132 break;
133 }
134 }
135 return t;
136 }
137
138 /* 2 */
139 static Matrix
140 mx_stupid_bool_op (Matrix_bool_op op, double s, ComplexMatrix& a)
141 {
142 int ar = a.rows ();
143 int ac = a.columns ();
144 Matrix t (ar, ac);
145 for (int j = 0; j < ac; j++)
146 for (int i = 0; i < ar; i++)
147 {
148 switch (op)
149 {
150 case Matrix_LT:
151 t.elem (i,j) = s < real (a.elem (i,j));
152 break;
153 case Matrix_LE:
154 t.elem (i,j) = s <= real (a.elem (i,j));
155 break;
156 case Matrix_EQ:
157 t.elem (i,j) = s == a.elem (i,j);
158 break;
159 case Matrix_GE:
160 t.elem (i,j) = s >= real (a.elem (i,j));
161 break;
162 case Matrix_GT:
163 t.elem (i,j) = s > real (a.elem (i,j));
164 break;
165 case Matrix_NE:
166 t.elem (i,j) = s != a.elem (i,j);
167 break;
168 case Matrix_AND:
169 t.elem (i,j) = s && (a.elem (i,j) != 0.0);
170 break;
171 case Matrix_OR:
172 t.elem (i,j) = s || (a.elem (i,j) != 0.0);
173 break;
174 default:
175 panic_impossible ();
176 break;
177 }
178 }
179 return t;
180 }
181
182 /* 3 */
183 static Matrix
184 mx_stupid_bool_op (Matrix_bool_op op, Matrix& a, double s)
185 {
186 int ar = a.rows ();
187 int ac = a.columns ();
188 Matrix t (ar, ac);
189 for (int j = 0; j < ac; j++)
190 for (int i = 0; i < ar; i++)
191 {
192 switch (op)
193 {
194 case Matrix_LT:
195 t.elem (i,j) = a.elem (i,j) < s;
196 break;
197 case Matrix_LE:
198 t.elem (i,j) = a.elem (i,j) <= s;
199 break;
200 case Matrix_EQ:
201 t.elem (i,j) = a.elem (i,j) == s;
202 break;
203 case Matrix_GE:
204 t.elem (i,j) = a.elem (i,j) >= s;
205 break;
206 case Matrix_GT:
207 t.elem (i,j) = a.elem (i,j) > s;
208 break;
209 case Matrix_NE:
210 t.elem (i,j) = a.elem (i,j) != s;
211 break;
212 case Matrix_AND:
213 t.elem (i,j) = a.elem (i,j) && s;
214 break;
215 case Matrix_OR:
216 t.elem (i,j) = a.elem (i,j) || s;
217 break;
218 default:
219 panic_impossible ();
220 break;
221 }
222 }
223 return t;
224 }
225
226 /* 4 */
227 static Matrix
228 mx_stupid_bool_op (Matrix_bool_op op, Matrix& a, Complex& s)
229 {
230 int ar = a.rows ();
231 int ac = a.columns ();
232 Matrix t (ar, ac);
233 for (int j = 0; j < ac; j++)
234 for (int i = 0; i < ar; i++)
235 {
236 switch (op)
237 {
238 case Matrix_LT:
239 t.elem (i,j) = a.elem (i,j) < real (s);
240 break;
241 case Matrix_LE:
242 t.elem (i,j) = a.elem (i,j) <= real (s);
243 break;
244 case Matrix_EQ:
245 t.elem (i,j) = a.elem (i,j) == s;
246 break;
247 case Matrix_GE:
248 t.elem (i,j) = a.elem (i,j) >= real (s);
249 break;
250 case Matrix_GT:
251 t.elem (i,j) = a.elem (i,j) > real (s);
252 break;
253 case Matrix_NE:
254 t.elem (i,j) = a.elem (i,j) != s;
255 break;
256 case Matrix_AND:
257 t.elem (i,j) = a.elem (i,j) && (s != 0.0);
258 break;
259 case Matrix_OR:
260 t.elem (i,j) = a.elem (i,j) || (s != 0.0);
261 break;
262 default:
263 panic_impossible ();
264 break;
265 }
266 }
267 return t;
268 }
269
270 /* 5 */
271 static Matrix
272 mx_stupid_bool_op (Matrix_bool_op op, Matrix& a, Matrix& b)
273 {
274 int ar = a.rows ();
275 int ac = a.columns ();
276
277 if (ar != b.rows () || ac != b.columns ())
278 {
279 gripe_nonconformant ();
280 jump_to_top_level ();
281 }
282
283 Matrix c (ar, ac);
284
285 for (int j = 0; j < ac; j++)
286 for (int i = 0; i < ar; i++)
287 {
288 switch (op)
289 {
290 case Matrix_LT:
291 c.elem (i, j) = a.elem (i, j) < b.elem (i, j);
292 break;
293 case Matrix_LE:
294 c.elem (i, j) = a.elem (i, j) <= b.elem (i, j);
295 break;
296 case Matrix_EQ:
297 c.elem (i, j) = a.elem (i, j) == b.elem (i, j);
298 break;
299 case Matrix_GE:
300 c.elem (i, j) = a.elem (i, j) >= b.elem (i, j);
301 break;
302 case Matrix_GT:
303 c.elem (i, j) = a.elem (i, j) > b.elem (i, j);
304 break;
305 case Matrix_NE:
306 c.elem (i, j) = a.elem (i, j) != b.elem (i, j);
307 break;
308 case Matrix_AND:
309 c.elem (i, j) = a.elem (i, j) && b.elem (i, j);
310 break;
311 case Matrix_OR:
312 c.elem (i, j) = a.elem (i, j) || b.elem (i, j);
313 break;
314 default:
315 panic_impossible ();
316 break;
317 }
318 }
319 return c;
320 }
321
322 /* 6 */
323 static Matrix
324 mx_stupid_bool_op (Matrix_bool_op op, Matrix& a, ComplexMatrix& b)
325 {
326 int ar = a.rows ();
327 int ac = a.columns ();
328
329 if (ar != b.rows () || ac != b.columns ())
330 {
331 gripe_nonconformant ();
332 jump_to_top_level ();
333 }
334
335 Matrix c (ar, ac);
336
337 for (int j = 0; j < ac; j++)
338 for (int i = 0; i < ar; i++)
339 {
340 switch (op)
341 {
342 case Matrix_LT:
343 c.elem (i, j) = a.elem (i, j) < real (b.elem (i, j));
344 break;
345 case Matrix_LE:
346 c.elem (i, j) = a.elem (i, j) <= real (b.elem (i, j));
347 break;
348 case Matrix_EQ:
349 c.elem (i, j) = a.elem (i, j) == b.elem (i, j);
350 break;
351 case Matrix_GE:
352 c.elem (i, j) = a.elem (i, j) >= real (b.elem (i, j));
353 break;
354 case Matrix_GT:
355 c.elem (i, j) = a.elem (i, j) > real (b.elem (i, j));
356 break;
357 case Matrix_NE:
358 c.elem (i, j) = a.elem (i, j) != b.elem (i, j);
359 break;
360 case Matrix_AND:
361 c.elem (i, j) = a.elem (i, j) && (b.elem (i, j) != 0.0);
362 break;
363 case Matrix_OR:
364 c.elem (i, j) = a.elem (i, j) || (b.elem (i, j) != 0.0);
365 break;
366 default:
367 panic_impossible ();
368 break;
369 }
370 }
371 return c;
372 }
373
374 /* 7 */
375 static Matrix
376 mx_stupid_bool_op (Matrix_bool_op op, Complex& s, Matrix& a)
377 {
378 int ar = a.rows ();
379 int ac = a.columns ();
380 Matrix t (ar, ac);
381 for (int j = 0; j < ac; j++)
382 for (int i = 0; i < ar; i++)
383 {
384 switch (op)
385 {
386 case Matrix_LT:
387 t.elem (i,j) = real (s) < a.elem (i,j);
388 break;
389 case Matrix_LE:
390 t.elem (i,j) = real (s) <= a.elem (i,j);
391 break;
392 case Matrix_EQ:
393 t.elem (i,j) = s == a.elem (i,j);
394 break;
395 case Matrix_GE:
396 t.elem (i,j) = real (s) >= a.elem (i,j);
397 break;
398 case Matrix_GT:
399 t.elem (i,j) = real (s) > a.elem (i,j);
400 break;
401 case Matrix_NE:
402 t.elem (i,j) = s != a.elem (i,j);
403 break;
404 case Matrix_AND:
405 t.elem (i,j) = (s != 0.0) && a.elem (i,j);
406 break;
407 case Matrix_OR:
408 t.elem (i,j) = (s != 0.0) || a.elem (i,j);
409 break;
410 default:
411 panic_impossible ();
412 break;
413 }
414 }
415 return t;
416 }
417
418 /* 8 */
419 static Matrix
420 mx_stupid_bool_op (Matrix_bool_op op, Complex& s, ComplexMatrix& a)
421 {
422 int ar = a.rows ();
423 int ac = a.columns ();
424 Matrix t (ar, ac);
425 for (int j = 0; j < ac; j++)
426 for (int i = 0; i < ar; i++)
427 {
428 switch (op)
429 {
430 case Matrix_LT:
431 t.elem (i,j) = real (s) < real (a.elem (i,j));
432 break;
433 case Matrix_LE:
434 t.elem (i,j) = real (s) <= real (a.elem (i,j));
435 break;
436 case Matrix_EQ:
437 t.elem (i,j) = s == a.elem (i,j);
438 break;
439 case Matrix_GE:
440 t.elem (i,j) = real (s) >= real (a.elem (i,j));
441 break;
442 case Matrix_GT:
443 t.elem (i,j) = real (s) > real (a.elem (i,j));
444 break;
445 case Matrix_NE:
446 t.elem (i,j) = s != a.elem (i,j);
447 break;
448 case Matrix_AND:
449 t.elem (i,j) = (s != 0.0) && (a.elem (i,j) != 0.0);
450 break;
451 case Matrix_OR:
452 t.elem (i,j) = (s != 0.0) || (a.elem (i,j) != 0.0);
453 break;
454 default:
455 panic_impossible ();
456 break;
457 }
458 }
459 return t;
460 }
461
462 /* 9 */
463 static Matrix
464 mx_stupid_bool_op (Matrix_bool_op op, ComplexMatrix& a, double s)
465 {
466 int ar = a.rows ();
467 int ac = a.columns ();
468 Matrix t (ar, ac);
469 for (int j = 0; j < ac; j++)
470 for (int i = 0; i < ar; i++)
471 {
472 switch (op)
473 {
474 case Matrix_LT:
475 t.elem (i,j) = real (a.elem (i,j)) < s;
476 break;
477 case Matrix_LE:
478 t.elem (i,j) = real (a.elem (i,j)) <= s;
479 break;
480 case Matrix_EQ:
481 t.elem (i,j) = a.elem (i,j) == s;
482 break;
483 case Matrix_GE:
484 t.elem (i,j) = real (a.elem (i,j)) >= s;
485 break;
486 case Matrix_GT:
487 t.elem (i,j) = real (a.elem (i,j)) > s;
488 break;
489 case Matrix_NE:
490 t.elem (i,j) = a.elem (i,j) != s;
491 break;
492 case Matrix_AND:
493 t.elem (i,j) = (a.elem (i,j) != 0.0) && s;
494 break;
495 case Matrix_OR:
496 t.elem (i,j) = (a.elem (i,j) != 0.0) || s;
497 break;
498 default:
499 panic_impossible ();
500 break;
501 }
502 }
503 return t;
504 }
505
506 /* 10 */
507 static Matrix
508 mx_stupid_bool_op (Matrix_bool_op op, ComplexMatrix& a, Complex& s)
509 {
510 int ar = a.rows ();
511 int ac = a.columns ();
512 Matrix t (ar, ac);
513 for (int j = 0; j < ac; j++)
514 for (int i = 0; i < ar; i++)
515 {
516 switch (op)
517 {
518 case Matrix_LT:
519 t.elem (i,j) = real (a.elem (i,j)) < real (s);
520 break;
521 case Matrix_LE:
522 t.elem (i,j) = real (a.elem (i,j)) <= real (s);
523 break;
524 case Matrix_EQ:
525 t.elem (i,j) = a.elem (i,j) == s;
526 break;
527 case Matrix_GE:
528 t.elem (i,j) = real (a.elem (i,j)) >= real (s);
529 break;
530 case Matrix_GT:
531 t.elem (i,j) = real (a.elem (i,j)) > real (s);
532 break;
533 case Matrix_NE:
534 t.elem (i,j) = a.elem (i,j) != s;
535 break;
536 case Matrix_AND:
537 t.elem (i,j) = (a.elem (i,j) != 0.0) && (s != 0.0);
538 break;
539 case Matrix_OR:
540 t.elem (i,j) = (a.elem (i,j) != 0.0) || (s != 0.0);
541 break;
542 default:
543 panic_impossible ();
544 break;
545 }
546 }
547 return t;
548 }
549
550 /* 11 */
551 static Matrix
552 mx_stupid_bool_op (Matrix_bool_op op, ComplexMatrix& a, Matrix& b)
553 {
554 int ar = a.rows ();
555 int ac = a.columns ();
556
557 if (ar != b.rows () || ac != b.columns ())
558 {
559 gripe_nonconformant ();
560 jump_to_top_level ();
561 }
562
563 Matrix c (ar, ac);
564
565 for (int j = 0; j < ac; j++)
566 for (int i = 0; i < ar; i++)
567 {
568 switch (op)
569 {
570 case Matrix_LT:
571 c.elem (i, j) = real (a.elem (i, j)) < b.elem (i, j);
572 break;
573 case Matrix_LE:
574 c.elem (i, j) = real (a.elem (i, j)) <= b.elem (i, j);
575 break;
576 case Matrix_EQ:
577 c.elem (i, j) = a.elem (i, j) == b.elem (i, j);
578 break;
579 case Matrix_GE:
580 c.elem (i, j) = real (a.elem (i, j)) >= b.elem (i, j);
581 break;
582 case Matrix_GT:
583 c.elem (i, j) = real (a.elem (i, j)) > b.elem (i, j);
584 break;
585 case Matrix_NE:
586 c.elem (i, j) = a.elem (i, j) != b.elem (i, j);
587 break;
588 case Matrix_AND:
589 c.elem (i, j) = (a.elem (i, j) != 0.0) && b.elem (i, j);
590 break;
591 case Matrix_OR:
592 c.elem (i, j) = (a.elem (i, j) != 0.0) || b.elem (i, j);
593 break;
594 default:
595 panic_impossible ();
596 break;
597 }
598 }
599 return c;
600 }
601
602 /* 12 */
603 static Matrix
604 mx_stupid_bool_op (Matrix_bool_op op, ComplexMatrix& a, ComplexMatrix& b)
605 {
606 int ar = a.rows ();
607 int ac = a.columns ();
608
609 if (ar != b.rows () || ac != b.columns ())
610 {
611 gripe_nonconformant ();
612 jump_to_top_level ();
613 }
614
615 Matrix c (ar, ac);
616
617 for (int j = 0; j < ac; j++)
618 for (int i = 0; i < ar; i++)
619 {
620 switch (op)
621 {
622 case Matrix_LT:
623 c.elem (i, j) = real (a.elem (i, j)) < real (b.elem (i, j));
624 break;
625 case Matrix_LE:
626 c.elem (i, j) = real (a.elem (i, j)) <= real (b.elem (i, j));
627 break;
628 case Matrix_EQ:
629 c.elem (i, j) = a.elem (i, j) == b.elem (i, j);
630 break;
631 case Matrix_GE:
632 c.elem (i, j) = real (a.elem (i, j)) >= real (b.elem (i, j));
633 break;
634 case Matrix_GT:
635 c.elem (i, j) = real (a.elem (i, j)) > real (b.elem (i, j));
636 break;
637 case Matrix_NE:
638 c.elem (i, j) = a.elem (i, j) != b.elem (i, j);
639 break;
640 case Matrix_AND:
641 c.elem (i, j) = (a.elem (i, j) != 0.0) && (b.elem (i, j) != 0.0);
642 break;
643 case Matrix_OR:
644 c.elem (i, j) = (a.elem (i, j) != 0.0) || (b.elem (i, j) != 0.0);
645 break;
646 default:
647 panic_impossible ();
648 break;
649 }
650 }
651 return c;
652 }
653
654 /*
655 * Check row and column dimensions for binary matrix operations.
656 */
657 static inline int
658 m_add_conform (Matrix& m1, Matrix& m2, int warn)
659 {
660 int ok = (m1.rows () == m2.rows () && m1.columns () == m2.columns ());
661 if (!ok && warn)
662 gripe_nonconformant ();
663 return ok;
664 }
665
666 static inline int
667 m_add_conform (Matrix& m1, ComplexMatrix& m2, int warn)
668 {
669 int ok = (m1.rows () == m2.rows () && m1.columns () == m2.columns ());
670 if (!ok && warn)
671 gripe_nonconformant ();
672 return ok;
673 }
674
675 static inline int
676 m_add_conform (ComplexMatrix& m1, Matrix& m2, int warn)
677 {
678 int ok = (m1.rows () == m2.rows () && m1.columns () == m2.columns ());
679 if (!ok && warn)
680 gripe_nonconformant ();
681 return ok;
682 }
683
684 static inline int
685 m_add_conform (ComplexMatrix& m1, ComplexMatrix& m2, int warn)
686 {
687 int ok = (m1.rows () == m2.rows () && m1.columns () == m2.columns ());
688 if (!ok && warn)
689 gripe_nonconformant ();
690 return ok;
691 }
692
693 static inline int
694 m_mul_conform (Matrix& m1, Matrix& m2, int warn)
695 {
696 int ok = (m1.columns () == m2.rows ());
697 if (!ok && warn)
698 gripe_nonconformant ();
699 return ok;
700 }
701
702 static inline int
703 m_mul_conform (Matrix& m1, ComplexMatrix& m2, int warn)
704 {
705 int ok = (m1.columns () == m2.rows ());
706 if (!ok && warn)
707 gripe_nonconformant ();
708 return ok;
709 }
710
711 static inline int
712 m_mul_conform (ComplexMatrix& m1, Matrix& m2, int warn)
713 {
714 int ok = (m1.columns () == m2.rows ());
715 if (!ok && warn)
716 gripe_nonconformant ();
717 return ok;
718 }
719
720 static inline int
721 m_mul_conform (ComplexMatrix& m1, ComplexMatrix& m2, int warn)
722 {
723 int ok = (m1.columns () == m2.rows ());
724 if (!ok && warn)
725 gripe_nonconformant ();
726 return ok;
727 }
728
729 /*
730 * Unary operations. One for each numeric data type:
731 *
732 * scalar
733 * complex_scalar
734 * matrix
735 * complex_matrix
736 *
737 */
738
739 tree_constant
740 do_unary_op (double d, tree::expression_type t)
741 {
742 double result = 0.0;
743 switch (t)
744 {
745 case tree::not:
746 result = (! d);
747 break;
748 case tree::uminus:
749 result = -d;
750 break;
751 case tree::hermitian:
752 case tree::transpose:
753 result = d;
754 break;
755 default:
756 panic_impossible ();
757 break;
758 }
759
760 return tree_constant (result);
761 }
762
763 tree_constant
764 do_unary_op (Matrix& a, tree::expression_type t)
765 {
766 Matrix result;
767 switch (t)
768 {
769 case tree::not:
770 result = (! a);
771 break;
772 case tree::uminus:
773 result = -a;
774 break;
775 case tree::hermitian:
776 case tree::transpose:
777 result = a.transpose ();
778 break;
779 default:
780 panic_impossible ();
781 break;
782 }
783
784 return tree_constant (result);
785 }
786
787 tree_constant
788 do_unary_op (Complex& c, tree::expression_type t)
789 {
790 Complex result = 0.0;
791 switch (t)
792 {
793 case tree::not:
794 result = (c == 0.0);
795 break;
796 case tree::uminus:
797 result = -c;
798 break;
799 case tree::hermitian:
800 result = conj (c);
801 break;
802 case tree::transpose:
803 result = c;
804 break;
805 default:
806 panic_impossible ();
807 break;
808 }
809
810 return tree_constant (result);
811 }
812
813 tree_constant
814 do_unary_op (ComplexMatrix& a, tree::expression_type t)
815 {
816 ComplexMatrix result;
817 switch (t)
818 {
819 case tree::not:
820 result = (! a);
821 break;
822 case tree::uminus:
823 result = -a;
824 break;
825 case tree::hermitian:
826 result = a.hermitian ();
827 break;
828 case tree::transpose:
829 result = a.transpose ();
830 break;
831 default:
832 panic_impossible ();
833 break;
834 }
835
836 return tree_constant (result);
837 }
838
839 /*
840 * Binary operations. One for each type combination, in the order
841 * given here:
842 *
843 * op2 \ op1: s m cs cm
844 * +-- +---+---+----+----+
845 * scalar | | 1 | 5 | 9 | 13 |
846 * +---+---+----+----+
847 * matrix | 2 | 6 | 10 | 14 |
848 * +---+---+----+----+
849 * complex_scalar | 3 | 7 | 11 | 15 |
850 * +---+---+----+----+
851 * complex_matrix | 4 | 8 | 12 | 16 |
852 * +---+---+----+----+
853 */
854
855 /* 1 */
856 tree_constant
857 do_binary_op (double a, double b, tree::expression_type t)
858 {
859 double result = 0.0;
860 switch (t)
861 {
862 case tree::add:
863 result = a + b;
864 break;
865 case tree::subtract:
866 result = a - b;
867 break;
868 case tree::multiply:
869 case tree::el_mul:
870 result = a * b;
871 break;
872 case tree::divide:
873 case tree::el_div:
874 if (b == 0.0)
875 DIVIDE_BY_ZERO_ERROR;
876 result = a / b;
877 break;
878 case tree::leftdiv:
879 case tree::el_leftdiv:
880 if (a == 0.0)
881 DIVIDE_BY_ZERO_ERROR;
882 result = b / a;
883 break;
884 case tree::power:
885 case tree::elem_pow:
886 return xpow (a, b);
887 break;
888 case tree::cmp_lt:
889 result = a < b;
890 break;
891 case tree::cmp_le:
892 result = a <= b;
893 break;
894 case tree::cmp_eq:
895 result = a == b;
896 break;
897 case tree::cmp_ge:
898 result = a >= b;
899 break;
900 case tree::cmp_gt:
901 result = a > b;
902 break;
903 case tree::cmp_ne:
904 result = a != b;
905 break;
906 case tree::and:
907 result = (a && b);
908 break;
909 case tree::or:
910 result = (a || b);
911 break;
912 default:
913 panic_impossible ();
914 break;
915 }
916
917 return tree_constant (result);
918 }
919
920 /* 2 */
921 tree_constant
922 do_binary_op (double a, Matrix& b, tree::expression_type t)
923 {
924 Matrix result;
925 switch (t)
926 {
927 case tree::add:
928 result = a + b;
929 break;
930 case tree::subtract:
931 result = a - b;
932 break;
933 case tree::el_leftdiv:
934 case tree::leftdiv:
935 if (a == 0.0)
936 DIVIDE_BY_ZERO_ERROR;
937 a = 1.0 / a;
938 // fall through...
939 case tree::multiply:
940 case tree::el_mul:
941 result = a * b;
942 break;
943 case tree::el_div:
944 return x_el_div (a, b);
945 break;
946 case tree::divide:
947 error ("nonconformant right division");
948 return tree_constant ();
949 break;
950 case tree::power:
951 return xpow (a, b);
952 break;
953 case tree::elem_pow:
954 return elem_xpow (a, b);
955 break;
956 case tree::cmp_lt:
957 result = mx_stupid_bool_op (Matrix_LT, a, b);
958 break;
959 case tree::cmp_le:
960 result = mx_stupid_bool_op (Matrix_LE, a, b);
961 break;
962 case tree::cmp_eq:
963 result = mx_stupid_bool_op (Matrix_EQ, a, b);
964 break;
965 case tree::cmp_ge:
966 result = mx_stupid_bool_op (Matrix_GE, a, b);
967 break;
968 case tree::cmp_gt:
969 result = mx_stupid_bool_op (Matrix_GT, a, b);
970 break;
971 case tree::cmp_ne:
972 result = mx_stupid_bool_op (Matrix_NE, a, b);
973 break;
974 case tree::and:
975 result = mx_stupid_bool_op (Matrix_AND, a, b);
976 break;
977 case tree::or:
978 result = mx_stupid_bool_op (Matrix_OR, a, b);
979 break;
980 default:
981 panic_impossible ();
982 break;
983 }
984
985 return tree_constant (result);
986 }
987
988 /* 3 */
989 tree_constant
990 do_binary_op (double a, Complex& b, tree::expression_type t)
991 {
992 enum RT { RT_unknown, RT_real, RT_complex };
993 RT result_type = RT_unknown;
994
995 double result = 0.0;
996 Complex complex_result;
997 switch (t)
998 {
999 case tree::add:
1000 result_type = RT_complex;
1001 complex_result = a + b;
1002 break;
1003 case tree::subtract:
1004 result_type = RT_complex;
1005 complex_result = a - b;
1006 break;
1007 case tree::multiply:
1008 case tree::el_mul:
1009 result_type = RT_complex;
1010 complex_result = a * b;
1011 break;
1012 case tree::divide:
1013 case tree::el_div:
1014 result_type = RT_complex;
1015 if (b == 0.0)
1016 DIVIDE_BY_ZERO_ERROR;
1017 complex_result = a / b;
1018 break;
1019 case tree::leftdiv:
1020 case tree::el_leftdiv:
1021 result_type = RT_complex;
1022 if (a == 0.0)
1023 DIVIDE_BY_ZERO_ERROR;
1024 complex_result = b / a;
1025 break;
1026 case tree::power:
1027 case tree::elem_pow:
1028 return xpow (a, b);
1029 break;
1030 case tree::cmp_lt:
1031 result_type = RT_real;
1032 result = a < real (b);
1033 break;
1034 case tree::cmp_le:
1035 result_type = RT_real;
1036 result = a <= real (b);
1037 break;
1038 case tree::cmp_eq:
1039 result_type = RT_real;
1040 result = a == b;
1041 break;
1042 case tree::cmp_ge:
1043 result_type = RT_real;
1044 result = a >= real (b);
1045 break;
1046 case tree::cmp_gt:
1047 result_type = RT_real;
1048 result = a > real (b);
1049 break;
1050 case tree::cmp_ne:
1051 result_type = RT_real;
1052 result = a != b;
1053 break;
1054 case tree::and:
1055 result_type = RT_real;
1056 result = (a && (b != 0.0));
1057 break;
1058 case tree::or:
1059 result_type = RT_real;
1060 result = (a || (b != 0.0));
1061 break;
1062 default:
1063 panic_impossible ();
1064 break;
1065 }
1066
1067 assert (result_type != RT_unknown);
1068 if (result_type == RT_real)
1069 return tree_constant (result);
1070 else
1071 return tree_constant (complex_result);
1072 }
1073
1074 /* 4 */
1075 tree_constant
1076 do_binary_op (double a, ComplexMatrix& b, tree::expression_type t)
1077 {
1078 enum RT { RT_unknown, RT_real, RT_complex };
1079 RT result_type = RT_unknown;
1080
1081 Matrix result;
1082 ComplexMatrix complex_result;
1083 switch (t)
1084 {
1085 case tree::add:
1086 result_type = RT_complex;
1087 complex_result = a + b;
1088 break;
1089 case tree::subtract:
1090 result_type = RT_complex;
1091 complex_result = a - b;
1092 break;
1093 case tree::el_leftdiv:
1094 case tree::leftdiv:
1095 if (a == 0.0)
1096 DIVIDE_BY_ZERO_ERROR;
1097 a = 1.0 / a;
1098 // fall through...
1099 case tree::multiply:
1100 case tree::el_mul:
1101 result_type = RT_complex;
1102 complex_result = a * b;
1103 break;
1104 case tree::el_div:
1105 return x_el_div (a, b);
1106 break;
1107 case tree::divide:
1108 error ("nonconformant right division");
1109 return tree_constant ();
1110 break;
1111 case tree::power:
1112 return xpow (a, b);
1113 break;
1114 case tree::elem_pow:
1115 return elem_xpow (a, b);
1116 break;
1117 case tree::cmp_lt:
1118 result_type = RT_real;
1119 result = mx_stupid_bool_op (Matrix_LT, a, b);
1120 break;
1121 case tree::cmp_le:
1122 result_type = RT_real;
1123 result = mx_stupid_bool_op (Matrix_LE, a, b);
1124 break;
1125 case tree::cmp_eq:
1126 result_type = RT_real;
1127 result = mx_stupid_bool_op (Matrix_EQ, a, b);
1128 break;
1129 case tree::cmp_ge:
1130 result_type = RT_real;
1131 result = mx_stupid_bool_op (Matrix_GE, a, b);
1132 break;
1133 case tree::cmp_gt:
1134 result_type = RT_real;
1135 result = mx_stupid_bool_op (Matrix_GT, a, b);
1136 break;
1137 case tree::cmp_ne:
1138 result_type = RT_real;
1139 result = mx_stupid_bool_op (Matrix_NE, a, b);
1140 break;
1141 case tree::and:
1142 result_type = RT_real;
1143 result = mx_stupid_bool_op (Matrix_AND, a, b);
1144 break;
1145 case tree::or:
1146 result_type = RT_real;
1147 result = mx_stupid_bool_op (Matrix_OR, a, b);
1148 break;
1149 default:
1150 panic_impossible ();
1151 break;
1152 }
1153
1154 assert (result_type != RT_unknown);
1155 if (result_type == RT_real)
1156 return tree_constant (result);
1157 else
1158 return tree_constant (complex_result);
1159 }
1160
1161 /* 5 */
1162 tree_constant
1163 do_binary_op (Matrix& a, double b, tree::expression_type t)
1164 {
1165 Matrix result;
1166 switch (t)
1167 {
1168 case tree::add:
1169 result = a + b;
1170 break;
1171 case tree::subtract:
1172 result = a - b;
1173 break;
1174 case tree::multiply:
1175 case tree::el_mul:
1176 result = a * b;
1177 break;
1178 case tree::divide:
1179 case tree::el_div:
1180 result = a / b;
1181 break;
1182 case tree::el_leftdiv:
1183 return x_el_div (b, a);
1184 break;
1185 case tree::leftdiv:
1186 error ("nonconformant left division");
1187 return tree_constant ();
1188 break;
1189 case tree::power:
1190 return xpow (a, b);
1191 break;
1192 case tree::elem_pow:
1193 return elem_xpow (a, b);
1194 break;
1195 case tree::cmp_lt:
1196 result = mx_stupid_bool_op (Matrix_LT, a, b);
1197 break;
1198 case tree::cmp_le:
1199 result = mx_stupid_bool_op (Matrix_LE, a, b);
1200 break;
1201 case tree::cmp_eq:
1202 result = mx_stupid_bool_op (Matrix_EQ, a, b);
1203 break;
1204 case tree::cmp_ge:
1205 result = mx_stupid_bool_op (Matrix_GE, a, b);
1206 break;
1207 case tree::cmp_gt:
1208 result = mx_stupid_bool_op (Matrix_GT, a, b);
1209 break;
1210 case tree::cmp_ne:
1211 result = mx_stupid_bool_op (Matrix_NE, a, b);
1212 break;
1213 case tree::and:
1214 result = mx_stupid_bool_op (Matrix_AND, a, b);
1215 break;
1216 case tree::or:
1217 result = mx_stupid_bool_op (Matrix_OR, a, b);
1218 break;
1219 default:
1220 panic_impossible ();
1221 break;
1222 }
1223
1224 return tree_constant (result);
1225 }
1226
1227 /* 6 */
1228 tree_constant
1229 do_binary_op (Matrix& a, Matrix& b, tree::expression_type t)
1230 {
1231 Matrix result;
1232
1233 int error_cond = 0;
1234
1235 switch (t)
1236 {
1237 case tree::add:
1238 if (m_add_conform (a, b, 1))
1239 result = a + b;
1240 else
1241 error_cond = 1;
1242 break;
1243 case tree::subtract:
1244 if (m_add_conform (a, b, 1))
1245 result = a - b;
1246 else
1247 error_cond = 1;
1248 break;
1249 case tree::el_mul:
1250 if (m_add_conform (a, b, 1))
1251 result = a.product (b);
1252 else
1253 error_cond = 1;
1254 break;
1255 case tree::multiply:
1256 if (m_mul_conform (a, b, 1))
1257 result = a * b;
1258 else
1259 error_cond = 1;
1260 break;
1261 case tree::el_div:
1262 if (m_add_conform (a, b, 1))
1263 result = a.quotient (b);
1264 else
1265 error_cond = 1;
1266 break;
1267 case tree::el_leftdiv:
1268 if (m_add_conform (a, b, 1))
1269 result = b.quotient (a);
1270 else
1271 error_cond = 1;
1272 break;
1273 case tree::leftdiv:
1274 return xleftdiv (a, b);
1275 break;
1276 case tree::divide:
1277 return xdiv (a, b);
1278 break;
1279 case tree::power:
1280 error ("can't do A ^ B for A and B both matrices");
1281 error_cond = 1;
1282 break;
1283 case tree::elem_pow:
1284 if (m_add_conform (a, b, 1))
1285 return elem_xpow (a, b);
1286 else
1287 error_cond = 1;
1288 break;
1289 case tree::cmp_lt:
1290 if (m_add_conform (a, b, 1))
1291 result = mx_stupid_bool_op (Matrix_LT, a, b);
1292 else
1293 error_cond = 1;
1294 break;
1295 case tree::cmp_le:
1296 if (m_add_conform (a, b, 1))
1297 result = mx_stupid_bool_op (Matrix_LE, a, b);
1298 else
1299 error_cond = 1;
1300 break;
1301 case tree::cmp_eq:
1302 if (m_add_conform (a, b, 1))
1303 result = mx_stupid_bool_op (Matrix_EQ, a, b);
1304 else
1305 error_cond = 1;
1306 break;
1307 case tree::cmp_ge:
1308 if (m_add_conform (a, b, 1))
1309 result = mx_stupid_bool_op (Matrix_GE, a, b);
1310 else
1311 error_cond = 1;
1312 break;
1313 case tree::cmp_gt:
1314 if (m_add_conform (a, b, 1))
1315 result = mx_stupid_bool_op (Matrix_GT, a, b);
1316 else
1317 error_cond = 1;
1318 break;
1319 case tree::cmp_ne:
1320 if (m_add_conform (a, b, 1))
1321 result = mx_stupid_bool_op (Matrix_NE, a, b);
1322 else
1323 error_cond = 1;
1324 break;
1325 case tree::and:
1326 if (m_add_conform (a, b, 1))
1327 result = mx_stupid_bool_op (Matrix_AND, a, b);
1328 else
1329 error_cond = 1;
1330 break;
1331 case tree::or:
1332 if (m_add_conform (a, b, 1))
1333 result = mx_stupid_bool_op (Matrix_OR, a, b);
1334 else
1335 error_cond = 1;
1336 break;
1337 default:
1338 panic_impossible ();
1339 break;
1340 }
1341
1342 if (error_cond)
1343 return tree_constant ();
1344 else
1345 return tree_constant (result);
1346 }
1347
1348 /* 7 */
1349 tree_constant
1350 do_binary_op (Matrix& a, Complex& b, tree::expression_type t)
1351 {
1352 enum RT { RT_unknown, RT_real, RT_complex };
1353 RT result_type = RT_unknown;
1354
1355 Matrix result;
1356 ComplexMatrix complex_result;
1357 switch (t)
1358 {
1359 case tree::add:
1360 result_type = RT_complex;
1361 complex_result = a + b;
1362 break;
1363 case tree::subtract:
1364 result_type = RT_complex;
1365 complex_result = a - b;
1366 break;
1367 case tree::multiply:
1368 case tree::el_mul:
1369 result_type = RT_complex;
1370 complex_result = a * b;
1371 break;
1372 case tree::divide:
1373 case tree::el_div:
1374 result_type = RT_complex;
1375 complex_result = a / b;
1376 break;
1377 case tree::el_leftdiv:
1378 return x_el_div (b, a);
1379 break;
1380 case tree::leftdiv:
1381 error ("nonconformant left division");
1382 return tree_constant ();
1383 break;
1384 case tree::power:
1385 return xpow (a, b);
1386 break;
1387 case tree::elem_pow:
1388 return elem_xpow (a, b);
1389 break;
1390 case tree::cmp_lt:
1391 result_type = RT_real;
1392 result = mx_stupid_bool_op (Matrix_LT, a, b);
1393 break;
1394 case tree::cmp_le:
1395 result_type = RT_real;
1396 result = mx_stupid_bool_op (Matrix_LE, a, b);
1397 break;
1398 case tree::cmp_eq:
1399 result_type = RT_real;
1400 result = mx_stupid_bool_op (Matrix_EQ, a, b);
1401 break;
1402 case tree::cmp_ge:
1403 result_type = RT_real;
1404 result = mx_stupid_bool_op (Matrix_GE, a, b);
1405 break;
1406 case tree::cmp_gt:
1407 result_type = RT_real;
1408 result = mx_stupid_bool_op (Matrix_GT, a, b);
1409 break;
1410 case tree::cmp_ne:
1411 result_type = RT_real;
1412 result = mx_stupid_bool_op (Matrix_NE, a, b);
1413 break;
1414 case tree::and:
1415 result_type = RT_real;
1416 result = mx_stupid_bool_op (Matrix_AND, a, b);
1417 break;
1418 case tree::or:
1419 result_type = RT_real;
1420 result = mx_stupid_bool_op (Matrix_OR, a, b);
1421 break;
1422 default:
1423 panic_impossible ();
1424 break;
1425 }
1426
1427 assert (result_type != RT_unknown);
1428 if (result_type == RT_real)
1429 return tree_constant (result);
1430 else
1431 return tree_constant (complex_result);
1432 }
1433
1434 /* 8 */
1435 tree_constant
1436 do_binary_op (Matrix& a, ComplexMatrix& b, tree::expression_type t)
1437 {
1438 enum RT { RT_unknown, RT_real, RT_complex };
1439 RT result_type = RT_unknown;
1440
1441 Matrix result;
1442 ComplexMatrix complex_result;
1443 switch (t)
1444 {
1445 case tree::add:
1446 result_type = RT_complex;
1447 if (m_add_conform (a, b, 1))
1448 complex_result = a + b;
1449 else
1450 return tree_constant ();
1451 break;
1452 case tree::subtract:
1453 result_type = RT_complex;
1454 if (m_add_conform (a, b, 1))
1455 complex_result = a - b;
1456 else
1457 return tree_constant ();
1458 break;
1459 case tree::el_mul:
1460 result_type = RT_complex;
1461 if (m_add_conform (a, b, 1))
1462 complex_result = a.product (b);
1463 else
1464 return tree_constant ();
1465 break;
1466 case tree::multiply:
1467 result_type = RT_complex;
1468 if (m_mul_conform (a, b, 1))
1469 complex_result = a * b;
1470 else
1471 return tree_constant ();
1472 break;
1473 case tree::el_div:
1474 result_type = RT_complex;
1475 if (m_add_conform (a, b, 1))
1476 complex_result = a.quotient (b);
1477 else
1478 return tree_constant ();
1479 break;
1480 case tree::el_leftdiv:
1481 result_type = RT_complex;
1482 if (m_add_conform (a, b, 1))
1483 complex_result = b.quotient (a);
1484 else
1485 return tree_constant ();
1486 break;
1487 case tree::leftdiv:
1488 return xleftdiv (a, b);
1489 break;
1490 case tree::divide:
1491 return xdiv (a, b);
1492 break;
1493 case tree::power:
1494 error ("can't do A ^ B for A and B both matrices");
1495 return tree_constant ();
1496 break;
1497 case tree::elem_pow:
1498 if (m_add_conform (a, b, 1))
1499 return elem_xpow (a, b);
1500 else
1501 return tree_constant ();
1502 break;
1503 case tree::cmp_lt:
1504 result_type = RT_real;
1505 if (m_add_conform (a, b, 1))
1506 result = mx_stupid_bool_op (Matrix_LT, a, b);
1507 else
1508 return tree_constant ();
1509 break;
1510 case tree::cmp_le:
1511 result_type = RT_real;
1512 if (m_add_conform (a, b, 1))
1513 result = mx_stupid_bool_op (Matrix_LE, a, b);
1514 else
1515 return tree_constant ();
1516 break;
1517 case tree::cmp_eq:
1518 result_type = RT_real;
1519 if (m_add_conform (a, b, 1))
1520 result = mx_stupid_bool_op (Matrix_EQ, a, b);
1521 else
1522 return tree_constant ();
1523 break;
1524 case tree::cmp_ge:
1525 result_type = RT_real;
1526 if (m_add_conform (a, b, 1))
1527 result = mx_stupid_bool_op (Matrix_GE, a, b);
1528 else
1529 return tree_constant ();
1530 break;
1531 case tree::cmp_gt:
1532 result_type = RT_real;
1533 if (m_add_conform (a, b, 1))
1534 result = mx_stupid_bool_op (Matrix_GT, a, b);
1535 else
1536 return tree_constant ();
1537 break;
1538 case tree::cmp_ne:
1539 result_type = RT_real;
1540 if (m_add_conform (a, b, 1))
1541 result = mx_stupid_bool_op (Matrix_NE, a, b);
1542 else
1543 return tree_constant ();
1544 break;
1545 case tree::and:
1546 result_type = RT_real;
1547 if (m_add_conform (a, b, 1))
1548 result = mx_stupid_bool_op (Matrix_AND, a, b);
1549 else
1550 return tree_constant ();
1551 break;
1552 case tree::or:
1553 result_type = RT_real;
1554 if (m_add_conform (a, b, 1))
1555 result = mx_stupid_bool_op (Matrix_OR, a, b);
1556 else
1557 return tree_constant ();
1558 break;
1559 default:
1560 panic_impossible ();
1561 break;
1562 }
1563
1564 assert (result_type != RT_unknown);
1565 if (result_type == RT_real)
1566 return tree_constant (result);
1567 else
1568 return tree_constant (complex_result);
1569 }
1570
1571 /* 9 */
1572 tree_constant
1573 do_binary_op (Complex& a, double b, tree::expression_type t)
1574 {
1575 enum RT { RT_unknown, RT_real, RT_complex };
1576 RT result_type = RT_unknown;
1577
1578 double result = 0.0;
1579 Complex complex_result;
1580 switch (t)
1581 {
1582 case tree::add:
1583 result_type = RT_complex;
1584 complex_result = a + b;
1585 break;
1586 case tree::subtract:
1587 result_type = RT_complex;
1588 complex_result = a - b;
1589 break;
1590 case tree::multiply:
1591 case tree::el_mul:
1592 result_type = RT_complex;
1593 complex_result = a * b;
1594 break;
1595 case tree::divide:
1596 case tree::el_div:
1597 result_type = RT_complex;
1598 if (b == 0.0)
1599 DIVIDE_BY_ZERO_ERROR;
1600 complex_result = a / b;
1601 break;
1602 case tree::leftdiv:
1603 case tree::el_leftdiv:
1604 result_type = RT_complex;
1605 if (a == 0.0)
1606 DIVIDE_BY_ZERO_ERROR;
1607 complex_result = b / a;
1608 break;
1609 case tree::power:
1610 case tree::elem_pow:
1611 return xpow (a, b);
1612 break;
1613 case tree::cmp_lt:
1614 result_type = RT_real;
1615 result = real (a) < b;
1616 break;
1617 case tree::cmp_le:
1618 result_type = RT_real;
1619 result = real (a) <= b;
1620 break;
1621 case tree::cmp_eq:
1622 result_type = RT_real;
1623 result = a == b;
1624 break;
1625 case tree::cmp_ge:
1626 result_type = RT_real;
1627 result = real (a) >= b;
1628 break;
1629 case tree::cmp_gt:
1630 result_type = RT_real;
1631 result = real (a) > b;
1632 break;
1633 case tree::cmp_ne:
1634 result_type = RT_real;
1635 result = a != b;
1636 break;
1637 case tree::and:
1638 result_type = RT_real;
1639 result = ((a != 0.0) && b);
1640 break;
1641 case tree::or:
1642 result_type = RT_real;
1643 result = ((a != 0.0) || b);
1644 break;
1645 default:
1646 panic_impossible ();
1647 break;
1648 }
1649
1650 assert (result_type != RT_unknown);
1651 if (result_type == RT_real)
1652 return tree_constant (result);
1653 else
1654 return tree_constant (complex_result);
1655 }
1656
1657 /* 10 */
1658 tree_constant
1659 do_binary_op (Complex& a, Matrix& b, tree::expression_type t)
1660 {
1661 enum RT { RT_unknown, RT_real, RT_complex };
1662 RT result_type = RT_unknown;
1663
1664 Matrix result;
1665 ComplexMatrix complex_result;
1666 switch (t)
1667 {
1668 case tree::add:
1669 result_type = RT_complex;
1670 complex_result = a + b;
1671 break;
1672 case tree::subtract:
1673 result_type = RT_complex;
1674 complex_result = a - b;
1675 break;
1676 case tree::el_leftdiv:
1677 case tree::leftdiv:
1678 if (a == 0.0)
1679 DIVIDE_BY_ZERO_ERROR;
1680 a = 1.0 / a;
1681 // fall through...
1682 case tree::multiply:
1683 case tree::el_mul:
1684 result_type = RT_complex;
1685 complex_result = a * b;
1686 break;
1687 case tree::el_div:
1688 return x_el_div (a, b);
1689 break;
1690 case tree::divide:
1691 error ("nonconformant right division");
1692 return tree_constant ();
1693 break;
1694 case tree::power:
1695 return xpow (a, b);
1696 break;
1697 case tree::elem_pow:
1698 return elem_xpow (a, b);
1699 break;
1700 case tree::cmp_lt:
1701 result_type = RT_real;
1702 result = mx_stupid_bool_op (Matrix_LT, a, b);
1703 break;
1704 case tree::cmp_le:
1705 result_type = RT_real;
1706 result = mx_stupid_bool_op (Matrix_LE, a, b);
1707 break;
1708 case tree::cmp_eq:
1709 result_type = RT_real;
1710 result = mx_stupid_bool_op (Matrix_EQ, a, b);
1711 break;
1712 case tree::cmp_ge:
1713 result_type = RT_real;
1714 result = mx_stupid_bool_op (Matrix_GE, a, b);
1715 break;
1716 case tree::cmp_gt:
1717 result_type = RT_real;
1718 result = mx_stupid_bool_op (Matrix_GT, a, b);
1719 break;
1720 case tree::cmp_ne:
1721 result_type = RT_real;
1722 result = mx_stupid_bool_op (Matrix_NE, a, b);
1723 break;
1724 case tree::and:
1725 result_type = RT_real;
1726 result = mx_stupid_bool_op (Matrix_AND, a, b);
1727 break;
1728 case tree::or:
1729 result_type = RT_real;
1730 result = mx_stupid_bool_op (Matrix_OR, a, b);
1731 break;
1732 default:
1733 panic_impossible ();
1734 break;
1735 }
1736
1737 assert (result_type != RT_unknown);
1738 if (result_type == RT_real)
1739 return tree_constant (result);
1740 else
1741 return tree_constant (complex_result);
1742 }
1743
1744 /* 11 */
1745 tree_constant
1746 do_binary_op (Complex& a, Complex& b, tree::expression_type t)
1747 {
1748 enum RT { RT_unknown, RT_real, RT_complex };
1749 RT result_type = RT_unknown;
1750
1751 double result = 0.0;
1752 Complex complex_result;
1753 switch (t)
1754 {
1755 case tree::add:
1756 result_type = RT_complex;
1757 complex_result = a + b;
1758 break;
1759 case tree::subtract:
1760 result_type = RT_complex;
1761 complex_result = a - b;
1762 break;
1763 case tree::multiply:
1764 case tree::el_mul:
1765 result_type = RT_complex;
1766 complex_result = a * b;
1767 break;
1768 case tree::divide:
1769 case tree::el_div:
1770 result_type = RT_complex;
1771 if (b == 0.0)
1772 DIVIDE_BY_ZERO_ERROR;
1773 complex_result = a / b;
1774 break;
1775 case tree::leftdiv:
1776 case tree::el_leftdiv:
1777 result_type = RT_complex;
1778 if (a == 0.0)
1779 DIVIDE_BY_ZERO_ERROR;
1780 complex_result = b / a;
1781 break;
1782 case tree::power:
1783 case tree::elem_pow:
1784 return xpow (a, b);
1785 break;
1786 case tree::cmp_lt:
1787 result_type = RT_real;
1788 result = real (a) < real (b);
1789 break;
1790 case tree::cmp_le:
1791 result_type = RT_real;
1792 result = real (a) <= real (b);
1793 break;
1794 case tree::cmp_eq:
1795 result_type = RT_real;
1796 result = a == b;
1797 break;
1798 case tree::cmp_ge:
1799 result_type = RT_real;
1800 result = real (a) >= real (b);
1801 break;
1802 case tree::cmp_gt:
1803 result_type = RT_real;
1804 result = real (a) > real (b);
1805 break;
1806 case tree::cmp_ne:
1807 result_type = RT_real;
1808 result = a != b;
1809 break;
1810 case tree::and:
1811 result_type = RT_real;
1812 result = ((a != 0.0) && (b != 0.0));
1813 break;
1814 case tree::or:
1815 result_type = RT_real;
1816 result = ((a != 0.0) || (b != 0.0));
1817 break;
1818 default:
1819 panic_impossible ();
1820 break;
1821 }
1822
1823 assert (result_type != RT_unknown);
1824 if (result_type == RT_real)
1825 return tree_constant (result);
1826 else
1827 return tree_constant (complex_result);
1828 }
1829
1830 /* 12 */
1831 tree_constant
1832 do_binary_op (Complex& a, ComplexMatrix& b, tree::expression_type t)
1833 {
1834 enum RT { RT_unknown, RT_real, RT_complex };
1835 RT result_type = RT_unknown;
1836
1837 Matrix result;
1838 ComplexMatrix complex_result;
1839 switch (t)
1840 {
1841 case tree::add:
1842 result_type = RT_complex;
1843 complex_result = a + b;
1844 break;
1845 case tree::subtract:
1846 result_type = RT_complex;
1847 complex_result = a - b;
1848 break;
1849 case tree::el_leftdiv:
1850 case tree::leftdiv:
1851 if (a == 0.0)
1852 DIVIDE_BY_ZERO_ERROR;
1853 a = 1.0 / a;
1854 // fall through...
1855 case tree::multiply:
1856 case tree::el_mul:
1857 result_type = RT_complex;
1858 complex_result = a * b;
1859 break;
1860 case tree::el_div:
1861 return x_el_div (a, b);
1862 break;
1863 case tree::divide:
1864 error ("nonconformant right division");
1865 return tree_constant ();
1866 break;
1867 case tree::power:
1868 return xpow (a, b);
1869 break;
1870 case tree::elem_pow:
1871 return elem_xpow (a, b);
1872 break;
1873 case tree::cmp_lt:
1874 result_type = RT_real;
1875 result = mx_stupid_bool_op (Matrix_LT, a, b);
1876 break;
1877 case tree::cmp_le:
1878 result_type = RT_real;
1879 result = mx_stupid_bool_op (Matrix_LE, a, b);
1880 break;
1881 case tree::cmp_eq:
1882 result_type = RT_real;
1883 result = mx_stupid_bool_op (Matrix_EQ, a, b);
1884 break;
1885 case tree::cmp_ge:
1886 result_type = RT_real;
1887 result = mx_stupid_bool_op (Matrix_GE, a, b);
1888 break;
1889 case tree::cmp_gt:
1890 result_type = RT_real;
1891 result = mx_stupid_bool_op (Matrix_GT, a, b);
1892 break;
1893 case tree::cmp_ne:
1894 result_type = RT_real;
1895 result = mx_stupid_bool_op (Matrix_NE, a, b);
1896 break;
1897 case tree::and:
1898 result_type = RT_real;
1899 result = mx_stupid_bool_op (Matrix_AND, a, b);
1900 break;
1901 case tree::or:
1902 result_type = RT_real;
1903 result = mx_stupid_bool_op (Matrix_OR, a, b);
1904 break;
1905 default:
1906 panic_impossible ();
1907 break;
1908 }
1909
1910 assert (result_type != RT_unknown);
1911 if (result_type == RT_real)
1912 return tree_constant (result);
1913 else
1914 return tree_constant (complex_result);
1915 }
1916
1917 /* 13 */
1918 tree_constant
1919 do_binary_op (ComplexMatrix& a, double b, tree::expression_type t)
1920 {
1921 enum RT { RT_unknown, RT_real, RT_complex };
1922 RT result_type = RT_unknown;
1923
1924 Matrix result;
1925 ComplexMatrix complex_result;
1926 switch (t)
1927 {
1928 case tree::add:
1929 result_type = RT_complex;
1930 complex_result = a + b;
1931 break;
1932 case tree::subtract:
1933 result_type = RT_complex;
1934 complex_result = a - b;
1935 break;
1936 case tree::multiply:
1937 case tree::el_mul:
1938 result_type = RT_complex;
1939 complex_result = a * b;
1940 break;
1941 case tree::divide:
1942 case tree::el_div:
1943 result_type = RT_complex;
1944 complex_result = a / b;
1945 break;
1946 case tree::el_leftdiv:
1947 return x_el_div (b, a);
1948 break;
1949 case tree::leftdiv:
1950 error ("nonconformant left division");
1951 return tree_constant ();
1952 break;
1953 case tree::power:
1954 return xpow (a, b);
1955 break;
1956 case tree::elem_pow:
1957 return elem_xpow (a, b);
1958 break;
1959 case tree::cmp_lt:
1960 result_type = RT_real;
1961 result = mx_stupid_bool_op (Matrix_LT, a, b);
1962 break;
1963 case tree::cmp_le:
1964 result_type = RT_real;
1965 result = mx_stupid_bool_op (Matrix_LE, a, b);
1966 break;
1967 case tree::cmp_eq:
1968 result_type = RT_real;
1969 result = mx_stupid_bool_op (Matrix_EQ, a, b);
1970 break;
1971 case tree::cmp_ge:
1972 result_type = RT_real;
1973 result = mx_stupid_bool_op (Matrix_GE, a, b);
1974 break;
1975 case tree::cmp_gt:
1976 result_type = RT_real;
1977 result = mx_stupid_bool_op (Matrix_GT, a, b);
1978 break;
1979 case tree::cmp_ne:
1980 result_type = RT_real;
1981 result = mx_stupid_bool_op (Matrix_NE, a, b);
1982 break;
1983 case tree::and:
1984 result_type = RT_real;
1985 result = mx_stupid_bool_op (Matrix_AND, a, b);
1986 break;
1987 case tree::or:
1988 result_type = RT_real;
1989 result = mx_stupid_bool_op (Matrix_OR, a, b);
1990 break;
1991 default:
1992 panic_impossible ();
1993 break;
1994 }
1995
1996 assert (result_type != RT_unknown);
1997 if (result_type == RT_real)
1998 return tree_constant (result);
1999 else
2000 return tree_constant (complex_result);
2001 }
2002
2003 /* 14 */
2004 tree_constant
2005 do_binary_op (ComplexMatrix& a, Matrix& b, tree::expression_type t)
2006 {
2007 enum RT { RT_unknown, RT_real, RT_complex };
2008 RT result_type = RT_unknown;
2009
2010 Matrix result;
2011 ComplexMatrix complex_result;
2012 switch (t)
2013 {
2014 case tree::add:
2015 result_type = RT_complex;
2016 if (m_add_conform (a, b, 1))
2017 complex_result = a + b;
2018 else
2019 return tree_constant ();
2020 break;
2021 case tree::subtract:
2022 result_type = RT_complex;
2023 if (m_add_conform (a, b, 1))
2024 complex_result = a - b;
2025 else
2026 return tree_constant ();
2027 break;
2028 case tree::el_mul:
2029 result_type = RT_complex;
2030 if (m_add_conform (a, b, 1))
2031 complex_result = a.product (b);
2032 else
2033 return tree_constant ();
2034 break;
2035 case tree::multiply:
2036 result_type = RT_complex;
2037 if (m_mul_conform (a, b, 1))
2038 complex_result = a * b;
2039 else
2040 return tree_constant ();
2041 break;
2042 case tree::el_div:
2043 result_type = RT_complex;
2044 if (m_add_conform (a, b, 1))
2045 complex_result = a.quotient (b);
2046 else
2047 return tree_constant ();
2048 break;
2049 case tree::el_leftdiv:
2050 result_type = RT_complex;
2051 if (m_add_conform (a, b, 1))
2052 complex_result = a.quotient (b);
2053 else
2054 return tree_constant ();
2055 break;
2056 case tree::leftdiv:
2057 return xleftdiv (a, b);
2058 break;
2059 case tree::divide:
2060 return xdiv (a, b);
2061 break;
2062 case tree::power:
2063 error ("can't do A ^ B for A and B both matrices");
2064 return tree_constant ();
2065 break;
2066 case tree::elem_pow:
2067 if (m_add_conform (a, b, 1))
2068 return elem_xpow (a, b);
2069 else
2070 return tree_constant ();
2071 break;
2072 case tree::cmp_lt:
2073 result_type = RT_real;
2074 if (m_add_conform (a, b, 1))
2075 result = mx_stupid_bool_op (Matrix_LT, a, b);
2076 else
2077 return tree_constant ();
2078 break;
2079 case tree::cmp_le:
2080 result_type = RT_real;
2081 if (m_add_conform (a, b, 1))
2082 result = mx_stupid_bool_op (Matrix_LE, a, b);
2083 else
2084 return tree_constant ();
2085 break;
2086 case tree::cmp_eq:
2087 result_type = RT_real;
2088 if (m_add_conform (a, b, 1))
2089 result = mx_stupid_bool_op (Matrix_EQ, a, b);
2090 else
2091 return tree_constant ();
2092 break;
2093 case tree::cmp_ge:
2094 result_type = RT_real;
2095 if (m_add_conform (a, b, 1))
2096 result = mx_stupid_bool_op (Matrix_GE, a, b);
2097 else
2098 return tree_constant ();
2099 break;
2100 case tree::cmp_gt:
2101 result_type = RT_real;
2102 if (m_add_conform (a, b, 1))
2103 result = mx_stupid_bool_op (Matrix_GT, a, b);
2104 else
2105 return tree_constant ();
2106 break;
2107 case tree::cmp_ne:
2108 result_type = RT_real;
2109 if (m_add_conform (a, b, 1))
2110 result = mx_stupid_bool_op (Matrix_NE, a, b);
2111 else
2112 return tree_constant ();
2113 break;
2114 case tree::and:
2115 result_type = RT_real;
2116 if (m_add_conform (a, b, 1))
2117 result = mx_stupid_bool_op (Matrix_AND, a, b);
2118 else
2119 return tree_constant ();
2120 break;
2121 case tree::or:
2122 result_type = RT_real;
2123 if (m_add_conform (a, b, 1))
2124 result = mx_stupid_bool_op (Matrix_OR, a, b);
2125 else
2126 return tree_constant ();
2127 break;
2128 default:
2129 panic_impossible ();
2130 break;
2131 }
2132
2133 assert (result_type != RT_unknown);
2134 if (result_type == RT_real)
2135 return tree_constant (result);
2136 else
2137 return tree_constant (complex_result);
2138 }
2139
2140 /* 15 */
2141 tree_constant
2142 do_binary_op (ComplexMatrix& a, Complex& b, tree::expression_type t)
2143 {
2144 enum RT { RT_unknown, RT_real, RT_complex };
2145 RT result_type = RT_unknown;
2146
2147 Matrix result;
2148 ComplexMatrix complex_result;
2149 switch (t)
2150 {
2151 case tree::add:
2152 result_type = RT_complex;
2153 complex_result = a + b;
2154 break;
2155 case tree::subtract:
2156 result_type = RT_complex;
2157 complex_result = a - b;
2158 break;
2159 case tree::multiply:
2160 case tree::el_mul:
2161 result_type = RT_complex;
2162 complex_result = a * b;
2163 break;
2164 case tree::divide:
2165 case tree::el_div:
2166 result_type = RT_complex;
2167 complex_result = a / b;
2168 break;
2169 case tree::el_leftdiv:
2170 return x_el_div (b, a);
2171 break;
2172 case tree::leftdiv:
2173 error ("nonconformant left division");
2174 return tree_constant ();
2175 break;
2176 case tree::power:
2177 return xpow (a, b);
2178 break;
2179 case tree::elem_pow:
2180 return elem_xpow (a, b);
2181 break;
2182 case tree::cmp_lt:
2183 result_type = RT_real;
2184 result = mx_stupid_bool_op (Matrix_LT, a, b);
2185 break;
2186 case tree::cmp_le:
2187 result_type = RT_real;
2188 result = mx_stupid_bool_op (Matrix_LE, a, b);
2189 break;
2190 case tree::cmp_eq:
2191 result_type = RT_real;
2192 result = mx_stupid_bool_op (Matrix_EQ, a, b);
2193 break;
2194 case tree::cmp_ge:
2195 result_type = RT_real;
2196 result = mx_stupid_bool_op (Matrix_GE, a, b);
2197 break;
2198 case tree::cmp_gt:
2199 result_type = RT_real;
2200 result = mx_stupid_bool_op (Matrix_GT, a, b);
2201 break;
2202 case tree::cmp_ne:
2203 result_type = RT_real;
2204 result = mx_stupid_bool_op (Matrix_NE, a, b);
2205 break;
2206 case tree::and:
2207 result_type = RT_real;
2208 result = mx_stupid_bool_op (Matrix_AND, a, b);
2209 break;
2210 case tree::or:
2211 result_type = RT_real;
2212 result = mx_stupid_bool_op (Matrix_OR, a, b);
2213 break;
2214 default:
2215 panic_impossible ();
2216 break;
2217 }
2218
2219 assert (result_type != RT_unknown);
2220 if (result_type == RT_real)
2221 return tree_constant (result);
2222 else
2223 return tree_constant (complex_result);
2224 }
2225
2226 /* 16 */
2227 tree_constant
2228 do_binary_op (ComplexMatrix& a, ComplexMatrix& b, tree::expression_type t)
2229 {
2230 enum RT { RT_unknown, RT_real, RT_complex };
2231 RT result_type = RT_unknown;
2232
2233 Matrix result;
2234 ComplexMatrix complex_result;
2235 switch (t)
2236 {
2237 case tree::add:
2238 result_type = RT_complex;
2239 if (m_add_conform (a, b, 1))
2240 complex_result = a + b;
2241 else
2242 return tree_constant ();
2243 break;
2244 case tree::subtract:
2245 result_type = RT_complex;
2246 if (m_add_conform (a, b, 1))
2247 complex_result = a - b;
2248 else
2249 return tree_constant ();
2250 break;
2251 case tree::el_mul:
2252 result_type = RT_complex;
2253 if (m_add_conform (a, b, 1))
2254 complex_result = a.product (b);
2255 else
2256 return tree_constant ();
2257 break;
2258 case tree::multiply:
2259 result_type = RT_complex;
2260 if (m_mul_conform (a, b, 1))
2261 complex_result = a * b;
2262 else
2263 return tree_constant ();
2264 break;
2265 case tree::el_div:
2266 result_type = RT_complex;
2267 if (m_add_conform (a, b, 1))
2268 complex_result = a.quotient (b);
2269 else
2270 return tree_constant ();
2271 break;
2272 case tree::el_leftdiv:
2273 result_type = RT_complex;
2274 if (m_add_conform (a, b, 1))
2275 complex_result = b.quotient (a);
2276 else
2277 return tree_constant ();
2278 break;
2279 case tree::leftdiv:
2280 return xleftdiv (a, b);
2281 break;
2282 case tree::divide:
2283 return xdiv (a, b);
2284 break;
2285 case tree::power:
2286 error ("can't do A ^ B for A and B both matrices");
2287 return tree_constant ();
2288 break;
2289 case tree::elem_pow:
2290 if (m_add_conform (a, b, 1))
2291 return elem_xpow (a, b);
2292 else
2293 return tree_constant ();
2294 break;
2295 case tree::cmp_lt:
2296 result_type = RT_real;
2297 if (m_add_conform (a, b, 1))
2298 result = mx_stupid_bool_op (Matrix_LT, a, b);
2299 else
2300 return tree_constant ();
2301 break;
2302 case tree::cmp_le:
2303 result_type = RT_real;
2304 if (m_add_conform (a, b, 1))
2305 result = mx_stupid_bool_op (Matrix_LE, a, b);
2306 else
2307 return tree_constant ();
2308 break;
2309 case tree::cmp_eq:
2310 result_type = RT_real;
2311 if (m_add_conform (a, b, 1))
2312 result = mx_stupid_bool_op (Matrix_EQ, a, b);
2313 else
2314 return tree_constant ();
2315 break;
2316 case tree::cmp_ge:
2317 result_type = RT_real;
2318 if (m_add_conform (a, b, 1))
2319 result = mx_stupid_bool_op (Matrix_GE, a, b);
2320 else
2321 return tree_constant ();
2322 break;
2323 case tree::cmp_gt:
2324 result_type = RT_real;
2325 if (m_add_conform (a, b, 1))
2326 result = mx_stupid_bool_op (Matrix_GT, a, b);
2327 else
2328 return tree_constant ();
2329 break;
2330 case tree::cmp_ne:
2331 result_type = RT_real;
2332 if (m_add_conform (a, b, 1))
2333 result = mx_stupid_bool_op (Matrix_NE, a, b);
2334 else
2335 return tree_constant ();
2336 break;
2337 case tree::and:
2338 result_type = RT_real;
2339 if (m_add_conform (a, b, 1))
2340 result = mx_stupid_bool_op (Matrix_AND, a, b);
2341 else
2342 return tree_constant ();
2343 break;
2344 case tree::or:
2345 result_type = RT_real;
2346 if (m_add_conform (a, b, 1))
2347 result = mx_stupid_bool_op (Matrix_OR, a, b);
2348 else
2349 return tree_constant ();
2350 break;
2351 default:
2352 panic_impossible ();
2353 break;
2354 }
2355
2356 assert (result_type != RT_unknown);
2357 if (result_type == RT_real)
2358 return tree_constant (result);
2359 else
2360 return tree_constant (complex_result);
2361 }
2362
2363 /*
2364 ;;; Local Variables: ***
2365 ;;; mode: C++ ***
2366 ;;; page-delimiter: "^/\\*" ***
2367 ;;; End: ***
2368 */