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