Mercurial > octave
changeset 493:1391e7ed65f6
[project @ 1994-07-06 14:54:58 by jwe]
| author | jwe |
|---|---|
| date | Wed, 06 Jul 1994 14:54:58 +0000 |
| parents | d45bdf960233 |
| children | 2c4d694b87e9 |
| files | src/pt-const.cc |
| diffstat | 1 files changed, 1092 insertions(+), 5811 deletions(-) [+] |
line wrap: on
line diff
--- a/src/pt-const.cc Wed Jul 06 14:54:58 1994 +0000 +++ b/src/pt-const.cc Wed Jul 06 14:54:58 1994 +0000 @@ -29,2482 +29,35 @@ #pragma implementation #endif -#include <ctype.h> -#include <string.h> -#include <fstream.h> +#include <strstream.h> #include <iostream.h> -#include <strstream.h> -#include <math.h> +#include <fstream.h> +#include <string.h> +#include <ctype.h> + +#include "EIG.h" -#include "mx-base.h" -#include "Range.h" - +#include "unwind-prot.h" +#include "tree-const.h" +#include "user-prefs.h" #include "variables.h" +#include "octave.h" #include "error.h" #include "gripes.h" -#include "user-prefs.h" -#include "utils.h" -#include "pager.h" -#include "mappers.h" -#include "pr-output.h" -#include "tree-const.h" -#include "arith-ops.h" -#include "idx-vector.h" -#include "unwind-prot.h" -#include "octave.h" #include "input.h" #include "octave-hist.h" +#include "pager.h" +#include "utils.h" #include "parse.h" #include "lex.h" -#include "tc-inlines.cc" - -// A couple of handy helper functions. - -static int -any_element_less_than (const Matrix& a, double val) -{ - int nr = a.rows (); - int nc = a.columns (); - for (int j = 0; j < nc; j++) - for (int i = 0; i < nr; i++) - if (a.elem (i, j) < val) - return 1; - return 0; -} - -static int -any_element_greater_than (const Matrix& a, double val) -{ - int nr = a.rows (); - int nc = a.columns (); - for (int j = 0; j < nc; j++) - for (int i = 0; i < nr; i++) - if (a.elem (i, j) > val) - return 1; - return 0; -} - -static int -any_element_is_complex (const ComplexMatrix& a) -{ - int nr = a.rows (); - int nc = a.columns (); - for (int j = 0; j < nc; j++) - for (int i = 0; i < nr; i++) - if (imag (a.elem (i, j)) != 0.0) - return 1; - return 0; -} - -// Now, the classes. - -/* - * The real representation of constants. - */ -tree_constant_rep::tree_constant_rep (void) -{ - type_tag = unknown_constant; -} - -tree_constant_rep::tree_constant_rep (double d) -{ - scalar = d; - type_tag = scalar_constant; -} - -tree_constant_rep::tree_constant_rep (const Matrix& m) -{ - if (m.rows () == 1 && m.columns () == 1) - { - scalar = m.elem (0, 0); - type_tag = scalar_constant; - } - else - { - matrix = new Matrix (m); - type_tag = matrix_constant; - } -} - -tree_constant_rep::tree_constant_rep (const DiagMatrix& d) -{ - if (d.rows () == 1 && d.columns () == 1) - { - scalar = d.elem (0, 0); - type_tag = scalar_constant; - } - else - { - matrix = new Matrix (d); - type_tag = matrix_constant; - } -} - -tree_constant_rep::tree_constant_rep (const RowVector& v, int - prefer_column_vector) -{ - int len = v.capacity (); - if (len == 1) - { - scalar = v.elem (0); - type_tag = scalar_constant; - } - else - { - int pcv = (prefer_column_vector < 0) - ? user_pref.prefer_column_vectors - : prefer_column_vector; - - if (pcv) - { - Matrix m (len, 1); - for (int i = 0; i < len; i++) - m.elem (i, 0) = v.elem (i); - matrix = new Matrix (m); - type_tag = matrix_constant; - } - else - { - Matrix m (1, len); - for (int i = 0; i < len; i++) - m.elem (0, i) = v.elem (i); - matrix = new Matrix (m); - type_tag = matrix_constant; - } - } -} - -tree_constant_rep::tree_constant_rep (const ColumnVector& v, - int prefer_column_vector) -{ - int len = v.capacity (); - if (len == 1) - { - scalar = v.elem (0); - type_tag = scalar_constant; - } - else - { - int pcv = (prefer_column_vector < 0) - ? user_pref.prefer_column_vectors - : prefer_column_vector; - - if (pcv) - { - Matrix m (len, 1); - for (int i = 0; i < len; i++) - m.elem (i, 0) = v.elem (i); - matrix = new Matrix (m); - type_tag = matrix_constant; - } - else - { - Matrix m (1, len); - for (int i = 0; i < len; i++) - m.elem (0, i) = v.elem (i); - matrix = new Matrix (m); - type_tag = matrix_constant; - } - } -} - -tree_constant_rep::tree_constant_rep (const Complex& c) -{ - complex_scalar = new Complex (c); - type_tag = complex_scalar_constant; -} - -tree_constant_rep::tree_constant_rep (const ComplexMatrix& m) -{ - if (m.rows () == 1 && m.columns () == 1) - { - complex_scalar = new Complex (m.elem (0, 0)); - type_tag = complex_scalar_constant; - } - else - { - complex_matrix = new ComplexMatrix (m); - type_tag = complex_matrix_constant; - } -} - -tree_constant_rep::tree_constant_rep (const ComplexDiagMatrix& d) -{ - if (d.rows () == 1 && d.columns () == 1) - { - complex_scalar = new Complex (d.elem (0, 0)); - type_tag = complex_scalar_constant; - } - else - { - complex_matrix = new ComplexMatrix (d); - type_tag = complex_matrix_constant; - } -} - -tree_constant_rep::tree_constant_rep (const ComplexRowVector& v, - int prefer_column_vector) -{ - int len = v.capacity (); - if (len == 1) - { - complex_scalar = new Complex (v.elem (0)); - type_tag = complex_scalar_constant; - } - else - { - int pcv = (prefer_column_vector < 0) - ? user_pref.prefer_column_vectors - : prefer_column_vector; - - if (pcv) - { - ComplexMatrix m (len, 1); - for (int i = 0; i < len; i++) - m.elem (i, 0) = v.elem (i); - complex_matrix = new ComplexMatrix (m); - type_tag = complex_matrix_constant; - } - else - { - ComplexMatrix m (1, len); - for (int i = 0; i < len; i++) - m.elem (0, i) = v.elem (i); - complex_matrix = new ComplexMatrix (m); - type_tag = complex_matrix_constant; - } - } -} - -tree_constant_rep::tree_constant_rep (const ComplexColumnVector& v, - int prefer_column_vector) -{ - int len = v.capacity (); - if (len == 1) - { - complex_scalar = new Complex (v.elem (0)); - type_tag = complex_scalar_constant; - } - else - { - int pcv = (prefer_column_vector < 0) - ? user_pref.prefer_column_vectors - : prefer_column_vector; - - if (pcv) - { - ComplexMatrix m (len, 1); - for (int i = 0; i < len; i++) - m.elem (i, 0) = v.elem (i); - complex_matrix = new ComplexMatrix (m); - type_tag = complex_matrix_constant; - } - else - { - ComplexMatrix m (1, len); - for (int i = 0; i < len; i++) - m.elem (0, i) = v.elem (i); - complex_matrix = new ComplexMatrix (m); - type_tag = complex_matrix_constant; - } - } -} - -tree_constant_rep::tree_constant_rep (const char *s) -{ - string = strsave (s); - type_tag = string_constant; -} - -tree_constant_rep::tree_constant_rep (double b, double l, double i) -{ - range = new Range (b, l, i); - int nel = range->nelem (); - if (nel < 0) - { - delete range; - type_tag = unknown_constant; - if (nel == -1) - ::error ("number of elements in range exceeds INT_MAX"); - else - ::error ("invalid range"); - } - else if (nel > 1) - type_tag = range_constant; - else - { - delete range; - if (nel == 1) - { - scalar = b; - type_tag = scalar_constant; - } - else if (nel == 0) - { - matrix = new Matrix (); - type_tag = matrix_constant; - } - else - panic_impossible (); - } -} - -tree_constant_rep::tree_constant_rep (const Range& r) -{ - if (r.nelem () > 1) - { - range = new Range (r); - type_tag = range_constant; - } - else if (r.nelem () == 1) - { - scalar = r.base (); - type_tag = scalar_constant; - } - else if (r.nelem () == 0) - { - matrix = new Matrix (); - type_tag = matrix_constant; - } - else - panic_impossible (); -} - -tree_constant_rep::tree_constant_rep (tree_constant_rep::constant_type t) -{ - assert (t == magic_colon); - - type_tag = magic_colon; -} - -tree_constant_rep::tree_constant_rep (const tree_constant_rep& t) -{ - type_tag = t.type_tag; - - switch (t.type_tag) - { - case unknown_constant: - break; - case scalar_constant: - scalar = t.scalar; - break; - case matrix_constant: - matrix = new Matrix (*(t.matrix)); - break; - case string_constant: - string = strsave (t.string); - break; - case complex_matrix_constant: - complex_matrix = new ComplexMatrix (*(t.complex_matrix)); - break; - case complex_scalar_constant: - complex_scalar = new Complex (*(t.complex_scalar)); - break; - case range_constant: - range = new Range (*(t.range)); - break; - case magic_colon: - break; - default: - panic_impossible (); - break; - } -} - -tree_constant_rep::~tree_constant_rep (void) -{ - switch (type_tag) - { - case unknown_constant: - break; - case scalar_constant: - break; - case matrix_constant: - delete matrix; - break; - case complex_scalar_constant: - delete complex_scalar; - break; - case complex_matrix_constant: - delete complex_matrix; - break; - case string_constant: - delete [] string; - break; - case range_constant: - delete range; - break; - case magic_colon: - break; - default: - panic_impossible (); - break; - } -} - -#if defined (MDEBUG) -void * -tree_constant_rep::operator new (size_t size) -{ - tree_constant_rep *p = ::new tree_constant_rep; - cerr << "tree_constant_rep::new(): " << p << "\n"; - return p; -} - -void -tree_constant_rep::operator delete (void *p, size_t size) -{ - cerr << "tree_constant_rep::delete(): " << p << "\n"; - ::delete p; -} +#ifndef MAX +#define MAX(a,b) ((a) > (b) ? (a) : (b)) #endif -void -tree_constant_rep::resize (int i, int j) -{ - switch (type_tag) - { - case matrix_constant: - matrix->resize (i, j); - break; - case complex_matrix_constant: - complex_matrix->resize (i, j); - break; - default: - panic_impossible (); - break; - } -} - -void -tree_constant_rep::resize (int i, int j, double val) -{ - switch (type_tag) - { - case matrix_constant: - matrix->resize (i, j, val); - break; - case complex_matrix_constant: - complex_matrix->resize (i, j, val); - break; - default: - panic_impossible (); - break; - } -} - -void -tree_constant_rep::maybe_resize (int i, int j) -{ - int nr = rows (); - int nc = columns (); - - i++; - j++; - - assert (i > 0 && j > 0); - - if (i > nr || j > nc) - { - if (user_pref.resize_on_range_error) - resize (MAX (i, nr), MAX (j, nc), 0.0); - else - { - if (i > nr) - ::error ("row index = %d exceeds max row dimension = %d", i, nr); - - if (j > nc) - ::error ("column index = %d exceeds max column dimension = %d", - j, nc); - } - } -} - -void -tree_constant_rep::maybe_resize (int i, force_orient f_orient = no_orient) -{ - int nr = rows (); - int nc = columns (); - - i++; - - assert (i >= 0 && (nr <= 1 || nc <= 1)); - -// This function never reduces the size of a vector, and all vectors -// have dimensions of at least 0x0. If i is 0, it is either because -// a vector has been indexed with a vector of all zeros (in which case -// the index vector is empty and nothing will happen) or a vector has -// been indexed with 0 (an error which will be caught elsewhere). - if (i == 0) - return; - - if (nr <= 1 && nc <= 1 && i >= 1) - { - if (user_pref.resize_on_range_error) - { - if (f_orient == row_orient) - resize (1, i, 0.0); - else if (f_orient == column_orient) - resize (i, 1, 0.0); - else if (user_pref.prefer_column_vectors) - resize (i, 1, 0.0); - else - resize (1, i, 0.0); - } - else - ::error ("matrix index = %d exceeds max dimension = %d", i, nc); - } - else if (nr == 1 && i > nc) - { - if (user_pref.resize_on_range_error) - resize (1, i, 0.0); - else - ::error ("matrix index = %d exceeds max dimension = %d", i, nc); - } - else if (nc == 1 && i > nr) - { - if (user_pref.resize_on_range_error) - resize (i, 1, 0.0); - else - ::error ("matrix index = %d exceeds max dimension = ", i, nc); - } -} - -double -tree_constant_rep::to_scalar (void) const -{ - tree_constant tmp = make_numeric (); - - double retval = 0.0; - - switch (tmp.const_type ()) - { - case tree_constant_rep::scalar_constant: - case tree_constant_rep::complex_scalar_constant: - retval = tmp.double_value (); - break; - case tree_constant_rep::matrix_constant: - if (user_pref.do_fortran_indexing) - { - Matrix m = tmp.matrix_value (); - retval = m (0, 0); - } - break; - case tree_constant_rep::complex_matrix_constant: - if (user_pref.do_fortran_indexing) - { - int flag = user_pref.ok_to_lose_imaginary_part; - if (flag == -1) - warning ("implicit conversion of complex value to real value"); - - if (flag != 0) - { - ComplexMatrix m = tmp.complex_matrix_value (); - return ::real (m (0, 0)); - } - else - jump_to_top_level (); - } - else - { - ::error ("complex matrix used in invalid context"); - jump_to_top_level (); - } - break; - default: - break; - } - return retval; -} - -ColumnVector -tree_constant_rep::to_vector (void) const -{ - tree_constant tmp = make_numeric (); - - ColumnVector retval; - - switch (tmp.const_type ()) - { - case tree_constant_rep::scalar_constant: - case tree_constant_rep::complex_scalar_constant: - retval.resize (1); - retval.elem (0) = tmp.double_value (); - break; - case tree_constant_rep::complex_matrix_constant: - case tree_constant_rep::matrix_constant: - { - Matrix m = tmp.matrix_value (); - int nr = m.rows (); - int nc = m.columns (); - if (nr == 1) - { - retval.resize (nc); - for (int i = 0; i < nc; i++) - retval.elem (i) = m (0, i); - } - else if (nc == 1) - { - retval.resize (nr); - for (int i = 0; i < nr; i++) - retval.elem (i) = m.elem (i, 0); - } - } - break; - default: - panic_impossible (); - break; - } - return retval; -} - -Matrix -tree_constant_rep::to_matrix (void) const -{ - tree_constant tmp = make_numeric (); - - Matrix retval; - - switch (tmp.const_type ()) - { - case tree_constant_rep::scalar_constant: - retval.resize (1, 1); - retval.elem (0, 0) = tmp.double_value (); - break; - case tree_constant_rep::matrix_constant: - retval = tmp.matrix_value (); - break; - default: - break; - } - return retval; -} - -tree_constant_rep::constant_type -tree_constant_rep::force_numeric (int force_str_conv = 0) -{ - switch (type_tag) - { - case scalar_constant: - case matrix_constant: - case complex_scalar_constant: - case complex_matrix_constant: - break; - case string_constant: - { - if (! force_str_conv && ! user_pref.implicit_str_to_num_ok) - { - ::error ("failed to convert `%s' to a numeric type --", string); - ::error ("default conversion turned off"); -// Abort! - jump_to_top_level (); - } - - int len = strlen (string); - if (len > 1) - { - type_tag = matrix_constant; - Matrix *tm = new Matrix (1, len); - for (int i = 0; i < len; i++) - tm->elem (0, i) = toascii ((int) string[i]); - matrix = tm; - } - else if (len == 1) - { - type_tag = scalar_constant; - scalar = toascii ((int) string[0]); - } - else if (len == 0) - { - type_tag = matrix_constant; - matrix = new Matrix (0, 0); - } - else - panic_impossible (); - } - break; - case range_constant: - { - int len = range->nelem (); - if (len > 1) - { - type_tag = matrix_constant; - Matrix *tm = new Matrix (1, len); - double b = range->base (); - double increment = range->inc (); - for (int i = 0; i < len; i++) - tm->elem (0, i) = b + i * increment; - matrix = tm; - } - else if (len == 1) - { - type_tag = scalar_constant; - scalar = range->base (); - } - } - break; - case magic_colon: - default: - panic_impossible (); - break; - } - return type_tag; -} - -tree_constant -tree_constant_rep::make_numeric (int force_str_conv = 0) const -{ - tree_constant retval; - switch (type_tag) - { - case scalar_constant: - retval = tree_constant (scalar); - break; - case matrix_constant: - retval = tree_constant (*matrix); - break; - case complex_scalar_constant: - retval = tree_constant (*complex_scalar); - break; - case complex_matrix_constant: - retval = tree_constant (*complex_matrix); - break; - case string_constant: - retval = tree_constant (string); - retval.force_numeric (force_str_conv); - break; - case range_constant: - retval = tree_constant (*range); - retval.force_numeric (force_str_conv); - break; - case magic_colon: - default: - panic_impossible (); - break; - } - return retval; -} - -tree_constant -do_binary_op (tree_constant& a, tree_constant& b, tree::expression_type t) -{ - tree_constant ans; - - int first_empty = (a.rows () == 0 || a.columns () == 0); - int second_empty = (b.rows () == 0 || b.columns () == 0); - - if (first_empty || second_empty) - { - int flag = user_pref.propagate_empty_matrices; - if (flag < 0) - warning ("binary operation on empty matrix"); - else if (flag == 0) - { - ::error ("invalid binary operation on empty matrix"); - return ans; - } - } - - tree_constant tmp_a = a.make_numeric (); - tree_constant tmp_b = b.make_numeric (); - - tree_constant_rep::constant_type a_type = tmp_a.const_type (); - tree_constant_rep::constant_type b_type = tmp_b.const_type (); - - double d1, d2; - Matrix m1, m2; - Complex c1, c2; - ComplexMatrix cm1, cm2; - - switch (a_type) - { - case tree_constant_rep::scalar_constant: - d1 = tmp_a.double_value (); - switch (b_type) - { - case tree_constant_rep::scalar_constant: - d2 = tmp_b.double_value (); - ans = do_binary_op (d1, d2, t); - break; - case tree_constant_rep::matrix_constant: - m2 = tmp_b.matrix_value (); - ans = do_binary_op (d1, m2, t); - break; - case tree_constant_rep::complex_scalar_constant: - c2 = tmp_b.complex_value (); - ans = do_binary_op (d1, c2, t); - break; - case tree_constant_rep::complex_matrix_constant: - cm2 = tmp_b.complex_matrix_value (); - ans = do_binary_op (d1, cm2, t); - break; - case tree_constant_rep::magic_colon: - default: - panic_impossible (); - break; - } - break; - case tree_constant_rep::matrix_constant: - m1 = tmp_a.matrix_value (); - switch (b_type) - { - case tree_constant_rep::scalar_constant: - d2 = tmp_b.double_value (); - ans = do_binary_op (m1, d2, t); - break; - case tree_constant_rep::matrix_constant: - m2 = tmp_b.matrix_value (); - ans = do_binary_op (m1, m2, t); - break; - case tree_constant_rep::complex_scalar_constant: - c2 = tmp_b.complex_value (); - ans = do_binary_op (m1, c2, t); - break; - case tree_constant_rep::complex_matrix_constant: - cm2 = tmp_b.complex_matrix_value (); - ans = do_binary_op (m1, cm2, t); - break; - case tree_constant_rep::magic_colon: - default: - panic_impossible (); - break; - } - break; - case tree_constant_rep::complex_scalar_constant: - c1 = tmp_a.complex_value (); - switch (b_type) - { - case tree_constant_rep::scalar_constant: - d2 = tmp_b.double_value (); - ans = do_binary_op (c1, d2, t); - break; - case tree_constant_rep::matrix_constant: - m2 = tmp_b.matrix_value (); - ans = do_binary_op (c1, m2, t); - break; - case tree_constant_rep::complex_scalar_constant: - c2 = tmp_b.complex_value (); - ans = do_binary_op (c1, c2, t); - break; - case tree_constant_rep::complex_matrix_constant: - cm2 = tmp_b.complex_matrix_value (); - ans = do_binary_op (c1, cm2, t); - break; - case tree_constant_rep::magic_colon: - default: - panic_impossible (); - break; - } - break; - case tree_constant_rep::complex_matrix_constant: - cm1 = tmp_a.complex_matrix_value (); - switch (b_type) - { - case tree_constant_rep::scalar_constant: - d2 = tmp_b.double_value (); - ans = do_binary_op (cm1, d2, t); - break; - case tree_constant_rep::matrix_constant: - m2 = tmp_b.matrix_value (); - ans = do_binary_op (cm1, m2, t); - break; - case tree_constant_rep::complex_scalar_constant: - c2 = tmp_b.complex_value (); - ans = do_binary_op (cm1, c2, t); - break; - case tree_constant_rep::complex_matrix_constant: - cm2 = tmp_b.complex_matrix_value (); - ans = do_binary_op (cm1, cm2, t); - break; - case tree_constant_rep::magic_colon: - default: - panic_impossible (); - break; - } - break; - case tree_constant_rep::magic_colon: - default: - panic_impossible (); - break; - } - return ans; -} - -tree_constant -do_unary_op (tree_constant& a, tree::expression_type t) -{ - tree_constant ans; - - if (a.rows () == 0 || a.columns () == 0) - { - int flag = user_pref.propagate_empty_matrices; - if (flag < 0) - warning ("unary operation on empty matrix"); - else if (flag == 0) - { - ::error ("invalid unary operation on empty matrix"); - return ans; - } - } - - tree_constant tmp_a = a.make_numeric (); - - switch (tmp_a.const_type ()) - { - case tree_constant_rep::scalar_constant: - ans = do_unary_op (tmp_a.double_value (), t); - break; - case tree_constant_rep::matrix_constant: - { - Matrix m = tmp_a.matrix_value (); - ans = do_unary_op (m, t); - } - break; - case tree_constant_rep::complex_scalar_constant: - ans = do_unary_op (tmp_a.complex_value (), t); - break; - case tree_constant_rep::complex_matrix_constant: - { - ComplexMatrix m = tmp_a.complex_matrix_value (); - ans = do_unary_op (m, t); - } - break; - case tree_constant_rep::magic_colon: - default: - panic_impossible (); - break; - } - return ans; -} - -void -tree_constant_rep::bump_value (tree::expression_type etype) -{ - switch (etype) - { - case tree::increment: - switch (type_tag) - { - case scalar_constant: - scalar++; - break; - case matrix_constant: - *matrix = *matrix + 1.0; - break; - case complex_scalar_constant: - *complex_scalar = *complex_scalar + 1.0; - break; - case complex_matrix_constant: - *complex_matrix = *complex_matrix + 1.0; - break; - case string_constant: - ::error ("string++ and ++string not implemented yet, ok?"); - break; - case range_constant: - range->set_base (range->base () + 1.0); - range->set_limit (range->limit () + 1.0); - break; - case magic_colon: - default: - panic_impossible (); - break; - } - break; - case tree::decrement: - switch (type_tag) - { - case scalar_constant: - scalar--; - break; - case matrix_constant: - *matrix = *matrix - 1.0; - break; - case string_constant: - ::error ("string-- and -- string not implemented yet, ok?"); - break; - case range_constant: - range->set_base (range->base () - 1.0); - range->set_limit (range->limit () - 1.0); - break; - case magic_colon: - default: - panic_impossible (); - break; - } - break; - default: - panic_impossible (); - break; - } -} - -void -tree_constant_rep::eval (int print) -{ - if (error_state) - return; - - switch (type_tag) - { - case complex_scalar_constant: - if (::imag (*complex_scalar) == 0.0) - { - double d = ::real (*complex_scalar); - delete complex_scalar; - scalar = d; - type_tag = scalar_constant; - } - break; - case complex_matrix_constant: - if (! any_element_is_complex (*complex_matrix)) - { - Matrix *m = new Matrix (::real (*complex_matrix)); - delete complex_matrix; - matrix = m; - type_tag = matrix_constant; - } - break; - case scalar_constant: - case matrix_constant: - case string_constant: - case range_constant: - case magic_colon: - break; - default: - panic_impossible (); - break; - } - -// Avoid calling rows() and columns() for things like magic_colon. - - int nr = 1; - int nc = 1; - if (type_tag == matrix_constant - || type_tag == complex_matrix_constant - || type_tag == range_constant) - { - nr = rows (); - nc = columns (); - } - - switch (type_tag) - { - case matrix_constant: - if (nr == 1 && nc == 1) - { - double d = matrix->elem (0, 0); - delete matrix; - scalar = d; - type_tag = scalar_constant; - } - break; - case complex_matrix_constant: - if (nr == 1 && nc == 1) - { - Complex c = complex_matrix->elem (0, 0); - delete complex_matrix; - complex_scalar = new Complex (c); - type_tag = complex_scalar_constant; - } - break; - case range_constant: - if (nr == 1 && nc == 1) - { - double d = range->base (); - delete range; - scalar = d; - type_tag = scalar_constant; - } - break; - default: - break; - } - - if (print) - { - ostrstream output_buf; - switch (type_tag) - { - case scalar_constant: - octave_print_internal (output_buf, scalar); - break; - case matrix_constant: - if (nr == 0 || nc == 0) - { - output_buf << "[]"; - if (user_pref.print_empty_dimensions) - output_buf << "(" << nr << "x" << nc << ")"; - output_buf << "\n"; - } - else - octave_print_internal (output_buf, *matrix); - break; - case complex_scalar_constant: - octave_print_internal (output_buf, *complex_scalar); - break; - case complex_matrix_constant: - if (nr == 0 || nc == 0) - { - output_buf << "[]"; - if (user_pref.print_empty_dimensions) - output_buf << "(" << nr << "x" << nc << ")"; - output_buf << "\n"; - } - else - octave_print_internal (output_buf, *complex_matrix); - break; - case string_constant: - output_buf << string << "\n"; - break; - case range_constant: - octave_print_internal (output_buf, *range); - break; - case magic_colon: - default: - panic_impossible (); - break; - } - - output_buf << ends; - maybe_page_output (output_buf); - } -} - -tree_constant * -tree_constant_rep::eval (const tree_constant *args, int nargin, int nargout, - int print) -{ - if (error_state) - return NULL_TREE_CONST; - - tree_constant *retval = new tree_constant [2]; - switch (type_tag) - { - case complex_scalar_constant: - case scalar_constant: - retval[0] = do_scalar_index (args, nargin); - break; - case complex_matrix_constant: - case matrix_constant: - retval[0] = do_matrix_index (args, nargin); - break; - case string_constant: - gripe_string_invalid (); -// retval[0] = do_string_index (args, nargin); - break; - case magic_colon: - case range_constant: -// This isn\'t great, but it\'s easier than implementing a lot of -// range indexing functions. - force_numeric (); - assert (type_tag != magic_colon && type_tag != range_constant); - return eval (args, nargin, nargout, print); - break; - default: - panic_impossible (); - break; - } - - if (retval[0].is_defined ()) - retval[0].eval (print); - return retval; -} - -int -tree_constant_rep::save (ostream& os, int mark_as_global, int precision) -{ - switch (type_tag) - { - case scalar_constant: - case matrix_constant: - case complex_scalar_constant: - case complex_matrix_constant: - case string_constant: - case range_constant: - if (mark_as_global) - os << "# type: global "; - else - os << "# type: "; - break; - case magic_colon: - default: - break; - } - - long old_precision = os.precision (); - os.precision (precision); - - switch (type_tag) - { - case scalar_constant: - os << "scalar\n" - << scalar << "\n"; - break; - case matrix_constant: - os << "matrix\n" - << "# rows: " << rows () << "\n" - << "# columns: " << columns () << "\n" - << *matrix ; - break; - case complex_scalar_constant: - os << "complex scalar\n" - << *complex_scalar << "\n"; - break; - case complex_matrix_constant: - os << "complex matrix\n" - << "# rows: " << rows () << "\n" - << "# columns: " << columns () << "\n" - << *complex_matrix ; - break; - case string_constant: - os << "string\n" - << "# length: " << strlen (string) << "\n" - << string << "\n"; - break; - case range_constant: - { - os << "range\n" - << "# base, limit, increment\n" - << range->base () << " " - << range->limit () << " " - << range->inc () << "\n"; - } - break; - case magic_colon: - default: - panic_impossible (); - break; - } - - os.precision (old_precision); - -// Really want to return 1 only if write is successful. - return 1; -} - -int -tree_constant_rep::save_three_d (ostream& os, int parametric) -{ - int nr = rows (); - int nc = columns (); - - switch (type_tag) - { - case matrix_constant: - os << "# 3D data...\n" - << "# type: matrix\n" - << "# total rows: " << nr << "\n" - << "# total columns: " << nc << "\n"; - - if (parametric) - { - int extras = nc % 3; - if (extras) - warning ("ignoring last %d columns", extras); - - for (int i = 0; i < nc-extras; i += 3) - { - os << matrix->extract (0, i, nr-1, i+2); - if (i+3 < nc-extras) - os << "\n"; - } - } - else - { - for (int i = 0; i < nc; i++) - { - os << matrix->extract (0, i, nr-1, i); - if (i+1 < nc) - os << "\n"; - } - } - break; - default: - ::error ("for now, I can only save real matrices in 3D format"); - return 0; - break; - } -// Really want to return 1 only if write is successful. - return 1; -} - -int -tree_constant_rep::load (istream& is) -{ - int is_global = 0; - - type_tag = unknown_constant; - -// Look for type keyword - - char *tag = extract_keyword (is, "type"); - - if (tag != (char *) NULL && *tag != '\0') - { - char *ptr = strchr (tag, ' '); - if (ptr != (char *) NULL) - { - *ptr = '\0'; - is_global = (strncmp (tag, "global", 6) == 0); - *ptr = ' '; - if (is_global) - ptr++; - else - ptr = tag; - } - else - ptr = tag; - - if (strncmp (ptr, "scalar", 6) == 0) - type_tag = load (is, scalar_constant); - else if (strncmp (ptr, "matrix", 6) == 0) - type_tag = load (is, matrix_constant); - else if (strncmp (ptr, "complex scalar", 14) == 0) - type_tag = load (is, complex_scalar_constant); - else if (strncmp (ptr, "complex matrix", 14) == 0) - type_tag = load (is, complex_matrix_constant); - else if (strncmp (ptr, "string", 6) == 0) - type_tag = load (is, string_constant); - else if (strncmp (ptr, "range", 5) == 0) - type_tag = load (is, range_constant); - else - ::error ("unknown constant type `%s'", tag); - } - else - ::error ("failed to extract keyword specifying value type"); - - delete [] tag; - - return is_global; -} - -tree_constant_rep::constant_type -tree_constant_rep::load (istream& is, tree_constant_rep::constant_type t) -{ - tree_constant_rep::constant_type status = unknown_constant; - - switch (t) - { - case scalar_constant: - is >> scalar; - if (is) - status = scalar_constant; - else - ::error ("failed to load scalar constant"); - break; - case matrix_constant: - { - int nr = 0, nc = 0; - - if (extract_keyword (is, "rows", nr) && nr > 0 - && extract_keyword (is, "columns", nc) && nc > 0) - { - matrix = new Matrix (nr, nc); - is >> *matrix; - if (is) - status = matrix_constant; - else - ::error ("failed to load matrix constant"); - } - else - ::error ("failed to extract number of rows and columns"); - } - break; - case complex_scalar_constant: - complex_scalar = new Complex; - is >> *complex_scalar; - if (is) - status = complex_scalar_constant; - else - ::error ("failed to load complex scalar constant"); - break; - case complex_matrix_constant: - { - int nr = 0, nc = 0; - - if (extract_keyword (is, "rows", nr) && nr > 0 - && extract_keyword (is, "columns", nc) && nc > 0) - { - complex_matrix = new ComplexMatrix (nr, nc); - is >> *complex_matrix; - if (is) - status = complex_matrix_constant; - else - ::error ("failed to load complex matrix constant"); - } - else - ::error ("failed to extract number of rows and columns"); - } - break; - case string_constant: - { - int len; - if (extract_keyword (is, "length", len) && len > 0) - { - string = new char [len+1]; - is.get (string, len+1, EOF); - if (is) - status = string_constant; - else - ::error ("failed to load string constant"); - } - else - ::error ("failed to extract string length"); - } - break; - case range_constant: - skip_comments (is); - range = new Range (); - is >> *range; - if (is) - status = range_constant; - else - ::error ("failed to load range constant"); - break; - default: - panic_impossible (); - break; - } - return status; -} - -double -tree_constant_rep::double_value (void) const -{ - switch (type_tag) - { - case scalar_constant: - return scalar; - case complex_scalar_constant: - { - int flag = user_pref.ok_to_lose_imaginary_part; - if (flag == -1) - warning ("implicit conversion of complex value to real value"); - - if (flag != 0) - return ::real (*complex_scalar); - - ::error ("implicit conversion of complex value to real value"); - ::error ("not allowed"); - jump_to_top_level (); - } - default: - panic_impossible (); - break; - } -} - -Matrix -tree_constant_rep::matrix_value (void) const -{ - switch (type_tag) - { - case scalar_constant: - return Matrix (1, 1, scalar); - case matrix_constant: - return *matrix; - case complex_scalar_constant: - case complex_matrix_constant: - { - int flag = user_pref.ok_to_lose_imaginary_part; - if (flag == -1) - warning ("implicit conversion of complex matrix to real matrix"); - - if (flag != 0) - { - if (type_tag == complex_scalar_constant) - return Matrix (1, 1, ::real (*complex_scalar)); - else if (type_tag == complex_matrix_constant) - return ::real (*complex_matrix); - else - panic_impossible (); - } - else - { - ::error ("implicit conversion of complex matrix to real matrix"); - ::error ("not allowed"); - } - jump_to_top_level (); - } - default: - panic_impossible (); - break; - } -} - -Complex -tree_constant_rep::complex_value (void) const -{ - switch (type_tag) - { - case complex_scalar_constant: - return *complex_scalar; - case scalar_constant: - return Complex (scalar); - default: - panic_impossible (); - break; - } -} - -ComplexMatrix -tree_constant_rep::complex_matrix_value (void) const -{ - switch (type_tag) - { - case scalar_constant: - { - return ComplexMatrix (1, 1, Complex (scalar)); - } - case complex_scalar_constant: - { - return ComplexMatrix (1, 1, *complex_scalar); - } - case matrix_constant: - { - return ComplexMatrix (*matrix); - } - case complex_matrix_constant: - return *complex_matrix; - break; - default: - panic_impossible (); - break; - } -} - -char * -tree_constant_rep::string_value (void) const -{ - assert (type_tag == string_constant); - return string; -} - -Range -tree_constant_rep::range_value (void) const -{ - assert (type_tag == range_constant); - return *range; -} - -int -tree_constant_rep::rows (void) const -{ - int retval = -1; - switch (type_tag) - { - case scalar_constant: - case complex_scalar_constant: - retval = 1; - break; - case string_constant: - case range_constant: - retval = (columns () > 0) ? 1 : 0; - break; - case matrix_constant: - retval = matrix->rows (); - break; - case complex_matrix_constant: - retval = complex_matrix->rows (); - break; - case magic_colon: - ::error ("invalid use of colon operator"); - break; - case unknown_constant: - retval = 0; - break; - default: - panic_impossible (); - break; - } - return retval; -} - -int -tree_constant_rep::columns (void) const -{ - int retval = -1; - switch (type_tag) - { - case scalar_constant: - case complex_scalar_constant: - retval = 1; - break; - case matrix_constant: - retval = matrix->columns (); - break; - case complex_matrix_constant: - retval = complex_matrix->columns (); - break; - case string_constant: - retval = strlen (string); - break; - case range_constant: - retval = range->nelem (); - break; - case magic_colon: - ::error ("invalid use of colon operator"); - break; - case unknown_constant: - retval = 0; - break; - default: - panic_impossible (); - break; - } - return retval; -} - -tree_constant -tree_constant_rep::all (void) const -{ - if (type_tag == string_constant || type_tag == range_constant) - { - tree_constant tmp = make_numeric (); - return tmp.all (); - } - - tree_constant retval; - switch (type_tag) - { - case scalar_constant: - { - double status = (scalar != 0.0); - retval = tree_constant (status); - } - break; - case matrix_constant: - { - Matrix m = matrix->all (); - retval = tree_constant (m); - } - break; - case complex_scalar_constant: - { - double status = (*complex_scalar != 0.0); - retval = tree_constant (status); - } - break; - case complex_matrix_constant: - { - Matrix m = complex_matrix->all (); - retval = tree_constant (m); - } - break; - case string_constant: - case range_constant: - case magic_colon: - default: - panic_impossible (); - break; - } - return retval; -} - -tree_constant -tree_constant_rep::any (void) const -{ - if (type_tag == string_constant || type_tag == range_constant) - { - tree_constant tmp = make_numeric (); - return tmp.any (); - } - - tree_constant retval; - switch (type_tag) - { - case scalar_constant: - { - double status = (scalar != 0.0); - retval = tree_constant (status); - } - break; - case matrix_constant: - { - Matrix m = matrix->any (); - retval = tree_constant (m); - } - break; - case complex_scalar_constant: - { - double status = (*complex_scalar != 0.0); - retval = tree_constant (status); - } - break; - case complex_matrix_constant: - { - Matrix m = complex_matrix->any (); - retval = tree_constant (m); - } - break; - case string_constant: - case range_constant: - case magic_colon: - default: - panic_impossible (); - break; - } - return retval; -} - -tree_constant -tree_constant_rep::isstr (void) const -{ - double status = 0.0; - if (type_tag == string_constant) - status = 1.0; - tree_constant retval (status); - return retval; -} - -tree_constant -tree_constant_rep::convert_to_str (void) -{ - tree_constant retval; - - switch (type_tag) - { - case complex_scalar_constant: - case scalar_constant: - { - double d = double_value (); - int i = NINT (d); -// Warn about out of range conversions? - char s[2]; - s[0] = (char) i; - s[1] = '\0'; - retval = tree_constant (s); - } - break; - case complex_matrix_constant: - case matrix_constant: - { - ColumnVector v = to_vector (); - int len = v.length (); - if (len == 0) - ::error ("can only convert vectors and scalars to strings"); - else - { - char *s = new char [len+1]; - s[len] = '\0'; - for (int i = 0; i < len; i++) - { - double d = v.elem (i); - int ival = NINT (d); -// Warn about out of range conversions? - s[i] = (char) ival; - } - retval = tree_constant (s); - delete [] s; - } - } - break; - case range_constant: - { - Range r = range_value (); - double b = r.base (); - double incr = r.inc (); - int nel = r.nelem (); - char *s = new char [nel+1]; - s[nel] = '\0'; - for (int i = 0; i < nel; i++) - { - double d = b + i * incr; - int ival = NINT (d); -// Warn about out of range conversions? - s[i] = (char) ival; - } - retval = tree_constant (s); - delete [] s; - } - break; - case string_constant: - retval = tree_constant (*this); - break; - case magic_colon: - default: - panic_impossible (); - break; - } - return retval; -} - -void -tree_constant_rep::convert_to_row_or_column_vector (void) -{ - assert (type_tag == matrix_constant || type_tag == complex_matrix_constant); - - int nr = rows (); - int nc = columns (); - - int len = nr * nc; - - assert (len > 0); - - int new_nr = 1; - int new_nc = 1; - - if (user_pref.prefer_column_vectors) - new_nr = len; - else - new_nc = len; - - if (type_tag == matrix_constant) - { - Matrix *m = new Matrix (new_nr, new_nc); - - double *cop_out = matrix->fortran_vec (); - - for (int i = 0; i < len; i++) - { - if (new_nr == 1) - m->elem (0, i) = *cop_out++; - else - m->elem (i, 0) = *cop_out++; - } - - delete matrix; - matrix = m; - } - else - { - ComplexMatrix *cm = new ComplexMatrix (new_nr, new_nc); - - Complex *cop_out = complex_matrix->fortran_vec (); - - for (int i = 0; i < len; i++) - { - if (new_nr == 1) - cm->elem (0, i) = *cop_out++; - else - cm->elem (i, 0) = *cop_out++; - } - - delete complex_matrix; - complex_matrix = cm; - } -} - -int -tree_constant_rep::is_true (void) const -{ - if (type_tag == string_constant || type_tag == range_constant) - { - tree_constant tmp = make_numeric (); - return tmp.is_true (); - } - - int retval; - switch (type_tag) - { - case scalar_constant: - retval = (scalar != 0.0); - break; - case matrix_constant: - { - Matrix m = (matrix->all ()) . all (); - retval = (m.rows () == 1 - && m.columns () == 1 - && m.elem (0, 0) != 0.0); - } - break; - case complex_scalar_constant: - retval = (*complex_scalar != 0.0); - break; - case complex_matrix_constant: - { - Matrix m = (complex_matrix->all ()) . all (); - retval = (m.rows () == 1 - && m.columns () == 1 - && m.elem (0, 0) != 0.0); - } - break; - case string_constant: - case range_constant: - case magic_colon: - default: - panic_impossible (); - break; - } - return retval; -} - -tree_constant -tree_constant_rep::cumprod (void) const -{ - if (type_tag == string_constant || type_tag == range_constant) - { - tree_constant tmp = make_numeric (); - return tmp.cumprod (); - } - - tree_constant retval; - switch (type_tag) - { - case scalar_constant: - retval = tree_constant (scalar); - break; - case matrix_constant: - { - Matrix m = matrix->cumprod (); - retval = tree_constant (m); - } - break; - case complex_scalar_constant: - retval = tree_constant (*complex_scalar); - break; - case complex_matrix_constant: - { - ComplexMatrix m = complex_matrix->cumprod (); - retval = tree_constant (m); - } - break; - case string_constant: - case range_constant: - case magic_colon: - default: - panic_impossible (); - break; - } - return retval; -} - -tree_constant -tree_constant_rep::cumsum (void) const -{ - if (type_tag == string_constant || type_tag == range_constant) - { - tree_constant tmp = make_numeric (); - return tmp.cumsum (); - } - - tree_constant retval; - switch (type_tag) - { - case scalar_constant: - retval = tree_constant (scalar); - break; - case matrix_constant: - { - Matrix m = matrix->cumsum (); - retval = tree_constant (m); - } - break; - case complex_scalar_constant: - retval = tree_constant (*complex_scalar); - break; - case complex_matrix_constant: - { - ComplexMatrix m = complex_matrix->cumsum (); - retval = tree_constant (m); - } - break; - case string_constant: - case range_constant: - case magic_colon: - default: - panic_impossible (); - break; - } - return retval; -} - -tree_constant -tree_constant_rep::prod (void) const -{ - if (type_tag == string_constant || type_tag == range_constant) - { - tree_constant tmp = make_numeric (); - return tmp.prod (); - } - - tree_constant retval; - switch (type_tag) - { - case scalar_constant: - retval = tree_constant (scalar); - break; - case matrix_constant: - { - Matrix m = matrix->prod (); - retval = tree_constant (m); - } - break; - case complex_scalar_constant: - retval = tree_constant (*complex_scalar); - break; - case complex_matrix_constant: - { - ComplexMatrix m = complex_matrix->prod (); - retval = tree_constant (m); - } - break; - case string_constant: - case range_constant: - case magic_colon: - default: - panic_impossible (); - break; - } - return retval; -} - -tree_constant -tree_constant_rep::sum (void) const -{ - if (type_tag == string_constant || type_tag == range_constant) - { - tree_constant tmp = make_numeric (); - return tmp.sum (); - } - - tree_constant retval; - switch (type_tag) - { - case scalar_constant: - retval = tree_constant (scalar); - break; - case matrix_constant: - { - Matrix m = matrix->sum (); - retval = tree_constant (m); - } - break; - case complex_scalar_constant: - retval = tree_constant (*complex_scalar); - break; - case complex_matrix_constant: - { - ComplexMatrix m = complex_matrix->sum (); - retval = tree_constant (m); - } - break; - case string_constant: - case range_constant: - case magic_colon: - default: - panic_impossible (); - break; - } - return retval; -} - -tree_constant -tree_constant_rep::sumsq (void) const -{ - if (type_tag == string_constant || type_tag == range_constant) - { - tree_constant tmp = make_numeric (); - return tmp.sumsq (); - } - - tree_constant retval; - switch (type_tag) - { - case scalar_constant: - retval = tree_constant (scalar * scalar); - break; - case matrix_constant: - { - Matrix m = matrix->sumsq (); - retval = tree_constant (m); - } - break; - case complex_scalar_constant: - { - Complex c (*complex_scalar); - retval = tree_constant (c * c); - } - break; - case complex_matrix_constant: - { - ComplexMatrix m = complex_matrix->sumsq (); - retval = tree_constant (m); - } - break; - case string_constant: - case range_constant: - case magic_colon: - default: - panic_impossible (); - break; - } - return retval; -} - -static tree_constant -make_diag (const Matrix& v, int k) -{ - int nr = v.rows (); - int nc = v.columns (); - assert (nc == 1 || nr == 1); - - tree_constant retval; - - int roff = 0; - int coff = 0; - if (k > 0) - { - roff = 0; - coff = k; - } - else if (k < 0) - { - roff = -k; - coff = 0; - } - - if (nr == 1) - { - int n = nc + ABS (k); - Matrix m (n, n, 0.0); - for (int i = 0; i < nc; i++) - m.elem (i+roff, i+coff) = v.elem (0, i); - retval = tree_constant (m); - } - else - { - int n = nr + ABS (k); - Matrix m (n, n, 0.0); - for (int i = 0; i < nr; i++) - m.elem (i+roff, i+coff) = v.elem (i, 0); - retval = tree_constant (m); - } - - return retval; -} - -static tree_constant -make_diag (const ComplexMatrix& v, int k) -{ - int nr = v.rows (); - int nc = v.columns (); - assert (nc == 1 || nr == 1); - - tree_constant retval; - - int roff = 0; - int coff = 0; - if (k > 0) - { - roff = 0; - coff = k; - } - else if (k < 0) - { - roff = -k; - coff = 0; - } - - if (nr == 1) - { - int n = nc + ABS (k); - ComplexMatrix m (n, n, 0.0); - for (int i = 0; i < nc; i++) - m.elem (i+roff, i+coff) = v.elem (0, i); - retval = tree_constant (m); - } - else - { - int n = nr + ABS (k); - ComplexMatrix m (n, n, 0.0); - for (int i = 0; i < nr; i++) - m.elem (i+roff, i+coff) = v.elem (i, 0); - retval = tree_constant (m); - } - - return retval; -} - -tree_constant -tree_constant_rep::diag (void) const -{ - if (type_tag == string_constant || type_tag == range_constant) - { - tree_constant tmp = make_numeric (); - return tmp.diag (); - } - - tree_constant retval; - switch (type_tag) - { - case scalar_constant: - retval = tree_constant (scalar); - break; - case matrix_constant: - { - int nr = rows (); - int nc = columns (); - if (nr == 0 || nc == 0) - { - Matrix mtmp; - retval = tree_constant (mtmp); - } - else if (nr == 1 || nc == 1) - retval = make_diag (matrix_value (), 0); - else - { - ColumnVector v = matrix->diag (); - if (v.capacity () > 0) - retval = tree_constant (v); - } - } - break; - case complex_scalar_constant: - retval = tree_constant (*complex_scalar); - break; - case complex_matrix_constant: - { - int nr = rows (); - int nc = columns (); - if (nr == 0 || nc == 0) - { - Matrix mtmp; - retval = tree_constant (mtmp); - } - else if (nr == 1 || nc == 1) - retval = make_diag (complex_matrix_value (), 0); - else - { - ComplexColumnVector v = complex_matrix->diag (); - if (v.capacity () > 0) - retval = tree_constant (v); - } - } - break; - case string_constant: - case range_constant: - case magic_colon: - default: - panic_impossible (); - break; - } - return retval; -} - -tree_constant -tree_constant_rep::diag (const tree_constant& a) const -{ - if (type_tag == string_constant || type_tag == range_constant) - { - tree_constant tmp = make_numeric (); - return tmp.diag (a); - } - - tree_constant tmp_a = a.make_numeric (); - - tree_constant_rep::constant_type a_type = tmp_a.const_type (); - - tree_constant retval; - - switch (type_tag) - { - case scalar_constant: - if (a_type == scalar_constant) - { - int k = NINT (tmp_a.double_value ()); - int n = ABS (k) + 1; - if (k == 0) - retval = tree_constant (scalar); - else if (k > 0) - { - Matrix m (n, n, 0.0); - m.elem (0, k) = scalar; - retval = tree_constant (m); - } - else if (k < 0) - { - Matrix m (n, n, 0.0); - m.elem (-k, 0) = scalar; - retval = tree_constant (m); - } - } - break; - case matrix_constant: - if (a_type == scalar_constant) - { - int k = NINT (tmp_a.double_value ()); - int nr = rows (); - int nc = columns (); - if (nr == 0 || nc == 0) - { - Matrix mtmp; - retval = tree_constant (mtmp); - } - else if (nr == 1 || nc == 1) - retval = make_diag (matrix_value (), k); - else - { - ColumnVector d = matrix->diag (k); - retval = tree_constant (d); - } - } - else - ::error ("diag: invalid second argument"); - - break; - case complex_scalar_constant: - if (a_type == scalar_constant) - { - int k = NINT (tmp_a.double_value ()); - int n = ABS (k) + 1; - if (k == 0) - retval = tree_constant (*complex_scalar); - else if (k > 0) - { - ComplexMatrix m (n, n, 0.0); - m.elem (0, k) = *complex_scalar; - retval = tree_constant (m); - } - else if (k < 0) - { - ComplexMatrix m (n, n, 0.0); - m.elem (-k, 0) = *complex_scalar; - retval = tree_constant (m); - } - } - break; - case complex_matrix_constant: - if (a_type == scalar_constant) - { - int k = NINT (tmp_a.double_value ()); - int nr = rows (); - int nc = columns (); - if (nr == 0 || nc == 0) - { - Matrix mtmp; - retval = tree_constant (mtmp); - } - else if (nr == 1 || nc == 1) - retval = make_diag (complex_matrix_value (), k); - else - { - ComplexColumnVector d = complex_matrix->diag (k); - retval = tree_constant (d); - } - } - else - ::error ("diag: invalid second argument"); - - break; - case string_constant: - case range_constant: - case magic_colon: - default: - panic_impossible (); - break; - } - return retval; -} - -tree_constant -tree_constant_rep::mapper (Mapper_fcn& m_fcn, int print) const -{ - tree_constant retval; - - if (type_tag == string_constant || type_tag == range_constant) - { - tree_constant tmp = make_numeric (); - return tmp.mapper (m_fcn, print); - } - - switch (type_tag) - { - case scalar_constant: - if (m_fcn.can_return_complex_for_real_arg - && (scalar < m_fcn.lower_limit - || scalar > m_fcn.upper_limit)) - { - if (m_fcn.c_c_mapper != NULL) - { - Complex c = m_fcn.c_c_mapper (Complex (scalar)); - retval = tree_constant (c); - } - else - panic_impossible (); - } - else - { - if (m_fcn.d_d_mapper != NULL) - { - double d = m_fcn.d_d_mapper (scalar); - retval = tree_constant (d); - } - else - panic_impossible (); - } - break; - case matrix_constant: - if (m_fcn.can_return_complex_for_real_arg - && (any_element_less_than (*matrix, m_fcn.lower_limit) - || any_element_greater_than (*matrix, m_fcn.upper_limit))) - { - if (m_fcn.c_c_mapper != NULL) - { - ComplexMatrix cm = map (m_fcn.c_c_mapper, - ComplexMatrix (*matrix)); - retval = tree_constant (cm); - } - else - panic_impossible (); - } - else - { - if (m_fcn.d_d_mapper != NULL) - { - Matrix m = map (m_fcn.d_d_mapper, *matrix); - retval = tree_constant (m); - } - else - panic_impossible (); - } - break; - case complex_scalar_constant: - if (m_fcn.d_c_mapper != NULL) - { - double d; - d = m_fcn.d_c_mapper (*complex_scalar); - retval = tree_constant (d); - } - else if (m_fcn.c_c_mapper != NULL) - { - Complex c; - c = m_fcn.c_c_mapper (*complex_scalar); - retval = tree_constant (c); - } - else - panic_impossible (); - break; - case complex_matrix_constant: - if (m_fcn.d_c_mapper != NULL) - { - Matrix m; - m = map (m_fcn.d_c_mapper, *complex_matrix); - retval = tree_constant (m); - } - else if (m_fcn.c_c_mapper != NULL) - { - ComplexMatrix cm; - cm = map (m_fcn.c_c_mapper, *complex_matrix); - retval = tree_constant (cm); - } - else - panic_impossible (); - break; - case string_constant: - case range_constant: - case magic_colon: - default: - panic_impossible (); - break; - } - - if (retval.is_defined ()) - return retval.eval (print); - else - return retval; -} +#ifndef MIN +#define MIN(a,b) ((a) < (b) ? (a) : (b)) +#endif tree_constant::~tree_constant (void) { @@ -2565,2567 +118,1154 @@ return retval; } - -/* - * Top-level tree-constant function that handles assignments. Only - * decide if the left-hand side is currently a scalar or a matrix and - * hand off to other functions to do the real work. - */ -void -tree_constant_rep::assign (tree_constant& rhs, tree_constant *args, int nargs) +Matrix +max (const Matrix& a, const Matrix& b) { - tree_constant rhs_tmp = rhs.make_numeric (); - -// This is easier than actually handling assignments to strings. -// An assignment to a range will normally require a conversion to a -// vector since it will normally destroy the equally-spaced property -// of the range elements. - - if (type_tag == string_constant || type_tag == range_constant) - force_numeric (); + int nr = a.rows (); + int nc = a.columns (); + if (nr != b.rows () || nc != b.columns ()) + { + error ("two-arg max expecting args of same size"); + return Matrix (); + } - switch (type_tag) - { - case complex_scalar_constant: - case scalar_constant: - case unknown_constant: - do_scalar_assignment (rhs_tmp, args, nargs); - break; - case complex_matrix_constant: - case matrix_constant: - do_matrix_assignment (rhs_tmp, args, nargs); - break; - case string_constant: - ::error ("invalid assignment to string type"); - break; - case range_constant: - case magic_colon: - default: - panic_impossible (); - break; - } + Matrix result (nr, nc); + + for (int j = 0; j < nc; j++) + for (int i = 0; i < nr; i++) + { + double a_elem = a.elem (i, j); + double b_elem = b.elem (i, j); + result.elem (i, j) = MAX (a_elem, b_elem); + } + + return result; } -/* - * Assignments to scalars. If resize_on_range_error is true, - * this can convert the left-hand size to a matrix. - */ -void -tree_constant_rep::do_scalar_assignment (tree_constant& rhs, - tree_constant *args, int nargs) +ComplexMatrix +max (const ComplexMatrix& a, const ComplexMatrix& b) { - assert (type_tag == unknown_constant - || type_tag == scalar_constant - || type_tag == complex_scalar_constant); - - if ((rhs.is_scalar_type () || rhs.is_zero_by_zero ()) - && valid_scalar_indices (args, nargs)) - { - if (rhs.is_zero_by_zero ()) - { - if (type_tag == complex_scalar_constant) - delete complex_scalar; - - matrix = new Matrix (0, 0); - type_tag = matrix_constant; - } - else if (type_tag == unknown_constant || type_tag == scalar_constant) - { - if (rhs.const_type () == scalar_constant) - { - scalar = rhs.double_value (); - type_tag = scalar_constant; - } - else if (rhs.const_type () == complex_scalar_constant) - { - complex_scalar = new Complex (rhs.complex_value ()); - type_tag = complex_scalar_constant; - } - else - { - ::error ("invalid assignment to scalar"); - return; - } - } - else - { - if (rhs.const_type () == scalar_constant) - { - delete complex_scalar; - scalar = rhs.double_value (); - type_tag = scalar_constant; - } - else if (rhs.const_type () == complex_scalar_constant) - { - *complex_scalar = rhs.complex_value (); - type_tag = complex_scalar_constant; - } - else - { - ::error ("invalid assignment to scalar"); - return; - } - } - } - else if (user_pref.resize_on_range_error) + int nr = a.rows (); + int nc = a.columns (); + if (nr != b.rows () || nc != b.columns ()) { - tree_constant_rep::constant_type old_type_tag = type_tag; - - if (type_tag == complex_scalar_constant) - { - Complex *old_complex = complex_scalar; - complex_matrix = new ComplexMatrix (1, 1, *complex_scalar); - type_tag = complex_matrix_constant; - delete old_complex; - } - else if (type_tag == scalar_constant) - { - matrix = new Matrix (1, 1, scalar); - type_tag = matrix_constant; - } - -// If there is an error, the call to do_matrix_assignment should not -// destroy the current value. tree_constant_rep::eval(int) will take -// care of converting single element matrices back to scalars. - - do_matrix_assignment (rhs, args, nargs); - -// I don't think there's any other way to revert back to unknown -// constant types, so here it is. - - if (old_type_tag == unknown_constant && error_state) - { - if (type_tag == matrix_constant) - delete matrix; - else if (type_tag == complex_matrix_constant) - delete complex_matrix; - - type_tag = unknown_constant; - } - } - else if (nargs > 3 || nargs < 2) - ::error ("invalid index expression for scalar type"); - else - ::error ("index invalid or out of range for scalar type"); -} - -/* - * Assignments to matrices (and vectors). - * - * For compatibility with Matlab, we allow assignment of an empty - * matrix to an expression with empty indices to do nothing. - */ -void -tree_constant_rep::do_matrix_assignment (tree_constant& rhs, - tree_constant *args, int nargs) -{ - assert (type_tag == unknown_constant - || type_tag == matrix_constant - || type_tag == complex_matrix_constant); - - if (type_tag == matrix_constant && rhs.is_complex_type ()) - { - Matrix *old_matrix = matrix; - complex_matrix = new ComplexMatrix (*matrix); - type_tag = complex_matrix_constant; - delete old_matrix; - } - else if (type_tag == unknown_constant) - { - if (rhs.is_complex_type ()) - { - complex_matrix = new ComplexMatrix (); - type_tag = complex_matrix_constant; - } - else - { - matrix = new Matrix (); - type_tag = matrix_constant; - } + error ("two-arg max expecting args of same size"); + return ComplexMatrix (); } -// The do_matrix_assignment functions can't handle empty matrices, so -// don't let any pass through here. - switch (nargs) - { - case 2: - if (args == NULL_TREE_CONST) - ::error ("matrix index is null"); - else if (args[1].is_undefined ()) - ::error ("matrix index is undefined"); - else - do_matrix_assignment (rhs, args[1]); - break; - case 3: - if (args == NULL_TREE_CONST) - ::error ("matrix indices are null"); - else if (args[1].is_undefined ()) - ::error ("first matrix index is undefined"); - else if (args[2].is_undefined ()) - ::error ("second matrix index is undefined"); - else if (args[1].is_empty () || args[2].is_empty ()) - { - if (! rhs.is_empty ()) - { - ::error ("in assignment expression, a matrix index is empty"); - ::error ("but hte right hand side is not an empty matrix"); - } -// XXX FIXME XXX -- to really be correct here, we should probably -// check to see if the assignment conforms, but that seems like more -// work than it's worth right now... - } - else - do_matrix_assignment (rhs, args[1], args[2]); - break; - default: - ::error ("too many indices for matrix expression"); - break; - } + ComplexMatrix result (nr, nc); + + for (int j = 0; j < nc; j++) + for (int i = 0; i < nr; i++) + { + double abs_a_elem = abs (a.elem (i, j)); + double abs_b_elem = abs (b.elem (i, j)); + if (abs_a_elem > abs_b_elem) + result.elem (i, j) = a.elem (i, j); + else + result.elem (i, j) = b.elem (i, j); + } + + return result; } -/* - * Matrix assignments indexed by a single value. - */ -void -tree_constant_rep::do_matrix_assignment (tree_constant& rhs, - tree_constant& i_arg) +Matrix +min (const Matrix& a, const Matrix& b) { - int nr = rows (); - int nc = columns (); - - if (user_pref.do_fortran_indexing || nr <= 1 || nc <= 1) + int nr = a.rows (); + int nc = a.columns (); + if (nr != b.rows () || nc != b.columns ()) { - if (i_arg.is_empty ()) - { - if (! rhs.is_empty ()) - { - ::error ("in assignment expression, matrix index is empty but"); - ::error ("right hand side is not an empty matrix"); - } -// XXX FIXME XXX -- to really be correct here, we should probably -// check to see if the assignment conforms, but that seems like more -// work than it's worth right now... + error ("two-arg min expecting args of same size"); + return Matrix (); + } + + Matrix result (nr, nc); -// The assignment functions can't handle empty matrices, so don't let -// any pass through here. - return; - } - -// We can't handle the case of assigning to a vector first, since even -// then, the two operations are not equivalent. For example, the -// expression V(:) = M is handled differently depending on whether the -// user specified do_fortran_indexing = "true". + for (int j = 0; j < nc; j++) + for (int i = 0; i < nr; i++) + { + double a_elem = a.elem (i, j); + double b_elem = b.elem (i, j); + result.elem (i, j) = MIN (a_elem, b_elem); + } - if (user_pref.do_fortran_indexing) - fortran_style_matrix_assignment (rhs, i_arg); - else if (nr <= 1 || nc <= 1) - vector_assignment (rhs, i_arg); - else - panic_impossible (); - } - else - ::error ("single index only valid for row or column vector"); + return result; } -/* - * Fortran-style assignments. Matrices are assumed to be stored in - * column-major order and it is ok to use a single index for - * multi-dimensional matrices. - */ -void -tree_constant_rep::fortran_style_matrix_assignment (tree_constant& rhs, - tree_constant& i_arg) +ComplexMatrix +min (const ComplexMatrix& a, const ComplexMatrix& b) { - tree_constant tmp_i = i_arg.make_numeric_or_magic (); - - tree_constant_rep::constant_type itype = tmp_i.const_type (); - - int nr = rows (); - int nc = columns (); - - int rhs_nr = rhs.rows (); - int rhs_nc = rhs.columns (); - - switch (itype) + int nr = a.rows (); + int nc = a.columns (); + if (nr != b.rows () || nc != b.columns ()) { - case complex_scalar_constant: - case scalar_constant: - { - int i = NINT (tmp_i.double_value ()); - int idx = i - 1; - - if (rhs_nr == 0 && rhs_nc == 0) - { - if (idx < nr * nc) - { - convert_to_row_or_column_vector (); - - nr = rows (); - nc = columns (); + error ("two-arg min expecting args of same size"); + return ComplexMatrix (); + } - if (nr == 1) - delete_column (idx); - else if (nc == 1) - delete_row (idx); - else - panic_impossible (); - } - return; - } - - if (index_check (idx, "") < 0) - return; - - if (nr <= 1 || nc <= 1) - { - maybe_resize (idx); - if (error_state) - return; - } - else if (range_max_check (idx, nr * nc) < 0) - return; - - nr = rows (); - nc = columns (); - - if (! indexed_assign_conforms (1, 1, rhs_nr, rhs_nc)) - { - ::error ("for A(int) = X: X must be a scalar"); - return; - } - int ii = fortran_row (i, nr) - 1; - int jj = fortran_column (i, nr) - 1; - do_matrix_assignment (rhs, ii, jj); - } - break; - case complex_matrix_constant: - case matrix_constant: - { - Matrix mi = tmp_i.matrix_value (); - int len = nr * nc; - idx_vector ii (mi, 1, "", len); // Always do fortran indexing here... - if (! ii) - return; + ComplexMatrix result (nr, nc); - if (rhs_nr == 0 && rhs_nc == 0) - { - ii.sort_uniq (); - int num_to_delete = 0; - for (int i = 0; i < ii.length (); i++) - { - if (ii.elem (i) < len) - num_to_delete++; - else - break; - } - - if (num_to_delete > 0) - { - if (num_to_delete != ii.length ()) - ii.shorten (num_to_delete); - - convert_to_row_or_column_vector (); - - nr = rows (); - nc = columns (); - - if (nr == 1) - delete_columns (ii); - else if (nc == 1) - delete_rows (ii); - else - panic_impossible (); - } - return; - } - - if (nr <= 1 || nc <= 1) - { - maybe_resize (ii.max ()); - if (error_state) - return; - } - else if (range_max_check (ii.max (), len) < 0) - return; + for (int j = 0; j < nc; j++) + for (int i = 0; i < nr; i++) + { + double abs_a_elem = abs (a.elem (i, j)); + double abs_b_elem = abs (b.elem (i, j)); + if (abs_a_elem < abs_b_elem) + result.elem (i, j) = a.elem (i, j); + else + result.elem (i, j) = b.elem (i, j); + } - int ilen = ii.capacity (); - - if (ilen != rhs_nr * rhs_nc) - { - ::error ("A(matrix) = X: X and matrix must have the same number"); - ::error ("of elements"); - } - else if (ilen == 1 && rhs.is_scalar_type ()) - { - int nr = rows (); - int idx = ii.elem (0); - int ii = fortran_row (idx + 1, nr) - 1; - int jj = fortran_column (idx + 1, nr) - 1; - - if (rhs.const_type () == scalar_constant) - matrix->elem (ii, jj) = rhs.double_value (); - else if (rhs.const_type () == complex_scalar_constant) - complex_matrix->elem (ii, jj) = rhs.complex_value (); - else - panic_impossible (); - } - else - fortran_style_matrix_assignment (rhs, ii); - } - break; - case string_constant: - gripe_string_invalid (); - break; - case range_constant: - gripe_range_invalid (); - break; - case magic_colon: -// a(:) = [] is equivalent to a(:,:) = foo. - if (rhs_nr == 0 && rhs_nc == 0) - do_matrix_assignment (rhs, magic_colon, magic_colon); - else - fortran_style_matrix_assignment (rhs, magic_colon); - break; - default: - panic_impossible (); - break; - } + return result; } -/* - * Fortran-style assignment for vector index. - */ -void -tree_constant_rep::fortran_style_matrix_assignment (tree_constant& rhs, - idx_vector& i) +static void +get_dimensions (const tree_constant& a, const char *warn_for, + int& nr, int& nc) { - assert (rhs.is_matrix_type ()); - - int ilen = i.capacity (); - - REP_RHS_MATRIX (rhs, rhs_m, rhs_cm, rhs_nr, rhs_nc); - - int len = rhs_nr * rhs_nc; - - if (len == ilen) - { - int nr = rows (); - if (rhs.const_type () == matrix_constant) - { - double *cop_out = rhs_m.fortran_vec (); - for (int k = 0; k < len; k++) - { - int ii = fortran_row (i.elem (k) + 1, nr) - 1; - int jj = fortran_column (i.elem (k) + 1, nr) - 1; - - matrix->elem (ii, jj) = *cop_out++; - } - } - else - { - Complex *cop_out = rhs_cm.fortran_vec (); - for (int k = 0; k < len; k++) - { - int ii = fortran_row (i.elem (k) + 1, nr) - 1; - int jj = fortran_column (i.elem (k) + 1, nr) - 1; + tree_constant tmpa = a.make_numeric (); - complex_matrix->elem (ii, jj) = *cop_out++; - } - } - } - else - ::error ("number of rows and columns must match for indexed assignment"); -} - -/* - * Fortran-style assignment for colon index. - */ -void -tree_constant_rep::fortran_style_matrix_assignment - (tree_constant& rhs, tree_constant_rep::constant_type mci) -{ - assert (rhs.is_matrix_type () && mci == tree_constant_rep::magic_colon); - - int nr = rows (); - int nc = columns (); - - REP_RHS_MATRIX (rhs, rhs_m, rhs_cm, rhs_nr, rhs_nc); - - int rhs_size = rhs_nr * rhs_nc; - if (rhs_size == 0) + if (tmpa.is_scalar_type ()) { - if (rhs.const_type () == matrix_constant) - { - delete matrix; - matrix = new Matrix (0, 0); - return; - } - else - panic_impossible (); - } - else if (nr*nc != rhs_size) - { - ::error ("A(:) = X: X and A must have the same number of elements"); - return; - } - - if (rhs.const_type () == matrix_constant) - { - double *cop_out = rhs_m.fortran_vec (); - for (int j = 0; j < nc; j++) - for (int i = 0; i < nr; i++) - matrix->elem (i, j) = *cop_out++; + double tmp = tmpa.double_value (); + nr = nc = NINT (tmp); } else { - Complex *cop_out = rhs_cm.fortran_vec (); - for (int j = 0; j < nc; j++) - for (int i = 0; i < nr; i++) - complex_matrix->elem (i, j) = *cop_out++; + nr = tmpa.rows (); + nc = tmpa.columns (); + + if ((nr == 1 && nc == 2) || (nr == 2 && nc == 1)) + { + ColumnVector v = tmpa.to_vector (); + + nr = NINT (v.elem (0)); + nc = NINT (v.elem (1)); + } + else + warning ("%s (A): use %s (size (A)) instead", warn_for, warn_for); + } + + check_dimensions (nr, nc, warn_for); // May set error_state. +} + +static void +get_dimensions (const tree_constant& a, const tree_constant& b, + const char *warn_for, int& nr, int& nc) +{ + tree_constant tmpa = a.make_numeric (); + tree_constant tmpb = b.make_numeric (); + + if (tmpa.is_scalar_type () && tmpb.is_scalar_type ()) + { + nr = NINT (tmpa.double_value ()); + nc = NINT (tmpb.double_value ()); + + check_dimensions (nr, nc, warn_for); // May set error_state. } + else + error ("%s: expecting two scalar arguments", warn_for); +} + +tree_constant +fill_matrix (const tree_constant& a, double val, const char *warn_for) +{ + int nr, nc; + get_dimensions (a, warn_for, nr, nc); + + if (error_state) + return tree_constant (); + + Matrix m (nr, nc, val); + + return tree_constant (m); +} + +tree_constant +fill_matrix (const tree_constant& a, const tree_constant& b, + double val, const char *warn_for) +{ + int nr, nc; + get_dimensions (a, b, warn_for, nr, nc); // May set error_state. + + if (error_state) + return tree_constant (); + + Matrix m (nr, nc, val); + + return tree_constant (m); } -/* - * Assignments to vectors. Hand off to other functions once we know - * what kind of index we have. For a colon, it is the same as - * assignment to a matrix indexed by two colons. - */ -void -tree_constant_rep::vector_assignment (tree_constant& rhs, tree_constant& i_arg) +tree_constant +identity_matrix (const tree_constant& a) { - int nr = rows (); - int nc = columns (); + int nr, nc; + get_dimensions (a, "eye", nr, nc); // May set error_state. + + if (error_state) + return tree_constant (); + + Matrix m (nr, nc, 0.0); - assert ((nr == 1 || nc == 1 || (nr == 0 && nc == 0)) - && ! user_pref.do_fortran_indexing); + if (nr > 0 && nc > 0) + { + int n = MIN (nr, nc); + for (int i = 0; i < n; i++) + m.elem (i, i) = 1.0; + } - tree_constant tmp_i = i_arg.make_numeric_or_range_or_magic (); - - tree_constant_rep::constant_type itype = tmp_i.const_type (); + return tree_constant (m); +} - switch (itype) +tree_constant +identity_matrix (const tree_constant& a, const tree_constant& b) +{ + int nr, nc; + get_dimensions (a, b, "eye", nr, nc); // May set error_state. + + if (error_state) + return tree_constant (); + + Matrix m (nr, nc, 0.0); + + if (nr > 0 && nc > 0) { - case complex_scalar_constant: - case scalar_constant: - { - int i = tree_to_mat_idx (tmp_i.double_value ()); - if (index_check (i, "") < 0) - return; - do_vector_assign (rhs, i); - } - break; - case complex_matrix_constant: - case matrix_constant: - { - Matrix mi = tmp_i.matrix_value (); - int len = nr * nc; - idx_vector iv (mi, user_pref.do_fortran_indexing, "", len); - if (! iv) - return; + int n = MIN (nr, nc); + for (int i = 0; i < n; i++) + m.elem (i, i) = 1.0; + } + + return tree_constant (m); +} - do_vector_assign (rhs, iv); - } - break; - case string_constant: - gripe_string_invalid (); - break; - case range_constant: - { - Range ri = tmp_i.range_value (); - int len = nr * nc; - if (len == 2 && is_zero_one (ri)) - { - do_vector_assign (rhs, 1); - } - else if (len == 2 && is_one_zero (ri)) - { - do_vector_assign (rhs, 0); - } - else +static tree_constant +find_nonzero_elem_idx (const Matrix& m) +{ + int count = 0; + int m_nr = m.rows (); + int m_nc = m.columns (); + + int i; + for (int j = 0; j < m_nc; j++) + for (i = 0; i < m_nr; i++) + if (m.elem (i, j) != 0) + count++; + + Matrix result; + + if (count == 0) + return result; + + if (m_nr == 1) + { + result.resize (1, count); + count = 0; + for (j = 0; j < m_nc; j++) + if (m.elem (0, j) != 0) { - if (index_check (ri, "") < 0) - return; - do_vector_assign (rhs, ri); + result (0, count) = j + 1; + count++; } - } - break; - case magic_colon: - { - int rhs_nr = rhs.rows (); - int rhs_nc = rhs.columns (); - - if (! indexed_assign_conforms (nr, nc, rhs_nr, rhs_nc)) - { - ::error ("A(:) = X: X and A must have the same dimensions"); - return; - } - do_matrix_assignment (rhs, magic_colon, magic_colon); - } - break; - default: - panic_impossible (); - break; + return tree_constant (result); + } + else + { + ColumnVector v (count); + count = 0; + for (j = 0; j < m_nc; j++) + for (i = 0; i < m_nr; i++) + if (m.elem (i, j) != 0) + { + v.elem (count) = m_nr * j + i + 1; + count++; + } + return tree_constant (v, 1); // Always make a column vector. } } -/* - * Check whether an indexed assignment to a vector is valid. - */ -void -tree_constant_rep::check_vector_assign (int rhs_nr, int rhs_nc, - int ilen, const char *rm) -{ - int nr = rows (); - int nc = columns (); - - if ((nr == 1 && nc == 1) || nr == 0 || nc == 0) // No orientation. - { - if (! (ilen == rhs_nr || ilen == rhs_nc)) - { - ::error ("A(%s) = X: X and %s must have the same number of elements", - rm, rm); - } - } - else if (nr == 1) // Preserve current row orientation. - { - if (! (rhs_nr == 1 && rhs_nc == ilen)) - { - ::error ("A(%s) = X: where A is a row vector, X must also be a", rm); - ::error ("row vector with the same number of elements as %s", rm); - } - } - else if (nc == 1) // Preserve current column orientation. - { - if (! (rhs_nc == 1 && rhs_nr == ilen)) - { - ::error ("A(%s) = X: where A is a column vector, X must also be", rm); - ::error ("a column vector with the same number of elements as %s", rm); - } - } - else - panic_impossible (); -} - -/* - * Assignment to a vector with an integer index. - */ -void -tree_constant_rep::do_vector_assign (tree_constant& rhs, int i) +static tree_constant +find_nonzero_elem_idx (const ComplexMatrix& m) { - int rhs_nr = rhs.rows (); - int rhs_nc = rhs.columns (); + int count = 0; + int m_nr = m.rows (); + int m_nc = m.columns (); - if (indexed_assign_conforms (1, 1, rhs_nr, rhs_nc)) + for (int j = 0; j < m_nc; j++) { - maybe_resize (i); - if (error_state) - return; + for (int i = 0; i < m_nr; i++) + if (m.elem (i, j) != 0) + count++; + } - int nr = rows (); - int nc = columns (); + Matrix result; - if (nr == 1) - { - REP_ELEM_ASSIGN (0, i, rhs.double_value (), rhs.complex_value (), - rhs.is_real_type ()); - } - else if (nc == 1) - { - REP_ELEM_ASSIGN (i, 0, rhs.double_value (), rhs.complex_value (), - rhs.is_real_type ()); - } - else - panic_impossible (); - } - else if (rhs_nr == 0 && rhs_nc == 0) + if (count == 0) + return result; + + if (m_nr == 1) { - int nr = rows (); - int nc = columns (); - - int len = nr > nc ? nr : nc; - - if (i < 0 || i >= len) - { - ::error ("A(int) = []: index out of range"); - return; - } - - if (nr == 1) - delete_column (i); - else if (nc == 1) - delete_row (i); - else - panic_impossible (); + result.resize (1, count); + count = 0; + for (j = 0; j < m_nc; j++) + if (m.elem (0, j) != 0) + { + result (0, count) = j + 1; + count++; + } + return tree_constant (result); } else { - ::error ("for A(int) = X: X must be a scalar"); - return; + ColumnVector v (count); + count = 0; + for (j = 0; j < m_nc; j++) + { + for (int i = 0; i < m_nr; i++) + if (m.elem (i, j) != 0) + { + v.elem (count) = m_nr * j + i + 1; + count++; + } + } + return tree_constant (v, 1); // Always make a column vector. } } -/* - * Assignment to a vector with a vector index. - */ -void -tree_constant_rep::do_vector_assign (tree_constant& rhs, idx_vector& iv) +tree_constant +find_nonzero_elem_idx (const tree_constant& a) { - if (rhs.is_zero_by_zero ()) - { - int nr = rows (); - int nc = columns (); - - int len = nr > nc ? nr : nc; + tree_constant retval; - if (iv.max () >= len) - { - ::error ("A(matrix) = []: index out of range"); - return; - } - - if (nr == 1) - delete_columns (iv); - else if (nc == 1) - delete_rows (iv); - else - panic_impossible (); - } - else if (rhs.is_scalar_type ()) - { - int nr = rows (); - int nc = columns (); - - if (iv.capacity () == 1) - { - int idx = iv.elem (0); + tree_constant tmp = a.make_numeric (); - if (nr == 1) - { - REP_ELEM_ASSIGN (0, idx, rhs.double_value (), - rhs.complex_value (), rhs.is_real_type ()); - } - else if (nc == 1) - { - REP_ELEM_ASSIGN (idx, 0, rhs.double_value (), - rhs.complex_value (), rhs.is_real_type ()); - } - else - panic_impossible (); - } - else - { - if (nr == 1) - { - ::error ("A(matrix) = X: where A is a row vector, X must also be a"); - ::error ("row vector with the same number of elements as matrix"); - } - else if (nc == 1) - { - ::error ("A(matrix) = X: where A is a column vector, X must also be a"); - ::error ("column vector with the same number of elements as matrix"); - } - else - panic_impossible (); - } - } - else if (rhs.is_matrix_type ()) + Matrix result; + + switch (tmp.const_type ()) { - REP_RHS_MATRIX (rhs, rhs_m, rhs_cm, rhs_nr, rhs_nc); - - int ilen = iv.capacity (); - check_vector_assign (rhs_nr, rhs_nc, ilen, "matrix"); - if (error_state) - return; - - force_orient f_orient = no_orient; - if (rhs_nr == 1 && rhs_nc != 1) - f_orient = row_orient; - else if (rhs_nc == 1 && rhs_nr != 1) - f_orient = column_orient; - - maybe_resize (iv.max (), f_orient); - if (error_state) - return; - - int nr = rows (); - int nc = columns (); - - if (nr == 1) - { - for (int i = 0; i < iv.capacity (); i++) - REP_ELEM_ASSIGN (0, iv.elem (i), rhs_m.elem (0, i), - rhs_cm.elem (0, i), rhs.is_real_type ()); - } - else if (nc == 1) - { - for (int i = 0; i < iv.capacity (); i++) - REP_ELEM_ASSIGN (iv.elem (i), 0, rhs_m.elem (i, 0), - rhs_cm.elem (i, 0), rhs.is_real_type ()); - } - else - panic_impossible (); + case tree_constant_rep::matrix_constant: + { + Matrix m = tmp.matrix_value (); + return find_nonzero_elem_idx (m); + } + break; + case tree_constant_rep::scalar_constant: + { + double d = tmp.double_value (); + if (d != 0.0) + return tree_constant (1.0); + else + return tree_constant (result); + } + break; + case tree_constant_rep::complex_matrix_constant: + { + ComplexMatrix m = tmp.complex_matrix_value (); + return find_nonzero_elem_idx (m); + } + break; + case tree_constant_rep::complex_scalar_constant: + { + Complex c = tmp.complex_value (); + if (c != 0.0) + return tree_constant (1.0); + else + return tree_constant (result); + } + break; + default: + break; } - else - panic_impossible (); + return retval; } -/* - * Assignment to a vector with a range index. - */ -void -tree_constant_rep::do_vector_assign (tree_constant& rhs, Range& ri) +// XXX FIXME XXX -- the next two functions (and expm) should really be just +// one... + +tree_constant * +matrix_log (const tree_constant& a) { - if (rhs.is_zero_by_zero ()) - { - int nr = rows (); - int nc = columns (); - - int len = nr > nc ? nr : nc; - - int b = tree_to_mat_idx (ri.min ()); - int l = tree_to_mat_idx (ri.max ()); - if (b < 0 || l >= len) - { - ::error ("A(range) = []: index out of range"); - return; - } + tree_constant *retval = new tree_constant [2]; - if (nr == 1) - delete_columns (ri); - else if (nc == 1) - delete_rows (ri); - else - panic_impossible (); - } - else if (rhs.is_scalar_type ()) + tree_constant tmp = a.make_numeric ();; + + if (tmp.rows () == 0 || tmp.columns () == 0) { - int nr = rows (); - int nc = columns (); - - if (nr == 1) + int flag = user_pref.propagate_empty_matrices; + if (flag != 0) { - ::error ("A(range) = X: where A is a row vector, X must also be a"); - ::error ("row vector with the same number of elements as range"); - } - else if (nc == 1) - { - ::error ("A(range) = X: where A is a column vector, X must also be a"); - ::error ("column vector with the same number of elements as range"); + if (flag < 0) + gripe_empty_arg ("logm", 0); + Matrix m; + retval = new tree_constant [2]; + retval[0] = tree_constant (m); + return retval; } else - panic_impossible (); + gripe_empty_arg ("logm", 1); } - else if (rhs.is_matrix_type ()) - { - REP_RHS_MATRIX (rhs, rhs_m, rhs_cm, rhs_nr, rhs_nc); - - int ilen = ri.nelem (); - check_vector_assign (rhs_nr, rhs_nc, ilen, "range"); - if (error_state) - return; - force_orient f_orient = no_orient; - if (rhs_nr == 1 && rhs_nc != 1) - f_orient = row_orient; - else if (rhs_nc == 1 && rhs_nr != 1) - f_orient = column_orient; + switch (tmp.const_type ()) + { + case tree_constant_rep::matrix_constant: + { + Matrix m = tmp.matrix_value (); - maybe_resize (tree_to_mat_idx (ri.max ()), f_orient); - if (error_state) - return; - - int nr = rows (); - int nc = columns (); - - double b = ri.base (); - double increment = ri.inc (); + int nr = m.rows (); + int nc = m.columns (); - if (nr == 1) - { - for (int i = 0; i < ri.nelem (); i++) - { - double tmp = b + i * increment; - int col = tree_to_mat_idx (tmp); - REP_ELEM_ASSIGN (0, col, rhs_m.elem (0, i), rhs_cm.elem (0, i), - rhs.is_real_type ()); - } - } - else if (nc == 1) - { - for (int i = 0; i < ri.nelem (); i++) - { - double tmp = b + i * increment; - int row = tree_to_mat_idx (tmp); - REP_ELEM_ASSIGN (row, 0, rhs_m.elem (i, 0), rhs_cm.elem (i, 0), - rhs.is_real_type ()); - } - } - else - panic_impossible (); - } - else - panic_impossible (); -} + if (nr == 0 || nc == 0 || nr != nc) + gripe_square_matrix_required ("logm"); + else + { + EIG m_eig (m); + ComplexColumnVector lambda (m_eig.eigenvalues ()); + ComplexMatrix Q (m_eig.eigenvectors ()); -/* - * Matrix assignment indexed by two values. This function determines - * the type of the first arugment, checks as much as possible, and - * then calls one of a set of functions to handle the specific cases: - * - * M (integer, arg2) = RHS (MA1) - * M (vector, arg2) = RHS (MA2) - * M (range, arg2) = RHS (MA3) - * M (colon, arg2) = RHS (MA4) - * - * Each of those functions determines the type of the second argument - * and calls another function to handle the real work of doing the - * assignment. - */ -void -tree_constant_rep::do_matrix_assignment (tree_constant& rhs, - tree_constant& i_arg, - tree_constant& j_arg) -{ - tree_constant tmp_i = i_arg.make_numeric_or_range_or_magic (); + for (int i = 0; i < nr; i++) + { + Complex elt = lambda.elem (i); + if (imag (elt) == 0.0 && real (elt) > 0.0) + lambda.elem (i) = log (real (elt)); + else + lambda.elem (i) = log (elt); + } - tree_constant_rep::constant_type itype = tmp_i.const_type (); + ComplexDiagMatrix D (lambda); + ComplexMatrix result = Q * D * Q.inverse (); - switch (itype) - { - case complex_scalar_constant: - case scalar_constant: - { - int i = tree_to_mat_idx (tmp_i.double_value ()); - if (index_check (i, "row") < 0) - return; - do_matrix_assignment (rhs, i, j_arg); + retval[0] = tree_constant (result); + } } break; - case complex_matrix_constant: - case matrix_constant: + case tree_constant_rep::complex_matrix_constant: { - Matrix mi = tmp_i.matrix_value (); - idx_vector iv (mi, user_pref.do_fortran_indexing, "row", rows ()); - if (! iv) - return; + ComplexMatrix m = tmp.complex_matrix_value (); - do_matrix_assignment (rhs, iv, j_arg); - } - break; - case string_constant: - gripe_string_invalid (); - break; - case range_constant: - { - Range ri = tmp_i.range_value (); - int nr = rows (); - if (nr == 2 && is_zero_one (ri)) - { - do_matrix_assignment (rhs, 1, j_arg); - } - else if (nr == 2 && is_one_zero (ri)) - { - do_matrix_assignment (rhs, 0, j_arg); - } + int nr = m.rows (); + int nc = m.columns (); + + if (nr == 0 || nc == 0 || nr != nc) + gripe_square_matrix_required ("logm"); else { - if (index_check (ri, "row") < 0) - return; - do_matrix_assignment (rhs, ri, j_arg); + EIG m_eig (m); + ComplexColumnVector lambda (m_eig.eigenvalues ()); + ComplexMatrix Q (m_eig.eigenvectors ()); + + for (int i = 0; i < nr; i++) + { + Complex elt = lambda.elem (i); + if (imag (elt) == 0.0 && real (elt) > 0.0) + lambda.elem (i) = log (real (elt)); + else + lambda.elem (i) = log (elt); + } + + ComplexDiagMatrix D (lambda); + ComplexMatrix result = Q * D * Q.inverse (); + + retval[0] = tree_constant (result); } } break; - case magic_colon: - do_matrix_assignment (rhs, magic_colon, j_arg); - break; - default: - panic_impossible (); - break; - } -} - -/* MA1 */ -void -tree_constant_rep::do_matrix_assignment (tree_constant& rhs, int i, - tree_constant& j_arg) -{ - tree_constant tmp_j = j_arg.make_numeric_or_range_or_magic (); - - tree_constant_rep::constant_type jtype = tmp_j.const_type (); - - int rhs_nr = rhs.rows (); - int rhs_nc = rhs.columns (); - - switch (jtype) - { - case complex_scalar_constant: - case scalar_constant: + case tree_constant_rep::scalar_constant: { - int j = tree_to_mat_idx (tmp_j.double_value ()); - if (index_check (j, "column") < 0) - return; - if (! indexed_assign_conforms (1, 1, rhs_nr, rhs_nc)) - { - ::error ("A(int,int) = X, X must be a scalar"); - return; - } - maybe_resize (i, j); - if (error_state) - return; - - do_matrix_assignment (rhs, i, j); - } - break; - case complex_matrix_constant: - case matrix_constant: - { - Matrix mj = tmp_j.matrix_value (); - idx_vector jv (mj, user_pref.do_fortran_indexing, "column", - columns ()); - if (! jv) - return; - - if (! indexed_assign_conforms (1, jv.capacity (), rhs_nr, rhs_nc)) - { - ::error ("A(int,matrix) = X: X must be a row vector with the same"); - ::error ("number of elements as matrix"); - return; - } - maybe_resize (i, jv.max ()); - if (error_state) - return; - - do_matrix_assignment (rhs, i, jv); - } - break; - case string_constant: - gripe_string_invalid (); - break; - case range_constant: - { - Range rj = tmp_j.range_value (); - if (! indexed_assign_conforms (1, rj.nelem (), rhs_nr, rhs_nc)) - { - ::error ("A(int,range) = X: X must be a row vector with the same"); - ::error ("number of elements as range"); - return; - } - - int nc = columns (); - if (nc == 2 && is_zero_one (rj) && rhs_nc == 1) - { - do_matrix_assignment (rhs, i, 1); - } - else if (nc == 2 && is_one_zero (rj) && rhs_nc == 1) - { - do_matrix_assignment (rhs, i, 0); - } + double d = tmp.double_value (); + if (d > 0.0) + retval[0] = tree_constant (log (d)); else { - if (index_check (rj, "column") < 0) - return; - maybe_resize (i, tree_to_mat_idx (rj.max ())); - if (error_state) - return; - - do_matrix_assignment (rhs, i, rj); + Complex dtmp (d); + retval[0] = tree_constant (log (dtmp)); } } break; - case magic_colon: + case tree_constant_rep::complex_scalar_constant: { - int nc = columns (); - int nr = rows (); - if (nc == 0 && nr == 0 && rhs_nr == 1) - { - if (rhs.is_complex_type ()) - { - complex_matrix = new ComplexMatrix (); - type_tag = complex_matrix_constant; - } - else - { - matrix = new Matrix (); - type_tag = matrix_constant; - } - maybe_resize (i, rhs_nc-1); - if (error_state) - return; - } - else if (indexed_assign_conforms (1, nc, rhs_nr, rhs_nc)) - { - maybe_resize (i, nc-1); - if (error_state) - return; - } - else if (rhs_nr == 0 && rhs_nc == 0) - { - if (i < 0 || i >= nr) - { - ::error ("A(int,:) = []: row index out of range"); - return; - } - } - else - { - ::error ("A(int,:) = X: X must be a row vector with the same"); - ::error ("number of columns as A"); - return; - } - - do_matrix_assignment (rhs, i, magic_colon); + Complex c = tmp.complex_value (); + retval[0] = tree_constant (log (c)); } break; default: - panic_impossible (); break; } + return retval; } -/* MA2 */ -void -tree_constant_rep::do_matrix_assignment (tree_constant& rhs, idx_vector& iv, - tree_constant& j_arg) +tree_constant * +matrix_sqrt (const tree_constant& a) { - tree_constant tmp_j = j_arg.make_numeric_or_range_or_magic (); + tree_constant *retval = new tree_constant [2]; - tree_constant_rep::constant_type jtype = tmp_j.const_type (); - - int rhs_nr = rhs.rows (); - int rhs_nc = rhs.columns (); - - switch (jtype) + tree_constant tmp = a.make_numeric ();; + + if (tmp.rows () == 0 || tmp.columns () == 0) { - case complex_scalar_constant: - case scalar_constant: - { - int j = tree_to_mat_idx (tmp_j.double_value ()); - if (index_check (j, "column") < 0) - return; - if (! indexed_assign_conforms (iv.capacity (), 1, rhs_nr, rhs_nc)) - { - ::error ("A(matrix,int) = X: X must be a column vector with the"); - ::error ("same number of elements as matrix"); - return; - } - maybe_resize (iv.max (), j); - if (error_state) - return; + int flag = user_pref.propagate_empty_matrices; + if (flag != 0) + { + if (flag < 0) + gripe_empty_arg ("sqrtm", 0); + Matrix m; + retval = new tree_constant [2]; + retval[0] = tree_constant (m); + return retval; + } + else + gripe_empty_arg ("sqrtm", 1); + } - do_matrix_assignment (rhs, iv, j); - } - break; - case complex_matrix_constant: - case matrix_constant: + switch (tmp.const_type ()) + { + case tree_constant_rep::matrix_constant: { - Matrix mj = tmp_j.matrix_value (); - idx_vector jv (mj, user_pref.do_fortran_indexing, "column", - columns ()); - if (! jv) - return; - - if (! indexed_assign_conforms (iv.capacity (), jv.capacity (), - rhs_nr, rhs_nc)) - { - ::error ("A(r_mat,c_mat) = X: the number of rows in X must match"); - ::error ("the number of elements in r_mat and the number of"); - ::error ("columns in X must match the number of elements in c_mat"); - return; - } - maybe_resize (iv.max (), jv.max ()); - if (error_state) - return; + Matrix m = tmp.matrix_value (); - do_matrix_assignment (rhs, iv, jv); - } - break; - case string_constant: - gripe_string_invalid (); - break; - case range_constant: - { - Range rj = tmp_j.range_value (); - if (! indexed_assign_conforms (iv.capacity (), rj.nelem (), - rhs_nr, rhs_nc)) - { - ::error ("A(matrix,range) = X: the number of rows in X must match"); - ::error ("the number of elements in matrix and the number of"); - ::error ("columns in X must match the number of elements in range"); - return; - } + int nr = m.rows (); + int nc = m.columns (); - int nc = columns (); - if (nc == 2 && is_zero_one (rj) && rhs_nc == 1) - { - do_matrix_assignment (rhs, iv, 1); - } - else if (nc == 2 && is_one_zero (rj) && rhs_nc == 1) - { - do_matrix_assignment (rhs, iv, 0); - } + if (nr == 0 || nc == 0 || nr != nc) + gripe_square_matrix_required ("sqrtm"); else { - if (index_check (rj, "column") < 0) - return; - maybe_resize (iv.max (), tree_to_mat_idx (rj.max ())); - if (error_state) - return; + EIG m_eig (m); + ComplexColumnVector lambda (m_eig.eigenvalues ()); + ComplexMatrix Q (m_eig.eigenvectors ()); - do_matrix_assignment (rhs, iv, rj); + for (int i = 0; i < nr; i++) + { + Complex elt = lambda.elem (i); + if (imag (elt) == 0.0 && real (elt) > 0.0) + lambda.elem (i) = sqrt (real (elt)); + else + lambda.elem (i) = sqrt (elt); + } + + ComplexDiagMatrix D (lambda); + ComplexMatrix result = Q * D * Q.inverse (); + + retval[0] = tree_constant (result); } } break; - case magic_colon: + case tree_constant_rep::complex_matrix_constant: { - int nc = columns (); - int new_nc = nc; - if (nc == 0) - new_nc = rhs_nc; + ComplexMatrix m = tmp.complex_matrix_value (); - if (indexed_assign_conforms (iv.capacity (), new_nc, - rhs_nr, rhs_nc)) + int nr = m.rows (); + int nc = m.columns (); + + if (nr == 0 || nc == 0 || nr != nc) + gripe_square_matrix_required ("sqrtm"); + else { - maybe_resize (iv.max (), new_nc-1); - if (error_state) - return; + EIG m_eig (m); + ComplexColumnVector lambda (m_eig.eigenvalues ()); + ComplexMatrix Q (m_eig.eigenvectors ()); + + for (int i = 0; i < nr; i++) + { + Complex elt = lambda.elem (i); + if (imag (elt) == 0.0 && real (elt) > 0.0) + lambda.elem (i) = sqrt (real (elt)); + else + lambda.elem (i) = sqrt (elt); + } + + ComplexDiagMatrix D (lambda); + ComplexMatrix result = Q * D * Q.inverse (); + + retval[0] = tree_constant (result); } - else if (rhs_nr == 0 && rhs_nc == 0) - { - if (iv.max () >= rows ()) - { - ::error ("A(matrix,:) = []: row index out of range"); - return; - } - } + } + break; + case tree_constant_rep::scalar_constant: + { + double d = tmp.double_value (); + if (d > 0.0) + retval[0] = tree_constant (sqrt (d)); else { - ::error ("A(matrix,:) = X: the number of rows in X must match the"); - ::error ("number of elements in matrix, and the number of columns"); - ::error ("in X must match the number of columns in A"); - return; + Complex dtmp (d); + retval[0] = tree_constant (sqrt (dtmp)); } + } + break; + case tree_constant_rep::complex_scalar_constant: + { + Complex c = tmp.complex_value (); + retval[0] = tree_constant (log (c)); + } + break; + default: + break; + } + return retval; +} - do_matrix_assignment (rhs, iv, magic_colon); - } +tree_constant * +column_max (const tree_constant *args, int nargin, int nargout) +{ + tree_constant *retval = NULL_TREE_CONST; + + tree_constant arg1; + tree_constant arg2; + tree_constant_rep::constant_type arg1_type = + tree_constant_rep::unknown_constant; + tree_constant_rep::constant_type arg2_type = + tree_constant_rep::unknown_constant; + + switch (nargin) + { + case 3: + arg2 = args[2].make_numeric (); + arg2_type = arg2.const_type (); +// Fall through... + case 2: + arg1 = args[1].make_numeric (); + arg1_type = arg1.const_type (); break; default: panic_impossible (); break; } -} -/* MA3 */ -void -tree_constant_rep::do_matrix_assignment (tree_constant& rhs, - Range& ri, tree_constant& j_arg) -{ - tree_constant tmp_j = j_arg.make_numeric_or_range_or_magic (); - - tree_constant_rep::constant_type jtype = tmp_j.const_type (); - - int rhs_nr = rhs.rows (); - int rhs_nc = rhs.columns (); - - switch (jtype) + if (nargin == 2 && nargout == 1) { - case complex_scalar_constant: - case scalar_constant: - { - int j = tree_to_mat_idx (tmp_j.double_value ()); - if (index_check (j, "column") < 0) - return; - if (! indexed_assign_conforms (ri.nelem (), 1, rhs_nr, rhs_nc)) + retval = new tree_constant [2]; + switch (arg1_type) + { + case tree_constant_rep::scalar_constant: + retval[0] = tree_constant (arg1.double_value ()); + break; + case tree_constant_rep::complex_scalar_constant: + retval[0] = tree_constant (arg1.complex_value ()); + break; + case tree_constant_rep::matrix_constant: + { + Matrix m = arg1.matrix_value (); + if (m.rows () == 1) + retval[0] = tree_constant (m.row_max ()); + else + retval[0] = tree_constant (m.column_max (), 0); + } + break; + case tree_constant_rep::complex_matrix_constant: + { + ComplexMatrix m = arg1.complex_matrix_value (); + if (m.rows () == 1) + retval[0] = tree_constant (m.row_max ()); + else + retval[0] = tree_constant (m.column_max (), 0); + } + break; + default: + panic_impossible (); + break; + } + } + else if (nargin == 2 && nargout == 2) + { + retval = new tree_constant [2]; + switch (arg1_type) + { + case tree_constant_rep::scalar_constant: { - ::error ("A(range,int) = X: X must be a column vector with the"); - ::error ("same number of elements as range"); - return; - } - maybe_resize (tree_to_mat_idx (ri.max ()), j); - if (error_state) - return; - - do_matrix_assignment (rhs, ri, j); - } - break; - case complex_matrix_constant: - case matrix_constant: - { - Matrix mj = tmp_j.matrix_value (); - idx_vector jv (mj, user_pref.do_fortran_indexing, "column", - columns ()); - if (! jv) - return; - - if (! indexed_assign_conforms (ri.nelem (), jv.capacity (), - rhs_nr, rhs_nc)) - { - ::error ("A(range,matrix) = X: the number of rows in X must match"); - ::error ("the number of elements in range and the number of"); - ::error ("columns in X must match the number of elements in matrix"); - return; + retval[0] = tree_constant (arg1.double_value ()); + retval[1] = tree_constant (1); } - maybe_resize (tree_to_mat_idx (ri.max ()), jv.max ()); - if (error_state) - return; - - do_matrix_assignment (rhs, ri, jv); - } - break; - case string_constant: - gripe_string_invalid (); - break; - case range_constant: - { - Range rj = tmp_j.range_value (); - if (! indexed_assign_conforms (ri.nelem (), rj.nelem (), - rhs_nr, rhs_nc)) + break; + case tree_constant_rep::complex_scalar_constant: { - ::error ("A(r_range,c_range) = X: the number of rows in X must"); - ::error ("match the number of elements in r_range and the number"); - ::error ("of columns in X must match the number of elements in"); - ::error ("c_range"); - return; + retval[0] = tree_constant (arg1.complex_value ()); + retval[1] = tree_constant (1); } - - int nc = columns (); - if (nc == 2 && is_zero_one (rj) && rhs_nc == 1) - { - do_matrix_assignment (rhs, ri, 1); - } - else if (nc == 2 && is_one_zero (rj) && rhs_nc == 1) - { - do_matrix_assignment (rhs, ri, 0); - } - else + break; + case tree_constant_rep::matrix_constant: { - if (index_check (rj, "column") < 0) - return; - - maybe_resize (tree_to_mat_idx (ri.max ()), - tree_to_mat_idx (rj.max ())); - - if (error_state) - return; - - do_matrix_assignment (rhs, ri, rj); - } - } - break; - case magic_colon: - { - int nc = columns (); - int new_nc = nc; - if (nc == 0) - new_nc = rhs_nc; - - if (indexed_assign_conforms (ri.nelem (), new_nc, rhs_nr, rhs_nc)) - { - maybe_resize (tree_to_mat_idx (ri.max ()), new_nc-1); - if (error_state) - return; - } - else if (rhs_nr == 0 && rhs_nc == 0) - { - int b = tree_to_mat_idx (ri.min ()); - int l = tree_to_mat_idx (ri.max ()); - if (b < 0 || l >= rows ()) + Matrix m = arg1.matrix_value (); + if (m.rows () == 1) { - ::error ("A(range,:) = []: row index out of range"); - return; + retval[0] = tree_constant (m.row_max ()); + retval[1] = tree_constant (m.row_max_loc ()); + } + else + { + retval[0] = tree_constant (m.column_max (), 0); + retval[1] = tree_constant (m.column_max_loc (), 0); } } - else + break; + case tree_constant_rep::complex_matrix_constant: { - ::error ("A(range,:) = X: the number of rows in X must match the"); - ::error ("number of elements in range, and the number of columns"); - ::error ("in X must match the number of columns in A"); - return; - } - - do_matrix_assignment (rhs, ri, magic_colon); - } - break; - default: - panic_impossible (); - break; - } -} - -/* MA4 */ -void -tree_constant_rep::do_matrix_assignment (tree_constant& rhs, - tree_constant_rep::constant_type i, - tree_constant& j_arg) -{ - tree_constant tmp_j = j_arg.make_numeric_or_range_or_magic (); - - tree_constant_rep::constant_type jtype = tmp_j.const_type (); - - int rhs_nr = rhs.rows (); - int rhs_nc = rhs.columns (); - - switch (jtype) - { - case complex_scalar_constant: - case scalar_constant: - { - int j = tree_to_mat_idx (tmp_j.double_value ()); - if (index_check (j, "column") < 0) - return; - int nr = rows (); - int nc = columns (); - if (nr == 0 && nc == 0 && rhs_nc == 1) - { - if (rhs.is_complex_type ()) + ComplexMatrix m = arg1.complex_matrix_value (); + if (m.rows () == 1) { - complex_matrix = new ComplexMatrix (); - type_tag = complex_matrix_constant; + retval[0] = tree_constant (m.row_max ()); + retval[1] = tree_constant (m.row_max_loc ()); } else { - matrix = new Matrix (); - type_tag = matrix_constant; - } - maybe_resize (rhs_nr-1, j); - if (error_state) - return; - } - else if (indexed_assign_conforms (nr, 1, rhs_nr, rhs_nc)) - { - maybe_resize (nr-1, j); - if (error_state) - return; - } - else if (rhs_nr == 0 && rhs_nc == 0) - { - if (j < 0 || j >= nc) - { - ::error ("A(:,int) = []: column index out of range"); - return; - } - } - else - { - ::error ("A(:,int) = X: X must be a column vector with the same"); - ::error ("number of rows as A"); - return; - } - - do_matrix_assignment (rhs, magic_colon, j); - } - break; - case complex_matrix_constant: - case matrix_constant: - { - Matrix mj = tmp_j.matrix_value (); - idx_vector jv (mj, user_pref.do_fortran_indexing, "column", - columns ()); - if (! jv) - return; - - int nr = rows (); - int new_nr = nr; - if (nr == 0) - new_nr = rhs_nr; - - if (indexed_assign_conforms (new_nr, jv.capacity (), - rhs_nr, rhs_nc)) - { - maybe_resize (new_nr-1, jv.max ()); - if (error_state) - return; - } - else if (rhs_nr == 0 && rhs_nc == 0) - { - if (jv.max () >= columns ()) - { - ::error ("A(:,matrix) = []: column index out of range"); - return; + retval[0] = tree_constant (m.column_max (), 0); + retval[1] = tree_constant (m.column_max_loc (), 0); } } - else - { - ::error ("A(:,matrix) = X: the number of rows in X must match the"); - ::error ("number of rows in A, and the number of columns in X must"); - ::error ("match the number of elements in matrix"); - return; - } - - do_matrix_assignment (rhs, magic_colon, jv); - } - break; - case string_constant: - gripe_string_invalid (); - break; - case range_constant: - { - Range rj = tmp_j.range_value (); - int nr = rows (); - int new_nr = nr; - if (nr == 0) - new_nr = rhs_nr; - - if (indexed_assign_conforms (new_nr, rj.nelem (), rhs_nr, rhs_nc)) - { - int nc = columns (); - if (nc == 2 && is_zero_one (rj) && rhs_nc == 1) + break; + default: + panic_impossible (); + break; + } + } + else if (nargin == 3) + { + if (arg1.rows () == arg2.rows () + && arg1.columns () == arg2.columns ()) + { + retval = new tree_constant [2]; + switch (arg1_type) + { + case tree_constant_rep::scalar_constant: { - do_matrix_assignment (rhs, magic_colon, 1); + double result; + double a_elem = arg1.double_value (); + double b_elem = arg2.double_value (); + result = MAX (a_elem, b_elem); + retval[0] = tree_constant (result); } - else if (nc == 2 && is_one_zero (rj) && rhs_nc == 1) + break; + case tree_constant_rep::complex_scalar_constant: { - do_matrix_assignment (rhs, magic_colon, 0); + Complex result; + Complex a_elem = arg1.complex_value (); + Complex b_elem = arg2.complex_value (); + if (abs (a_elem) > abs (b_elem)) + result = a_elem; + else + result = b_elem; + retval[0] = tree_constant (result); + } + break; + case tree_constant_rep::matrix_constant: + { + Matrix result; + result = max (arg1.matrix_value (), arg2.matrix_value ()); + retval[0] = tree_constant (result); } - else - { - if (index_check (rj, "column") < 0) - return; - maybe_resize (new_nr-1, tree_to_mat_idx (rj.max ())); - if (error_state) - return; - } - } - else if (rhs_nr == 0 && rhs_nc == 0) - { - int b = tree_to_mat_idx (rj.min ()); - int l = tree_to_mat_idx (rj.max ()); - if (b < 0 || l >= columns ()) + break; + case tree_constant_rep::complex_matrix_constant: { - ::error ("A(:,range) = []: column index out of range"); - return; + ComplexMatrix result; + result = max (arg1.complex_matrix_value (), + arg2.complex_matrix_value ()); + retval[0] = tree_constant (result); } - } - else - { - ::error ("A(:,range) = X: the number of rows in X must match the"); - ::error ("number of rows in A, and the number of columns in X"); - ::error ("must match the number of elements in range"); - return; - } + break; + default: + panic_impossible (); + break; + } + } + else + error ("max: nonconformant matrices"); + } + else + panic_impossible (); + + return retval; +} - do_matrix_assignment (rhs, magic_colon, rj); - } - break; - case magic_colon: -// a(:,:) = foo is equivalent to a = foo. - do_matrix_assignment (rhs, magic_colon, magic_colon); +tree_constant * +column_min (const tree_constant *args, int nargin, int nargout) +{ + tree_constant *retval = NULL_TREE_CONST; + + tree_constant arg1; + tree_constant arg2; + tree_constant_rep::constant_type arg1_type = + tree_constant_rep::unknown_constant; + tree_constant_rep::constant_type arg2_type = + tree_constant_rep::unknown_constant; + + switch (nargin) + { + case 3: + arg2 = args[2].make_numeric (); + arg2_type = arg2.const_type (); +// Fall through... + case 2: + arg1 = args[1].make_numeric (); + arg1_type = arg1.const_type (); break; default: panic_impossible (); break; } -} -/* - * Functions that actually handle assignment to a matrix using two - * index values. - * - * idx2 - * +---+---+----+----+ - * idx1 | i | v | r | c | - * ---------+---+---+----+----+ - * integer | 1 | 5 | 9 | 13 | - * ---------+---+---+----+----+ - * vector | 2 | 6 | 10 | 14 | - * ---------+---+---+----+----+ - * range | 3 | 7 | 11 | 15 | - * ---------+---+---+----+----+ - * colon | 4 | 8 | 12 | 16 | - * ---------+---+---+----+----+ - */ - -/* 1 */ -void -tree_constant_rep::do_matrix_assignment (tree_constant& rhs, int i, int j) -{ - REP_ELEM_ASSIGN (i, j, rhs.double_value (), rhs.complex_value (), - rhs.is_real_type ()); -} - -/* 2 */ -void -tree_constant_rep::do_matrix_assignment (tree_constant& rhs, int i, - idx_vector& jv) -{ - REP_RHS_MATRIX (rhs, rhs_m, rhs_cm, rhs_nr, rhs_nc); - - for (int j = 0; j < jv.capacity (); j++) - REP_ELEM_ASSIGN (i, jv.elem (j), rhs_m.elem (0, j), - rhs_cm.elem (0, j), rhs.is_real_type ()); -} - -/* 3 */ -void -tree_constant_rep::do_matrix_assignment (tree_constant& rhs, int i, Range& rj) -{ - REP_RHS_MATRIX (rhs, rhs_m, rhs_cm, rhs_nr, rhs_nc); - - double b = rj.base (); - double increment = rj.inc (); - - for (int j = 0; j < rj.nelem (); j++) - { - double tmp = b + j * increment; - int col = tree_to_mat_idx (tmp); - REP_ELEM_ASSIGN (i, col, rhs_m.elem (0, j), rhs_cm.elem (0, j), - rhs.is_real_type ()); - } -} - -/* 4 */ -void -tree_constant_rep::do_matrix_assignment (tree_constant& rhs, int i, - tree_constant_rep::constant_type mcj) -{ - assert (mcj == magic_colon); - - int nc = columns (); - - if (rhs.is_zero_by_zero ()) - { - delete_row (i); - } - else if (rhs.is_matrix_type ()) - { - REP_RHS_MATRIX (rhs, rhs_m, rhs_cm, rhs_nr, rhs_nc); - - for (int j = 0; j < nc; j++) - REP_ELEM_ASSIGN (i, j, rhs_m.elem (0, j), rhs_cm.elem (0, j), - rhs.is_real_type ()); - } - else if (rhs.is_scalar_type () && nc == 1) - { - REP_ELEM_ASSIGN (i, 0, rhs.double_value (), - rhs.complex_value (), rhs.is_real_type ()); - } - else - panic_impossible (); -} - -/* 5 */ -void -tree_constant_rep::do_matrix_assignment (tree_constant& rhs, - idx_vector& iv, int j) -{ - REP_RHS_MATRIX (rhs, rhs_m, rhs_cm, rhs_nr, rhs_nc); - - for (int i = 0; i < iv.capacity (); i++) + if (nargin == 2 && nargout == 1) { - int row = iv.elem (i); - REP_ELEM_ASSIGN (row, j, rhs_m.elem (i, 0), - rhs_cm.elem (i, 0), rhs.is_real_type ()); - } -} - -/* 6 */ -void -tree_constant_rep::do_matrix_assignment (tree_constant& rhs, - idx_vector& iv, idx_vector& jv) -{ - REP_RHS_MATRIX (rhs, rhs_m, rhs_cm, rhs_nr, rhs_nc); - - for (int i = 0; i < iv.capacity (); i++) - { - int row = iv.elem (i); - for (int j = 0; j < jv.capacity (); j++) - { - int col = jv.elem (j); - REP_ELEM_ASSIGN (row, col, rhs_m.elem (i, j), - rhs_cm.elem (i, j), rhs.is_real_type ()); - } - } -} - -/* 7 */ -void -tree_constant_rep::do_matrix_assignment (tree_constant& rhs, - idx_vector& iv, Range& rj) -{ - REP_RHS_MATRIX (rhs, rhs_m, rhs_cm, rhs_nr, rhs_nc); - - double b = rj.base (); - double increment = rj.inc (); - - for (int i = 0; i < iv.capacity (); i++) - { - int row = iv.elem (i); - for (int j = 0; j < rj.nelem (); j++) + retval = new tree_constant [2]; + switch (arg1_type) { - double tmp = b + j * increment; - int col = tree_to_mat_idx (tmp); - REP_ELEM_ASSIGN (row, col, rhs_m.elem (i, j), - rhs_cm.elem (i, j), rhs.is_real_type ()); - } - } -} - -/* 8 */ -void -tree_constant_rep::do_matrix_assignment (tree_constant& rhs, idx_vector& iv, - tree_constant_rep::constant_type mcj) -{ - assert (mcj == magic_colon); - - if (rhs.is_zero_by_zero ()) - { - delete_rows (iv); - } - else - { - REP_RHS_MATRIX (rhs, rhs_m, rhs_cm, rhs_nr, rhs_nc); - - int nc = columns (); - - for (int j = 0; j < nc; j++) - { - for (int i = 0; i < iv.capacity (); i++) - { - int row = iv.elem (i); - REP_ELEM_ASSIGN (row, j, rhs_m.elem (i, j), - rhs_cm.elem (i, j), rhs.is_real_type ()); - } - } - } -} - -/* 9 */ -void -tree_constant_rep::do_matrix_assignment (tree_constant& rhs, Range& ri, int j) -{ - REP_RHS_MATRIX (rhs, rhs_m, rhs_cm, rhs_nr, rhs_nc); - - double b = ri.base (); - double increment = ri.inc (); - - for (int i = 0; i < ri.nelem (); i++) - { - double tmp = b + i * increment; - int row = tree_to_mat_idx (tmp); - REP_ELEM_ASSIGN (row, j, rhs_m.elem (i, 0), - rhs_cm.elem (i, 0), rhs.is_real_type ()); - } -} - -/* 10 */ -void -tree_constant_rep::do_matrix_assignment (tree_constant& rhs, Range& ri, - idx_vector& jv) -{ - REP_RHS_MATRIX (rhs, rhs_m, rhs_cm, rhs_nr, rhs_nc); - - double b = ri.base (); - double increment = ri.inc (); - - for (int j = 0; j < jv.capacity (); j++) - { - int col = jv.elem (j); - for (int i = 0; i < ri.nelem (); i++) - { - double tmp = b + i * increment; - int row = tree_to_mat_idx (tmp); - REP_ELEM_ASSIGN (row, col, rhs_m.elem (i, j), - rhs_m.elem (i, j), rhs.is_real_type ()); + case tree_constant_rep::scalar_constant: + retval[0] = tree_constant (arg1.double_value ()); + break; + case tree_constant_rep::complex_scalar_constant: + retval[0] = tree_constant (arg1.complex_value ()); + break; + case tree_constant_rep::matrix_constant: + { + Matrix m = arg1.matrix_value (); + if (m.rows () == 1) + retval[0] = tree_constant (m.row_min ()); + else + retval[0] = tree_constant (m.column_min (), 0); + } + break; + case tree_constant_rep::complex_matrix_constant: + { + ComplexMatrix m = arg1.complex_matrix_value (); + if (m.rows () == 1) + retval[0] = tree_constant (m.row_min ()); + else + retval[0] = tree_constant (m.column_min (), 0); + } + break; + default: + panic_impossible (); + break; } } -} - -/* 11 */ -void -tree_constant_rep::do_matrix_assignment (tree_constant& rhs, Range& ri, - Range& rj) -{ - double ib = ri.base (); - double iinc = ri.inc (); - double jb = rj.base (); - double jinc = rj.inc (); - - REP_RHS_MATRIX (rhs, rhs_m, rhs_cm, rhs_nr, rhs_nc); - - for (int i = 0; i < ri.nelem (); i++) + else if (nargin == 2 && nargout == 2) { - double itmp = ib + i * iinc; - int row = tree_to_mat_idx (itmp); - for (int j = 0; j < rj.nelem (); j++) - { - double jtmp = jb + j * jinc; - int col = tree_to_mat_idx (jtmp); - REP_ELEM_ASSIGN (row, col, rhs_m.elem (i, j), - rhs_cm.elem (i, j), rhs.is_real_type ()); - } - } -} - -/* 12 */ -void -tree_constant_rep::do_matrix_assignment (tree_constant& rhs, Range& ri, - tree_constant_rep::constant_type mcj) -{ - assert (mcj == magic_colon); - - if (rhs.is_zero_by_zero ()) - { - delete_rows (ri); - } - else - { - REP_RHS_MATRIX (rhs, rhs_m, rhs_cm, rhs_nr, rhs_nc); - - double ib = ri.base (); - double iinc = ri.inc (); - - int nc = columns (); - - for (int i = 0; i < ri.nelem (); i++) - { - double itmp = ib + i * iinc; - int row = tree_to_mat_idx (itmp); - for (int j = 0; j < nc; j++) - REP_ELEM_ASSIGN (row, j, rhs_m.elem (i, j), - rhs_cm.elem (i, j), rhs.is_real_type ()); - } + retval = new tree_constant [2]; + switch (arg1_type) + { + case tree_constant_rep::scalar_constant: + { + retval[0] = tree_constant (arg1.double_value ()); + retval[1] = tree_constant (1); + } + break; + case tree_constant_rep::complex_scalar_constant: + { + retval[0] = tree_constant (arg1.complex_value ()); + retval[1] = tree_constant (1); + } + break; + case tree_constant_rep::matrix_constant: + { + Matrix m = arg1.matrix_value (); + if (m.rows () == 1) + { + retval[0] = tree_constant (m.row_min ()); + retval[1] = tree_constant (m.row_min_loc ()); + } + else + { + retval[0] = tree_constant (m.column_min (), 0); + retval[1] = tree_constant (m.column_min_loc (), 0); + } + } + break; + case tree_constant_rep::complex_matrix_constant: + { + ComplexMatrix m = arg1.complex_matrix_value (); + if (m.rows () == 1) + { + retval[0] = tree_constant (m.row_min ()); + retval[1] = tree_constant (m.row_min_loc ()); + } + else + { + retval[0] = tree_constant (m.column_min (), 0); + retval[1] = tree_constant (m.column_min_loc (), 0); + } + } + break; + default: + panic_impossible (); + break; + } } -} - -/* 13 */ -void -tree_constant_rep::do_matrix_assignment (tree_constant& rhs, - tree_constant_rep::constant_type mci, - int j) -{ - assert (mci == magic_colon); - - int nr = rows (); - - if (rhs.is_zero_by_zero ()) - { - delete_column (j); - } - else if (rhs.is_matrix_type ()) - { - REP_RHS_MATRIX (rhs, rhs_m, rhs_cm, rhs_nr, rhs_nc); - - for (int i = 0; i < nr; i++) - REP_ELEM_ASSIGN (i, j, rhs_m.elem (i, 0), - rhs_cm.elem (i, 0), rhs.is_real_type ()); - } - else if (rhs.is_scalar_type () && nr == 1) + else if (nargin == 3) { - REP_ELEM_ASSIGN (0, j, rhs.double_value (), - rhs.complex_value (), rhs.is_real_type ()); - } - else - panic_impossible (); -} - -/* 14 */ -void -tree_constant_rep::do_matrix_assignment (tree_constant& rhs, - tree_constant_rep::constant_type mci, - idx_vector& jv) -{ - assert (mci == magic_colon); - - if (rhs.is_zero_by_zero ()) - { - delete_columns (jv); - } - else - { - REP_RHS_MATRIX (rhs, rhs_m, rhs_cm, rhs_nr, rhs_nc); - - int nr = rows (); - - for (int i = 0; i < nr; i++) + if (arg1.rows () == arg2.rows () + && arg1.columns () == arg2.columns ()) { - for (int j = 0; j < jv.capacity (); j++) - { - int col = jv.elem (j); - REP_ELEM_ASSIGN (i, col, rhs_m.elem (i, j), - rhs_cm.elem (i, j), rhs.is_real_type ()); + retval = new tree_constant [2]; + switch (arg1_type) + { + case tree_constant_rep::scalar_constant: + { + double result; + double a_elem = arg1.double_value (); + double b_elem = arg2.double_value (); + result = MIN (a_elem, b_elem); + retval[0] = tree_constant (result); + } + break; + case tree_constant_rep::complex_scalar_constant: + { + Complex result; + Complex a_elem = arg1.complex_value (); + Complex b_elem = arg2.complex_value (); + if (abs (a_elem) < abs (b_elem)) + result = a_elem; + else + result = b_elem; + retval[0] = tree_constant (result); + } + break; + case tree_constant_rep::matrix_constant: + { + Matrix result; + result = min (arg1.matrix_value (), arg2.matrix_value ()); + retval[0] = tree_constant (result); + } + break; + case tree_constant_rep::complex_matrix_constant: + { + ComplexMatrix result; + result = min (arg1.complex_matrix_value (), + arg2.complex_matrix_value ()); + retval[0] = tree_constant (result); + } + break; + default: + panic_impossible (); + break; } } - } -} - -/* 15 */ -void -tree_constant_rep::do_matrix_assignment (tree_constant& rhs, - tree_constant_rep::constant_type mci, - Range& rj) -{ - assert (mci == magic_colon); - - if (rhs.is_zero_by_zero ()) - { - delete_columns (rj); - } - else - { - REP_RHS_MATRIX (rhs, rhs_m, rhs_cm, rhs_nr, rhs_nc); - - int nr = rows (); - - double jb = rj.base (); - double jinc = rj.inc (); - - for (int j = 0; j < rj.nelem (); j++) - { - double jtmp = jb + j * jinc; - int col = tree_to_mat_idx (jtmp); - for (int i = 0; i < nr; i++) - { - REP_ELEM_ASSIGN (i, col, rhs_m.elem (i, j), - rhs_cm.elem (i, j), rhs.is_real_type ()); - } - } - } -} - -/* 16 */ -void -tree_constant_rep::do_matrix_assignment (tree_constant& rhs, - tree_constant_rep::constant_type mci, - tree_constant_rep::constant_type mcj) -{ - assert (mci == magic_colon && mcj == magic_colon); - - switch (type_tag) - { - case scalar_constant: - break; - case matrix_constant: - delete matrix; - break; - case complex_scalar_constant: - delete complex_scalar; - break; - case complex_matrix_constant: - delete complex_matrix; - break; - case string_constant: - delete [] string; - break; - case range_constant: - delete range; - break; - case magic_colon: - default: - panic_impossible (); - break; - } - - type_tag = rhs.const_type (); - - switch (type_tag) - { - case scalar_constant: - scalar = rhs.double_value (); - break; - case matrix_constant: - matrix = new Matrix (rhs.matrix_value ()); - break; - case string_constant: - string = strsave (rhs.string_value ()); - break; - case complex_matrix_constant: - complex_matrix = new ComplexMatrix (rhs.complex_matrix_value ()); - break; - case complex_scalar_constant: - complex_scalar = new Complex (rhs.complex_value ()); - break; - case range_constant: - range = new Range (rhs.range_value ()); - break; - case magic_colon: - default: - panic_impossible (); - break; - } -} - -/* - * Functions for deleting rows or columns of a matrix. These are used - * to handle statements like - * - * M (i, j) = [] - */ -void -tree_constant_rep::delete_row (int idx) -{ - if (type_tag == matrix_constant) - { - int nr = matrix->rows (); - int nc = matrix->columns (); - Matrix *new_matrix = new Matrix (nr-1, nc); - int ii = 0; - for (int i = 0; i < nr; i++) - { - if (i != idx) - { - for (int j = 0; j < nc; j++) - new_matrix->elem (ii, j) = matrix->elem (i, j); - ii++; - } - } - delete matrix; - matrix = new_matrix; - } - else if (type_tag == complex_matrix_constant) - { - int nr = complex_matrix->rows (); - int nc = complex_matrix->columns (); - ComplexMatrix *new_matrix = new ComplexMatrix (nr-1, nc); - int ii = 0; - for (int i = 0; i < nr; i++) - { - if (i != idx) - { - for (int j = 0; j < nc; j++) - new_matrix->elem (ii, j) = complex_matrix->elem (i, j); - ii++; - } - } - delete complex_matrix; - complex_matrix = new_matrix; + else + error ("min: nonconformant matrices"); } else panic_impossible (); -} - -void -tree_constant_rep::delete_rows (idx_vector& iv) -{ - iv.sort_uniq (); - int num_to_delete = iv.length (); - - int nr = rows (); - int nc = columns (); - -// If deleting all rows of a column vector, make result 0x0. - if (nc == 1 && num_to_delete == nr) - nc = 0; - - if (type_tag == matrix_constant) - { - Matrix *new_matrix = new Matrix (nr-num_to_delete, nc); - if (nr > num_to_delete) - { - int ii = 0; - int idx = 0; - for (int i = 0; i < nr; i++) - { - if (i == iv.elem (idx)) - idx++; - else - { - for (int j = 0; j < nc; j++) - new_matrix->elem (ii, j) = matrix->elem (i, j); - ii++; - } - } - } - delete matrix; - matrix = new_matrix; - } - else if (type_tag == complex_matrix_constant) - { - ComplexMatrix *new_matrix = new ComplexMatrix (nr-num_to_delete, nc); - if (nr > num_to_delete) - { - int ii = 0; - int idx = 0; - for (int i = 0; i < nr; i++) - { - if (i == iv.elem (idx)) - idx++; - else - { - for (int j = 0; j < nc; j++) - new_matrix->elem (ii, j) = complex_matrix->elem (i, j); - ii++; - } - } - } - delete complex_matrix; - complex_matrix = new_matrix; - } - else - panic_impossible (); -} - -void -tree_constant_rep::delete_rows (Range& ri) -{ - ri.sort (); - int num_to_delete = ri.nelem (); - - int nr = rows (); - int nc = columns (); - -// If deleting all rows of a column vector, make result 0x0. - if (nc == 1 && num_to_delete == nr) - nc = 0; - - double ib = ri.base (); - double iinc = ri.inc (); - - int max_idx = tree_to_mat_idx (ri.max ()); - - if (type_tag == matrix_constant) - { - Matrix *new_matrix = new Matrix (nr-num_to_delete, nc); - if (nr > num_to_delete) - { - int ii = 0; - int idx = 0; - for (int i = 0; i < nr; i++) - { - double itmp = ib + idx * iinc; - int row = tree_to_mat_idx (itmp); - - if (i == row && row <= max_idx) - idx++; - else - { - for (int j = 0; j < nc; j++) - new_matrix->elem (ii, j) = matrix->elem (i, j); - ii++; - } - } - } - delete matrix; - matrix = new_matrix; - } - else if (type_tag == complex_matrix_constant) - { - ComplexMatrix *new_matrix = new ComplexMatrix (nr-num_to_delete, nc); - if (nr > num_to_delete) - { - int ii = 0; - int idx = 0; - for (int i = 0; i < nr; i++) - { - double itmp = ib + idx * iinc; - int row = tree_to_mat_idx (itmp); - - if (i == row && row <= max_idx) - idx++; - else - { - for (int j = 0; j < nc; j++) - new_matrix->elem (ii, j) = complex_matrix->elem (i, j); - ii++; - } - } - } - delete complex_matrix; - complex_matrix = new_matrix; - } - else - panic_impossible (); -} - -void -tree_constant_rep::delete_column (int idx) -{ - if (type_tag == matrix_constant) - { - int nr = matrix->rows (); - int nc = matrix->columns (); - Matrix *new_matrix = new Matrix (nr, nc-1); - int jj = 0; - for (int j = 0; j < nc; j++) - { - if (j != idx) - { - for (int i = 0; i < nr; i++) - new_matrix->elem (i, jj) = matrix->elem (i, j); - jj++; - } - } - delete matrix; - matrix = new_matrix; - } - else if (type_tag == complex_matrix_constant) - { - int nr = complex_matrix->rows (); - int nc = complex_matrix->columns (); - ComplexMatrix *new_matrix = new ComplexMatrix (nr, nc-1); - int jj = 0; - for (int j = 0; j < nc; j++) - { - if (j != idx) - { - for (int i = 0; i < nr; i++) - new_matrix->elem (i, jj) = complex_matrix->elem (i, j); - jj++; - } - } - delete complex_matrix; - complex_matrix = new_matrix; - } - else - panic_impossible (); -} - -void -tree_constant_rep::delete_columns (idx_vector& jv) -{ - jv.sort_uniq (); - int num_to_delete = jv.length (); - - int nr = rows (); - int nc = columns (); - -// If deleting all columns of a row vector, make result 0x0. - if (nr == 1 && num_to_delete == nc) - nr = 0; - - if (type_tag == matrix_constant) - { - Matrix *new_matrix = new Matrix (nr, nc-num_to_delete); - if (nc > num_to_delete) - { - int jj = 0; - int idx = 0; - for (int j = 0; j < nc; j++) - { - if (j == jv.elem (idx)) - idx++; - else - { - for (int i = 0; i < nr; i++) - new_matrix->elem (i, jj) = matrix->elem (i, j); - jj++; - } - } - } - delete matrix; - matrix = new_matrix; - } - else if (type_tag == complex_matrix_constant) - { - ComplexMatrix *new_matrix = new ComplexMatrix (nr, nc-num_to_delete); - if (nc > num_to_delete) - { - int jj = 0; - int idx = 0; - for (int j = 0; j < nc; j++) - { - if (j == jv.elem (idx)) - idx++; - else - { - for (int i = 0; i < nr; i++) - new_matrix->elem (i, jj) = complex_matrix->elem (i, j); - jj++; - } - } - } - delete complex_matrix; - complex_matrix = new_matrix; - } - else - panic_impossible (); -} - -void -tree_constant_rep::delete_columns (Range& rj) -{ - rj.sort (); - int num_to_delete = rj.nelem (); - - int nr = rows (); - int nc = columns (); - -// If deleting all columns of a row vector, make result 0x0. - if (nr == 1 && num_to_delete == nc) - nr = 0; - - double jb = rj.base (); - double jinc = rj.inc (); - - int max_idx = tree_to_mat_idx (rj.max ()); - - if (type_tag == matrix_constant) - { - Matrix *new_matrix = new Matrix (nr, nc-num_to_delete); - if (nc > num_to_delete) - { - int jj = 0; - int idx = 0; - for (int j = 0; j < nc; j++) - { - double jtmp = jb + idx * jinc; - int col = tree_to_mat_idx (jtmp); - - if (j == col && col <= max_idx) - idx++; - else - { - for (int i = 0; i < nr; i++) - new_matrix->elem (i, jj) = matrix->elem (i, j); - jj++; - } - } - } - delete matrix; - matrix = new_matrix; - } - else if (type_tag == complex_matrix_constant) - { - ComplexMatrix *new_matrix = new ComplexMatrix (nr, nc-num_to_delete); - if (nc > num_to_delete) - { - int jj = 0; - int idx = 0; - for (int j = 0; j < nc; j++) - { - double jtmp = jb + idx * jinc; - int col = tree_to_mat_idx (jtmp); - - if (j == col && col <= max_idx) - idx++; - else - { - for (int i = 0; i < nr; i++) - new_matrix->elem (i, jj) = complex_matrix->elem (i, j); - jj++; - } - } - } - delete complex_matrix; - complex_matrix = new_matrix; - } - else - panic_impossible (); -} - - -int -tree_constant_rep::valid_as_scalar_index (void) const -{ - int valid = type_tag == magic_colon - || (type_tag == scalar_constant && NINT (scalar) == 1) - || (type_tag == range_constant - && range->nelem () == 1 && NINT (range->base ()) == 1); - - return valid; -} - -tree_constant -tree_constant_rep::do_scalar_index (const tree_constant *args, - int nargs) const -{ - if (valid_scalar_indices (args, nargs)) - { - if (type_tag == scalar_constant) - return tree_constant (scalar); - else if (type_tag == complex_scalar_constant) - return tree_constant (*complex_scalar); - else - panic_impossible (); - } - else - { - int rows = 0; - int cols = 0; - - switch (nargs) - { - case 3: - { - if (args[2].is_matrix_type ()) - { - Matrix mj = args[2].matrix_value (); - - idx_vector j (mj, user_pref.do_fortran_indexing, ""); - if (! j) - return tree_constant (); - - int len = j.length (); - if (len == j.ones_count ()) - cols = len; - } - else if (args[2].const_type () == magic_colon - || (args[2].is_scalar_type () - && NINT (args[2].double_value ()) == 1)) - { - cols = 1; - } - else - break; - } -// Fall through... - case 2: - { - if (args[1].is_matrix_type ()) - { - Matrix mi = args[1].matrix_value (); - - idx_vector i (mi, user_pref.do_fortran_indexing, ""); - if (! i) - return tree_constant (); - - int len = i.length (); - if (len == i.ones_count ()) - rows = len; - } - else if (args[1].const_type () == magic_colon - || (args[1].is_scalar_type () - && NINT (args[1].double_value ()) == 1)) - { - rows = 1; - } - else if (args[1].is_scalar_type () - && NINT (args[1].double_value ()) == 0) - { - Matrix m (0, 0); - return tree_constant (m); - } - else - break; - - if (cols == 0) - { - if (user_pref.prefer_column_vectors) - cols = 1; - else - { - cols = rows; - rows = 1; - } - } - - if (type_tag == scalar_constant) - { - Matrix m (rows, cols, scalar); - return tree_constant (m); - } - else if (type_tag == complex_scalar_constant) - { - ComplexMatrix cm (rows, cols, *complex_scalar); - return tree_constant (cm); - } - else - panic_impossible (); - } - break; - default: - ::error ("illegal number of arguments for scalar type"); - return tree_constant (); - break; - } - } - - ::error ("index invalid or out of range for scalar type"); - return tree_constant (); -} - -tree_constant -tree_constant_rep::do_matrix_index (const tree_constant *args, - int nargin) const -{ - tree_constant retval; - - switch (nargin) - { - case 2: - if (args == NULL_TREE_CONST) - ::error ("matrix index is null"); - else if (args[1].is_undefined ()) - ::error ("matrix index is a null expression"); - else - retval = do_matrix_index (args[1]); - break; - case 3: - if (args == NULL_TREE_CONST) - ::error ("matrix indices are null"); - else if (args[1].is_undefined ()) - ::error ("first matrix index is a null expression"); - else if (args[2].is_undefined ()) - ::error ("second matrix index is a null expression"); - else - retval = do_matrix_index (args[1], args[2]); - break; - default: - ::error ("too many indices for matrix expression"); - break; - } - - return retval; -} - -tree_constant -tree_constant_rep::do_matrix_index (const tree_constant& i_arg) const -{ - tree_constant retval; - - int nr = rows (); - int nc = columns (); - - if (user_pref.do_fortran_indexing) - retval = fortran_style_matrix_index (i_arg); - else if (nr <= 1 || nc <= 1) - retval = do_vector_index (i_arg); - else - ::error ("single index only valid for row or column vector"); return retval; } -tree_constant -tree_constant_rep::fortran_style_matrix_index - (const tree_constant& i_arg) const +static void +mx_sort (Matrix& m, Matrix& idx, int return_idx) { - tree_constant retval; + int nr = m.rows (); + int nc = m.columns (); + idx.resize (nr, nc); + int i, j; + + if (return_idx) + { + for (j = 0; j < nc; j++) + for (i = 0; i < nr; i++) + idx.elem (i, j) = i+1; + } + + for (j = 0; j < nc; j++) + { + for (int gap = nr/2; gap > 0; gap /= 2) + for (i = gap; i < nr; i++) + for (int k = i - gap; + k >= 0 && m.elem (k, j) > m.elem (k+gap, j); + k -= gap) + { + double tmp = m.elem (k, j); + m.elem (k, j) = m.elem (k+gap, j); + m.elem (k+gap, j) = tmp; + + if (return_idx) + { + double tmp = idx.elem (k, j); + idx.elem (k, j) = idx.elem (k+gap, j); + idx.elem (k+gap, j) = tmp; + } + } + } +} - tree_constant tmp_i = i_arg.make_numeric_or_magic (); +static void +mx_sort (RowVector& v, RowVector& idx, int return_idx) +{ + int n = v.capacity (); + idx.resize (n); + int i; + + if (return_idx) + for (i = 0; i < n; i++) + idx.elem (i) = i+1; - tree_constant_rep::constant_type itype = tmp_i.const_type (); + for (int gap = n/2; gap > 0; gap /= 2) + for (i = gap; i < n; i++) + for (int k = i - gap; + k >= 0 && v.elem (k) > v.elem (k+gap); + k -= gap) + { + double tmp = v.elem (k); + v.elem (k) = v.elem (k+gap); + v.elem (k+gap) = tmp; - int nr = rows (); - int nc = columns (); + if (return_idx) + { + double tmp = idx.elem (k); + idx.elem (k) = idx.elem (k+gap); + idx.elem (k+gap) = tmp; + } + } +} + +static void +mx_sort (ComplexMatrix& cm, Matrix& idx, int return_idx) +{ + int nr = cm.rows (); + int nc = cm.columns (); + idx.resize (nr, nc); + int i, j; - switch (itype) + if (return_idx) + { + for (j = 0; j < nc; j++) + for (i = 0; i < nr; i++) + idx.elem (i, j) = i+1; + } + + for (j = 0; j < nc; j++) { - case complex_scalar_constant: - case scalar_constant: + for (int gap = nr/2; gap > 0; gap /= 2) + for (i = gap; i < nr; i++) + for (int k = i - gap; + k >= 0 && abs (cm.elem (k, j)) > abs (cm.elem (k+gap, j)); + k -= gap) + { + Complex ctmp = cm.elem (k, j); + cm.elem (k, j) = cm.elem (k+gap, j); + cm.elem (k+gap, j) = ctmp; + + if (return_idx) + { + double tmp = idx.elem (k, j); + idx.elem (k, j) = idx.elem (k+gap, j); + idx.elem (k+gap, j) = tmp; + } + } + } +} + +static void +mx_sort (ComplexRowVector& cv, RowVector& idx, int return_idx) +{ + int n = cv.capacity (); + idx.resize (n); + int i; + + if (return_idx) + for (i = 0; i < n; i++) + idx.elem (i) = i+1; + + for (int gap = n/2; gap > 0; gap /= 2) + for (i = gap; i < n; i++) + for (int k = i - gap; + k >= 0 && abs (cv.elem (k)) > abs (cv.elem (k+gap)); + k -= gap) + { + Complex tmp = cv.elem (k); + cv.elem (k) = cv.elem (k+gap); + cv.elem (k+gap) = tmp; + + if (return_idx) + { + double tmp = idx.elem (k); + idx.elem (k) = idx.elem (k+gap); + idx.elem (k+gap) = tmp; + } + } +} + +tree_constant * +sort (const tree_constant *args, int nargin, int nargout) +{ +// Assumes that we have been given the correct number of arguments. + + tree_constant *retval = NULL_TREE_CONST; + + int return_idx = nargout > 1; + if (return_idx) + retval = new tree_constant [3]; + else + retval = new tree_constant [2]; + + switch (args[1].const_type ()) + { + case tree_constant_rep::scalar_constant: { - int i = NINT (tmp_i.double_value ()); - int ii = fortran_row (i, nr) - 1; - int jj = fortran_column (i, nr) - 1; - if (index_check (i-1, "") < 0) - return tree_constant (); - if (range_max_check (i-1, nr * nc) < 0) - return tree_constant (); - retval = do_matrix_index (ii, jj); + retval [0] = tree_constant (args[1].double_value ()); + if (return_idx) + retval [1] = tree_constant (1.0); } break; - case complex_matrix_constant: - case matrix_constant: + case tree_constant_rep::complex_scalar_constant: + { + retval [0] = tree_constant (args[1].complex_value ()); + if (return_idx) + retval [1] = tree_constant (1.0); + } + break; + case tree_constant_rep::string_constant: + case tree_constant_rep::range_constant: + case tree_constant_rep::matrix_constant: { - Matrix mi = tmp_i.matrix_value (); - if (mi.rows () == 0 || mi.columns () == 0) + Matrix m = args[1].to_matrix (); + if (m.rows () == 1) { - Matrix mtmp; - retval = tree_constant (mtmp); + int nc = m.columns (); + RowVector v (nc); + for (int i = 0; i < nc; i++) + v.elem (i) = m.elem (0, i); + RowVector idx; + mx_sort (v, idx, return_idx); + + retval [0] = tree_constant (v, 0); + if (return_idx) + retval [1] = tree_constant (idx, 0); } else { -// Yes, we really do want to call this with mi. - retval = fortran_style_matrix_index (mi); +// Sorts m in place, optionally computes index Matrix. + Matrix idx; + mx_sort (m, idx, return_idx); + + retval [0] = tree_constant (m); + if (return_idx) + retval [1] = tree_constant (idx); } } break; - case string_constant: - gripe_string_invalid (); - break; - case range_constant: - gripe_range_invalid (); - break; - case magic_colon: - retval = do_matrix_index (magic_colon); + case tree_constant_rep::complex_matrix_constant: + { + ComplexMatrix cm = args[1].complex_matrix_value (); + if (cm.rows () == 1) + { + int nc = cm.columns (); + ComplexRowVector cv (nc); + for (int i = 0; i < nc; i++) + cv.elem (i) = cm.elem (0, i); + RowVector idx; + mx_sort (cv, idx, return_idx); + + retval [0] = tree_constant (cv, 0); + if (return_idx) + retval [1] = tree_constant (idx, 0); + } + else + { +// Sorts cm in place, optionally computes index Matrix. + Matrix idx; + mx_sort (cm, idx, return_idx); + + retval [0] = tree_constant (cm); + if (return_idx) + retval [1] = tree_constant (idx); + } + } break; default: panic_impossible (); @@ -5135,1033 +1275,174 @@ return retval; } -tree_constant -tree_constant_rep::fortran_style_matrix_index (const Matrix& mi) const +tree_constant * +feval (const tree_constant *args, int nargin, int nargout) { - assert (is_matrix_type ()); - - tree_constant retval; - - int nr = rows (); - int nc = columns (); - - int len = nr * nc; - - int index_nr = mi.rows (); - int index_nc = mi.columns (); - - if (index_nr >= 1 && index_nc >= 1) - { - const double *cop_out = (const double *) NULL; - const Complex *c_cop_out = (const Complex *) NULL; - int real_type = type_tag == matrix_constant; - if (real_type) - cop_out = matrix->data (); - else - c_cop_out = complex_matrix->data (); - - const double *cop_out_index = mi.data (); - - idx_vector iv (mi, 1, "", len); - if (! iv) - return tree_constant (); - - int result_size = iv.length (); - - if (nc == 1 || (nr != 1 && iv.one_zero_only ())) - { - CRMATRIX (m, cm, result_size, 1); - - for (int i = 0; i < result_size; i++) - { - int idx = iv.elem (i); - CRMATRIX_ASSIGN_ELEM (m, cm, i, 0, cop_out [idx], - c_cop_out [idx], real_type); - } - - ASSIGN_CRMATRIX_TO (retval, m, cm); - } - else if (nr == 1) - { - CRMATRIX (m, cm, 1, result_size); +// Assumes that we have been given the correct number of arguments. - for (int i = 0; i < result_size; i++) - { - int idx = iv.elem (i); - CRMATRIX_ASSIGN_ELEM (m, cm, 0, i, cop_out [idx], - c_cop_out [idx], real_type); - } - - ASSIGN_CRMATRIX_TO (retval, m, cm); - } - else - { - CRMATRIX (m, cm, index_nr, index_nc); - - for (int j = 0; j < index_nc; j++) - for (int i = 0; i < index_nr; i++) - { - double tmp = *cop_out_index++; - int idx = tree_to_mat_idx (tmp); - CRMATRIX_ASSIGN_ELEM (m, cm, i, j, cop_out [idx], - c_cop_out [idx], real_type); - } + tree_constant *retval = NULL_TREE_CONST; - ASSIGN_CRMATRIX_TO (retval, m, cm); - } - } - else - { - if (index_nr == 0 || index_nc == 0) - ::error ("empty matrix invalid as index"); - else - ::error ("invalid matrix index"); - return tree_constant (); - } - - return retval; -} - -tree_constant -tree_constant_rep::do_vector_index (const tree_constant& i_arg) const -{ - tree_constant retval; - - tree_constant tmp_i = i_arg.make_numeric_or_range_or_magic (); - - tree_constant_rep::constant_type itype = tmp_i.const_type (); - - int nr = rows (); - int nc = columns (); - - int len = nr > nc ? nr : nc; - - if (nr == 0 || nc == 0) + tree_fvc *fcn = is_valid_function (args[1], "feval", 1); + if (fcn != (tree_fvc *) NULL) { - ::error ("attempt to index empty matrix"); - return retval; - } - - assert ((nr == 1 || nc == 1) && ! user_pref.do_fortran_indexing); - - int swap_indices = (nr == 1); - - switch (itype) - { - case complex_scalar_constant: - case scalar_constant: - { - int i = tree_to_mat_idx (tmp_i.double_value ()); - if (index_check (i, "") < 0) - return tree_constant (); - if (swap_indices) - { - if (range_max_check (i, nc) < 0) - return tree_constant (); - retval = do_matrix_index (0, i); - } - else - { - if (range_max_check (i, nr) < 0) - return tree_constant (); - retval = do_matrix_index (i, 0); - } - } - break; - case complex_matrix_constant: - case matrix_constant: - { - Matrix mi = tmp_i.matrix_value (); - if (mi.rows () == 0 || mi.columns () == 0) - { - Matrix mtmp; - retval = tree_constant (mtmp); - } - else - { - idx_vector iv (mi, user_pref.do_fortran_indexing, "", len); - if (! iv) - return tree_constant (); - - if (swap_indices) - { - if (range_max_check (iv.max (), nc) < 0) - return tree_constant (); - retval = do_matrix_index (0, iv); - } - else - { - if (range_max_check (iv.max (), nr) < 0) - return tree_constant (); - retval = do_matrix_index (iv, 0); - } - } - } - break; - case string_constant: - gripe_string_invalid (); - break; - case range_constant: - { - Range ri = tmp_i.range_value (); - if (len == 2 && is_zero_one (ri)) - { - if (swap_indices) - retval = do_matrix_index (0, 1); - else - retval = do_matrix_index (1, 0); - } - else if (len == 2 && is_one_zero (ri)) - { - retval = do_matrix_index (0, 0); - } - else - { - if (index_check (ri, "") < 0) - return tree_constant (); - if (swap_indices) - { - if (range_max_check (tree_to_mat_idx (ri.max ()), nc) < 0) - return tree_constant (); - retval = do_matrix_index (0, ri); - } - else - { - if (range_max_check (tree_to_mat_idx (ri.max ()), nr) < 0) - return tree_constant (); - retval = do_matrix_index (ri, 0); - } - } - } - break; - case magic_colon: - if (swap_indices) - retval = do_matrix_index (0, magic_colon); - else - retval = do_matrix_index (magic_colon, 0); - break; - default: - panic_impossible (); - break; + args++; + nargin--; + retval = fcn->eval (0, nargout, args, nargin); } return retval; } tree_constant -tree_constant_rep::do_matrix_index (const tree_constant& i_arg, - const tree_constant& j_arg) const +eval_string (const char *string, int print, int ans_assign, + int& parse_status) { - tree_constant retval; - - tree_constant tmp_i = i_arg.make_numeric_or_range_or_magic (); - - tree_constant_rep::constant_type itype = tmp_i.const_type (); + begin_unwind_frame ("eval_string"); - switch (itype) - { - case complex_scalar_constant: - case scalar_constant: - { - int i = tree_to_mat_idx (tmp_i.double_value ()); - if (index_check (i, "row") < 0) - return tree_constant (); - retval = do_matrix_index (i, j_arg); - } - break; - case complex_matrix_constant: - case matrix_constant: - { - Matrix mi = tmp_i.matrix_value (); - idx_vector iv (mi, user_pref.do_fortran_indexing, "row", rows ()); - if (! iv) - return tree_constant (); + unwind_protect_int (get_input_from_eval_string); + unwind_protect_ptr (global_command); + unwind_protect_ptr (current_eval_string); + + get_input_from_eval_string = 1; + current_eval_string = string; + + YY_BUFFER_STATE old_buf = current_buffer (); + YY_BUFFER_STATE new_buf = create_buffer ((FILE *) NULL); - if (iv.length () == 0) - { - Matrix mtmp; - retval = tree_constant (mtmp); - } - else - retval = do_matrix_index (iv, j_arg); - } - break; - case string_constant: - gripe_string_invalid (); - break; - case range_constant: - { - Range ri = tmp_i.range_value (); - int nr = rows (); - if (nr == 2 && is_zero_one (ri)) - { - retval = do_matrix_index (1, j_arg); - } - else if (nr == 2 && is_one_zero (ri)) - { - retval = do_matrix_index (0, j_arg); - } - else - { - if (index_check (ri, "row") < 0) - return tree_constant (); - retval = do_matrix_index (ri, j_arg); - } - } - break; - case magic_colon: - retval = do_matrix_index (magic_colon, j_arg); - break; - default: - panic_impossible (); - break; - } + add_unwind_protect (restore_input_buffer, (void *) old_buf); + add_unwind_protect (delete_input_buffer, (void *) new_buf); + + switch_to_buffer (new_buf); + + unwind_protect_ptr (curr_sym_tab); + + reset_parser (); - return retval; -} + parse_status = yyparse (); -tree_constant -tree_constant_rep::do_matrix_index (int i, const tree_constant& j_arg) const -{ +// Important to reset the idea of where input is coming from before +// trying to eval the command we just parsed -- it might contain the +// name of an function file that still needs to be parsed! + + tree *command = global_command; + + run_unwind_frame ("eval_string"); + tree_constant retval; - tree_constant tmp_j = j_arg.make_numeric_or_range_or_magic (); - - tree_constant_rep::constant_type jtype = tmp_j.const_type (); - - int nr = rows (); - int nc = columns (); - - switch (jtype) + if (parse_status == 0 && command != NULL_TREE) { - case complex_scalar_constant: - case scalar_constant: - { - int j = tree_to_mat_idx (tmp_j.double_value ()); - if (index_check (j, "column") < 0) - return tree_constant (); - if (range_max_check (i, j, nr, nc) < 0) - return tree_constant (); - retval = do_matrix_index (i, j); - } - break; - case complex_matrix_constant: - case matrix_constant: - { - Matrix mj = tmp_j.matrix_value (); - idx_vector jv (mj, user_pref.do_fortran_indexing, "column", nc); - if (! jv) - return tree_constant (); - - if (jv.length () == 0) - { - Matrix mtmp; - retval = tree_constant (mtmp); - } - else - { - if (range_max_check (i, jv.max (), nr, nc) < 0) - return tree_constant (); - retval = do_matrix_index (i, jv); - } - } - break; - case string_constant: - gripe_string_invalid (); - break; - case range_constant: - { - Range rj = tmp_j.range_value (); - if (nc == 2 && is_zero_one (rj)) - { - retval = do_matrix_index (i, 1); - } - else if (nc == 2 && is_one_zero (rj)) - { - retval = do_matrix_index (i, 0); - } - else - { - if (index_check (rj, "column") < 0) - return tree_constant (); - if (range_max_check (i, tree_to_mat_idx (rj.max ()), nr, nc) < 0) - return tree_constant (); - retval = do_matrix_index (i, rj); - } - } - break; - case magic_colon: - if (range_max_check (i, 0, nr, nc) < 0) - return tree_constant (); - retval = do_matrix_index (i, magic_colon); - break; - default: - panic_impossible (); - break; - } - - return retval; -} - -tree_constant -tree_constant_rep::do_matrix_index (const idx_vector& iv, - const tree_constant& j_arg) const -{ - tree_constant retval; - - tree_constant tmp_j = j_arg.make_numeric_or_range_or_magic (); - - tree_constant_rep::constant_type jtype = tmp_j.const_type (); - - int nr = rows (); - int nc = columns (); - - switch (jtype) - { - case complex_scalar_constant: - case scalar_constant: - { - int j = tree_to_mat_idx (tmp_j.double_value ()); - if (index_check (j, "column") < 0) - return tree_constant (); - if (range_max_check (iv.max (), j, nr, nc) < 0) - return tree_constant (); - retval = do_matrix_index (iv, j); - } - break; - case complex_matrix_constant: - case matrix_constant: - { - Matrix mj = tmp_j.matrix_value (); - idx_vector jv (mj, user_pref.do_fortran_indexing, "column", nc); - if (! jv) - return tree_constant (); - - if (jv.length () == 0) - { - Matrix mtmp; - retval = tree_constant (mtmp); - } - else - { - if (range_max_check (iv.max (), jv.max (), nr, nc) < 0) - return tree_constant (); - retval = do_matrix_index (iv, jv); - } - } - break; - case string_constant: - gripe_string_invalid (); - break; - case range_constant: - { - Range rj = tmp_j.range_value (); - if (nc == 2 && is_zero_one (rj)) - { - retval = do_matrix_index (iv, 1); - } - else if (nc == 2 && is_one_zero (rj)) - { - retval = do_matrix_index (iv, 0); - } - else - { - if (index_check (rj, "column") < 0) - return tree_constant (); - if (range_max_check (iv.max (), tree_to_mat_idx (rj.max ()), - nr, nc) < 0) - return tree_constant (); - retval = do_matrix_index (iv, rj); - } - } - break; - case magic_colon: - if (range_max_check (iv.max (), 0, nr, nc) < 0) - return tree_constant (); - retval = do_matrix_index (iv, magic_colon); - break; - default: - panic_impossible (); - break; + retval = command->eval (print); + delete command; } return retval; } tree_constant -tree_constant_rep::do_matrix_index (const Range& ri, - const tree_constant& j_arg) const +eval_string (const tree_constant& arg, int& parse_status) { - tree_constant retval; - - tree_constant tmp_j = j_arg.make_numeric_or_range_or_magic (); - - tree_constant_rep::constant_type jtype = tmp_j.const_type (); - - int nr = rows (); - int nc = columns (); - - switch (jtype) + if (! arg.is_string_type ()) { - case complex_scalar_constant: - case scalar_constant: - { - int j = tree_to_mat_idx (tmp_j.double_value ()); - if (index_check (j, "column") < 0) - return tree_constant (); - if (range_max_check (tree_to_mat_idx (ri.max ()), j, nr, nc) < 0) - return tree_constant (); - retval = do_matrix_index (ri, j); - } - break; - case complex_matrix_constant: - case matrix_constant: - { - Matrix mj = tmp_j.matrix_value (); - idx_vector jv (mj, user_pref.do_fortran_indexing, "column", nc); - if (! jv) - return tree_constant (); - - if (jv.length () == 0) - { - Matrix mtmp; - retval = tree_constant (mtmp); - } - else - { - if (range_max_check (tree_to_mat_idx (ri.max ()), - jv.max (), nr, nc) < 0) - return tree_constant (); - retval = do_matrix_index (ri, jv); - } - } - break; - case string_constant: - gripe_string_invalid (); - break; - case range_constant: - { - Range rj = tmp_j.range_value (); - if (nc == 2 && is_zero_one (rj)) - { - retval = do_matrix_index (ri, 1); - } - else if (nc == 2 && is_one_zero (rj)) - { - retval = do_matrix_index (ri, 0); - } - else - { - if (index_check (rj, "column") < 0) - return tree_constant (); - if (range_max_check (tree_to_mat_idx (ri.max ()), - tree_to_mat_idx (rj.max ()), nr, nc) < 0) - return tree_constant (); - retval = do_matrix_index (ri, rj); - } - } - break; - case magic_colon: - retval = do_matrix_index (ri, magic_colon); - break; - default: - panic_impossible (); - break; + error ("eval: expecting string argument"); + return -1; } - return retval; + char *string = arg.string_value (); + +// Yes Virginia, we always print here... + + return eval_string (string, 1, 1, parse_status); +} + +static int +match_sans_spaces (const char *standard, const char *test) +{ + const char *tp = test; + while (*tp == ' ' || *tp == '\t') + tp++; + + const char *ep = test + strlen (test) - 1; + while (*ep == ' ' || *ep == '\t') + ep--; + + int len = ep - tp + 1; + + return (strncmp (standard, tp, len) == 0); } tree_constant -tree_constant_rep::do_matrix_index (tree_constant_rep::constant_type mci, - const tree_constant& j_arg) const +get_user_input (const tree_constant *args, int nargin, int nargout, + int debug = 0) { tree_constant retval; - tree_constant tmp_j = j_arg.make_numeric_or_range_or_magic (); - - tree_constant_rep::constant_type jtype = tmp_j.const_type (); - - int nr = rows (); - int nc = columns (); - - switch (jtype) + int read_as_string = 0; + if (nargin == 3) { - case complex_scalar_constant: - case scalar_constant: - { - int j = tree_to_mat_idx (tmp_j.double_value ()); - if (index_check (j, "column") < 0) - return tree_constant (); - if (range_max_check (0, j, nr, nc) < 0) - return tree_constant (); - retval = do_matrix_index (magic_colon, j); - } - break; - case complex_matrix_constant: - case matrix_constant: - { - Matrix mj = tmp_j.matrix_value (); - idx_vector jv (mj, user_pref.do_fortran_indexing, "column", nc); - if (! jv) - return tree_constant (); - - if (jv.length () == 0) - { - Matrix mtmp; - retval = tree_constant (mtmp); - } - else - { - if (range_max_check (0, jv.max (), nr, nc) < 0) - return tree_constant (); - retval = do_matrix_index (magic_colon, jv); - } - } - break; - case string_constant: - gripe_string_invalid (); - break; - case range_constant: - { - Range rj = tmp_j.range_value (); - if (nc == 2 && is_zero_one (rj)) - { - retval = do_matrix_index (magic_colon, 1); - } - else if (nc == 2 && is_one_zero (rj)) - { - retval = do_matrix_index (magic_colon, 0); - } - else - { - if (index_check (rj, "column") < 0) - return tree_constant (); - if (range_max_check (0, tree_to_mat_idx (rj.max ()), nr, nc) < 0) - return tree_constant (); - retval = do_matrix_index (magic_colon, rj); - } - } - break; - case magic_colon: - retval = do_matrix_index (magic_colon, magic_colon); - break; - default: - panic_impossible (); - break; - } - - return retval; -} - -tree_constant -tree_constant_rep::do_matrix_index (int i, int j) const -{ - tree_constant retval; - - if (type_tag == matrix_constant) - retval = tree_constant (matrix->elem (i, j)); - else - retval = tree_constant (complex_matrix->elem (i, j)); - - return retval; -} - -tree_constant -tree_constant_rep::do_matrix_index (int i, const idx_vector& jv) const -{ - tree_constant retval; - - int jlen = jv.capacity (); - - CRMATRIX (m, cm, 1, jlen); - - for (int j = 0; j < jlen; j++) - { - int col = jv.elem (j); - CRMATRIX_ASSIGN_REP_ELEM (m, cm, 0, j, i, col); - } - ASSIGN_CRMATRIX_TO (retval, m, cm); - - return retval; -} - -tree_constant -tree_constant_rep::do_matrix_index (int i, const Range& rj) const -{ - tree_constant retval; - - int jlen = rj.nelem (); - - CRMATRIX (m, cm, 1, jlen); - - double b = rj.base (); - double increment = rj.inc (); - for (int j = 0; j < jlen; j++) - { - double tmp = b + j * increment; - int col = tree_to_mat_idx (tmp); - CRMATRIX_ASSIGN_REP_ELEM (m, cm, 0, j, i, col); - } - - ASSIGN_CRMATRIX_TO (retval, m, cm); - - return retval; -} - -tree_constant -tree_constant_rep::do_matrix_index - (int i, tree_constant_rep::constant_type mcj) const -{ - assert (mcj == magic_colon); - - tree_constant retval; - - int nc = columns (); - - CRMATRIX (m, cm, 1, nc); - - for (int j = 0; j < nc; j++) - { - CRMATRIX_ASSIGN_REP_ELEM (m, cm, 0, j, i, j); - } - - ASSIGN_CRMATRIX_TO (retval, m, cm); - - return retval; -} - -tree_constant -tree_constant_rep::do_matrix_index (const idx_vector& iv, int j) const -{ - tree_constant retval; - - int ilen = iv.capacity (); - - CRMATRIX (m, cm, ilen, 1); - - for (int i = 0; i < ilen; i++) - { - int row = iv.elem (i); - CRMATRIX_ASSIGN_REP_ELEM (m, cm, i, 0, row, j); - } - - ASSIGN_CRMATRIX_TO (retval, m, cm); - - return retval; -} - -tree_constant -tree_constant_rep::do_matrix_index (const idx_vector& iv, - const idx_vector& jv) const -{ - tree_constant retval; - - int ilen = iv.capacity (); - int jlen = jv.capacity (); - - CRMATRIX (m, cm, ilen, jlen); - - for (int i = 0; i < ilen; i++) - { - int row = iv.elem (i); - for (int j = 0; j < jlen; j++) + if (args[2].is_string_type () + && strcmp ("s", args[2].string_value ()) == 0) + read_as_string++; + else { - int col = jv.elem (j); - CRMATRIX_ASSIGN_REP_ELEM (m, cm, i, j, row, col); + error ("input: unrecognized second argument"); + return retval; } } - ASSIGN_CRMATRIX_TO (retval, m, cm); - - return retval; -} - -tree_constant -tree_constant_rep::do_matrix_index (const idx_vector& iv, - const Range& rj) const -{ - tree_constant retval; - - int ilen = iv.capacity (); - int jlen = rj.nelem (); - - CRMATRIX (m, cm, ilen, jlen); - - double b = rj.base (); - double increment = rj.inc (); - - for (int i = 0; i < ilen; i++) - { - int row = iv.elem (i); - for (int j = 0; j < jlen; j++) - { - double tmp = b + j * increment; - int col = tree_to_mat_idx (tmp); - CRMATRIX_ASSIGN_REP_ELEM (m, cm, i, j, row, col); - } - } - - ASSIGN_CRMATRIX_TO (retval, m, cm); - - return retval; -} - -tree_constant -tree_constant_rep::do_matrix_index - (const idx_vector& iv, tree_constant_rep::constant_type mcj) const -{ - assert (mcj == magic_colon); - - tree_constant retval; - - int nc = columns (); - int ilen = iv.capacity (); - - CRMATRIX (m, cm, ilen, nc); - - for (int j = 0; j < nc; j++) - { - for (int i = 0; i < ilen; i++) + char *prompt = "debug> "; + if (nargin > 1) + { + if (args[1].is_string_type ()) + prompt = args[1].string_value (); + else { - int row = iv.elem (i); - CRMATRIX_ASSIGN_REP_ELEM (m, cm, i, j, row, j); - } - } - - ASSIGN_CRMATRIX_TO (retval, m, cm); - - return retval; -} - -tree_constant -tree_constant_rep::do_matrix_index (const Range& ri, int j) const -{ - tree_constant retval; - - int ilen = ri.nelem (); - - CRMATRIX (m, cm, ilen, 1); - - double b = ri.base (); - double increment = ri.inc (); - for (int i = 0; i < ilen; i++) - { - double tmp = b + i * increment; - int row = tree_to_mat_idx (tmp); - CRMATRIX_ASSIGN_REP_ELEM (m, cm, i, 0, row, j); - } - - ASSIGN_CRMATRIX_TO (retval, m, cm); - - return retval; -} - -tree_constant -tree_constant_rep::do_matrix_index (const Range& ri, - const idx_vector& jv) const -{ - tree_constant retval; - - int ilen = ri.nelem (); - int jlen = jv.capacity (); - - CRMATRIX (m, cm, ilen, jlen); - - double b = ri.base (); - double increment = ri.inc (); - for (int i = 0; i < ilen; i++) - { - double tmp = b + i * increment; - int row = tree_to_mat_idx (tmp); - for (int j = 0; j < jlen; j++) - { - int col = jv.elem (j); - CRMATRIX_ASSIGN_REP_ELEM (m, cm, i, j, row, col); - } - } - - ASSIGN_CRMATRIX_TO (retval, m, cm); - - return retval; -} - -tree_constant -tree_constant_rep::do_matrix_index (const Range& ri, const Range& rj) const -{ - tree_constant retval; - - int ilen = ri.nelem (); - int jlen = rj.nelem (); - - CRMATRIX (m, cm, ilen, jlen); - - double ib = ri.base (); - double iinc = ri.inc (); - double jb = rj.base (); - double jinc = rj.inc (); - - for (int i = 0; i < ilen; i++) - { - double itmp = ib + i * iinc; - int row = tree_to_mat_idx (itmp); - for (int j = 0; j < jlen; j++) - { - double jtmp = jb + j * jinc; - int col = tree_to_mat_idx (jtmp); - - CRMATRIX_ASSIGN_REP_ELEM (m, cm, i, j, row, col); + error ("input: unrecognized argument"); + return retval; } } - ASSIGN_CRMATRIX_TO (retval, m, cm); + again: - return retval; -} + flush_output_to_pager (); + + char *input_buf = gnu_readline (prompt); -tree_constant -tree_constant_rep::do_matrix_index - (const Range& ri, tree_constant_rep::constant_type mcj) const -{ - assert (mcj == magic_colon); + if (input_buf != (char *) NULL) + { + if (input_buf) + maybe_save_history (input_buf); - tree_constant retval; - - int nc = columns (); + int len = strlen (input_buf); - int ilen = ri.nelem (); - - CRMATRIX (m, cm, ilen, nc); - - double ib = ri.base (); - double iinc = ri.inc (); + if (len < 1) + { + if (debug) + goto again; + else + return retval; + } - for (int i = 0; i < ilen; i++) - { - double itmp = ib + i * iinc; - int row = tree_to_mat_idx (itmp); - for (int j = 0; j < nc; j++) + if (match_sans_spaces ("exit", input_buf) + || match_sans_spaces ("quit", input_buf) + || match_sans_spaces ("return", input_buf)) + return tree_constant (); + else if (read_as_string) + retval = tree_constant (input_buf); + else { - CRMATRIX_ASSIGN_REP_ELEM (m, cm, i, j, row, j); + int parse_status; + retval = eval_string (input_buf, 0, 0, parse_status); + if (debug && retval.is_defined ()) + retval.eval (1); } } - - ASSIGN_CRMATRIX_TO (retval, m, cm); - - return retval; -} - -tree_constant -tree_constant_rep::do_matrix_index (tree_constant_rep::constant_type mci, - int j) const -{ - assert (mci == magic_colon); - - tree_constant retval; - - int nr = rows (); - - CRMATRIX (m, cm, nr, 1); - - for (int i = 0; i < nr; i++) - { - CRMATRIX_ASSIGN_REP_ELEM (m, cm, i, 0, i, j); - } - - ASSIGN_CRMATRIX_TO (retval, m, cm); - - return retval; -} - -tree_constant -tree_constant_rep::do_matrix_index (tree_constant_rep::constant_type mci, - const idx_vector& jv) const -{ - assert (mci == magic_colon); - - tree_constant retval; - - int nr = rows (); - int jlen = jv.capacity (); - - CRMATRIX (m, cm, nr, jlen); - - for (int i = 0; i < nr; i++) - { - for (int j = 0; j < jlen; j++) - { - int col = jv.elem (j); - CRMATRIX_ASSIGN_REP_ELEM (m, cm, i, j, i, col); - } - } - - ASSIGN_CRMATRIX_TO (retval, m, cm); - - return retval; -} + else + error ("input: reading user-input failed!"); -tree_constant -tree_constant_rep::do_matrix_index (tree_constant_rep::constant_type mci, - const Range& rj) const -{ - assert (mci == magic_colon); - - tree_constant retval; - - int nr = rows (); - int jlen = rj.nelem (); - - CRMATRIX (m, cm, nr, jlen); - - double jb = rj.base (); - double jinc = rj.inc (); - - for (int j = 0; j < jlen; j++) - { - double jtmp = jb + j * jinc; - int col = tree_to_mat_idx (jtmp); - for (int i = 0; i < nr; i++) - { - CRMATRIX_ASSIGN_REP_ELEM (m, cm, i, j, i, col); - } - } - - ASSIGN_CRMATRIX_TO (retval, m, cm); - - return retval; -} + if (debug) + goto again; -tree_constant -tree_constant_rep::do_matrix_index (tree_constant_rep::constant_type mci, - tree_constant_rep::constant_type mcj) const -{ - assert (mci == magic_colon && mcj == magic_colon); - - return tree_constant (*this); -} - -tree_constant -tree_constant_rep::do_matrix_index - (tree_constant_rep::constant_type mci) const -{ - assert (mci == magic_colon); - - tree_constant retval; - int nr = rows (); - int nc = columns (); - int size = nr * nc; - if (size > 0) - { - CRMATRIX (m, cm, size, 1); - int idx = 0; - for (int j = 0; j < nc; j++) - for (int i = 0; i < nr; i++) - { - CRMATRIX_ASSIGN_REP_ELEM (m, cm, idx, 0, i, j); - idx++; - } - ASSIGN_CRMATRIX_TO (retval, m, cm); - } return retval; }