X-Git-Url: https://www.ginac.de/ginac.git//ginac.git?p=ginac.git;a=blobdiff_plain;f=ginac%2Fmatrix.cpp;h=710190376a1c28b59c787d902655aa24d97ccf73;hp=c9b3521a4665cdd8f210782809e244288de7efa1;hb=0160f9ab1da453641e30539abcf0eaa4162582eb;hpb=b66548802c56b34d6b212a0196d622937841ca61 diff --git a/ginac/matrix.cpp b/ginac/matrix.cpp index c9b3521a..71019037 100644 --- a/ginac/matrix.cpp +++ b/ginac/matrix.cpp @@ -3,7 +3,7 @@ * Implementation of symbolic matrices */ /* - * GiNaC Copyright (C) 1999-2003 Johannes Gutenberg University Mainz, Germany + * GiNaC Copyright (C) 1999-2005 Johannes Gutenberg University Mainz, Germany * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by @@ -17,7 +17,7 @@ * * You should have received a copy of the GNU General Public License * along with this program; if not, write to the Free Software - * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA + * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA */ #include @@ -45,7 +45,7 @@ namespace GiNaC { GINAC_IMPLEMENT_REGISTERED_CLASS_OPT(matrix, basic, print_func(&matrix::do_print). print_func(&matrix::do_print_latex). - print_func(&basic::do_print_tree). + print_func(&matrix::do_print_tree). print_func(&matrix::do_print_python_repr)) ////////// @@ -53,9 +53,9 @@ GINAC_IMPLEMENT_REGISTERED_CLASS_OPT(matrix, basic, ////////// /** Default ctor. Initializes to 1 x 1-dimensional zero-matrix. */ -matrix::matrix() : inherited(TINFO_matrix), row(1), col(1) +matrix::matrix() : inherited(TINFO_matrix), row(1), col(1), m(1, _ex0) { - m.push_back(_ex0); + setflag(status_flags::not_shareable); } ////////// @@ -69,25 +69,28 @@ matrix::matrix() : inherited(TINFO_matrix), row(1), col(1) * @param r number of rows * @param c number of cols */ matrix::matrix(unsigned r, unsigned c) - : inherited(TINFO_matrix), row(r), col(c) + : inherited(TINFO_matrix), row(r), col(c), m(r*c, _ex0) { - m.resize(r*c, _ex0); + setflag(status_flags::not_shareable); } // protected /** Ctor from representation, for internal use only. */ matrix::matrix(unsigned r, unsigned c, const exvector & m2) - : inherited(TINFO_matrix), row(r), col(c), m(m2) {} + : inherited(TINFO_matrix), row(r), col(c), m(m2) +{ + setflag(status_flags::not_shareable); +} /** Construct matrix from (flat) list of elements. If the list has fewer * elements than the matrix, the remaining matrix elements are set to zero. * If the list has more elements than the matrix, the excessive elements are * thrown away. */ matrix::matrix(unsigned r, unsigned c, const lst & l) - : inherited(TINFO_matrix), row(r), col(c) + : inherited(TINFO_matrix), row(r), col(c), m(r*c, _ex0) { - m.resize(r*c, _ex0); + setflag(status_flags::not_shareable); size_t i = 0; for (lst::const_iterator it = l.begin(); it != l.end(); ++it, ++i) { @@ -105,6 +108,8 @@ matrix::matrix(unsigned r, unsigned c, const lst & l) matrix::matrix(const archive_node &n, lst &sym_lst) : inherited(n, sym_lst) { + setflag(status_flags::not_shareable); + if (!(n.find_unsigned("row", row)) || !(n.find_unsigned("col", col))) throw (std::runtime_error("unknown matrix dimensions in archive")); m.reserve(row * col); @@ -137,7 +142,7 @@ DEFAULT_UNARCHIVE(matrix) // public -void matrix::print_elements(const print_context & c, const std::string & row_start, const std::string & row_end, const std::string & row_sep, const std::string & col_sep) const +void matrix::print_elements(const print_context & c, const char *row_start, const char *row_end, const char *row_sep, const char *col_sep) const { for (unsigned ro=0; roconjugate(); + if (ev) { + ev->push_back(x); + continue; + } + if (are_ex_trivially_equal(x, *i)) { + continue; + } + ev = new exvector; + ev->reserve(m.size()); + for (exvector::const_iterator j=m.begin(); j!=i; ++j) { + ev->push_back(*j); + } + ev->push_back(x); + } + if (ev) { + ex result = matrix(row, col, *ev); + delete ev; + return result; + } + return *this; +} + // protected int matrix::compare_same_type(const basic & other) const @@ -614,12 +647,12 @@ matrix matrix::pow(const ex & expn) const // that this is not entirely optimal but close to optimal and // "better" algorithms are much harder to implement. (See Knuth, // TAoCP2, section "Evaluation of Powers" for a good discussion.) - while (b!=_num1) { + while (b!=*_num1_p) { if (b.is_odd()) { C = C.mul(A); --b; } - b /= _num2; // still integer. + b /= *_num2_p; // still integer. A = A.mul(A); } return A.mul(C); @@ -697,12 +730,12 @@ ex matrix::determinant(unsigned algo) const unsigned sparse_count = 0; // counts non-zero elements exvector::const_iterator r = m.begin(), rend = m.end(); while (r != rend) { - lst srl; // symbol replacement list + if (!r->info(info_flags::numeric)) + numeric_flag = false; + exmap srl; // symbol replacement list ex rtest = r->to_rational(srl); if (!rtest.is_zero()) ++sparse_count; - if (!rtest.info(info_flags::numeric)) - numeric_flag = false; if (!rtest.info(info_flags::crational_polynomial) && rtest.info(info_flags::rational_function)) normal_flag = true; @@ -731,7 +764,7 @@ ex matrix::determinant(unsigned algo) const else return m[0].expand(); } - + // Compute the determinant switch(algo) { case determinant_algo::gauss: { @@ -827,7 +860,7 @@ ex matrix::trace() const tr += m[r*col+r]; if (tr.info(info_flags::rational_function) && - !tr.info(info_flags::crational_polynomial)) + !tr.info(info_flags::crational_polynomial)) return tr.normal(); else return tr.expand(); @@ -845,7 +878,7 @@ ex matrix::trace() const * @return characteristic polynomial as new expression * @exception logic_error (matrix not square) * @see matrix::determinant() */ -ex matrix::charpoly(const symbol & lambda) const +ex matrix::charpoly(const ex & lambda) const { if (row != col) throw (std::logic_error("matrix::charpoly(): matrix not square")); @@ -865,13 +898,13 @@ ex matrix::charpoly(const symbol & lambda) const matrix B(*this); ex c = B.trace(); - ex poly = power(lambda,row)-c*power(lambda,row-1); + ex poly = power(lambda, row) - c*power(lambda, row-1); for (unsigned i=1; imul(B); c = B.trace() / ex(i+1); - poly -= c*power(lambda,row-i-1); + poly -= c*power(lambda, row-i-1); } if (row%2) return -poly; @@ -933,14 +966,15 @@ matrix matrix::inverse() const * * @param vars n x p matrix, all elements must be symbols * @param rhs m x p matrix + * @param algo selects the solving algorithm * @return n x p solution matrix * @exception logic_error (incompatible matrices) * @exception invalid_argument (1st argument must be matrix of symbols) * @exception runtime_error (inconsistent linear system) * @see solve_algo */ matrix matrix::solve(const matrix & vars, - const matrix & rhs, - unsigned algo) const + const matrix & rhs, + unsigned algo) const { const unsigned m = this->rows(); const unsigned n = this->cols(); @@ -1033,6 +1067,29 @@ matrix matrix::solve(const matrix & vars, } +/** Compute the rank of this matrix. */ +unsigned matrix::rank() const +{ + // Method: + // Transform this matrix into upper echelon form and then count the + // number of non-zero rows. + + GINAC_ASSERT(row*col==m.capacity()); + + // Actually, any elimination scheme will do since we are only + // interested in the echelon matrix' zeros. + matrix to_eliminate = *this; + to_eliminate.fraction_free_elimination(); + + unsigned r = row*col; // index of last non-zero element + while (r--) { + if (!to_eliminate.m[r].is_zero()) + return 1+r/col; + } + return 0; +} + + // protected /** Recursive determinant for small matrices having at least one symbolic @@ -1148,7 +1205,7 @@ ex matrix::determinant_minor() const Pkey[j] = Pkey[j-1]+1; } while(fc); // next column, so change the role of A and B: - A = B; + A.swap(B); B.clear(); } @@ -1174,8 +1231,8 @@ int matrix::gauss_elimination(const bool det) int sign = 1; unsigned r0 = 0; - for (unsigned r1=0; (r1 0) sign = -sign; for (unsigned r2=r0+1; r2m[r2*n+r1].is_zero()) { + if (!this->m[r2*n+c0].is_zero()) { // yes, there is something to do in this row - ex piv = this->m[r2*n+r1] / this->m[r0*n+r1]; - for (unsigned c=r1+1; cm[r2*n+c0] / this->m[r0*n+c0]; + for (unsigned c=c0+1; cm[r2*n+c] -= piv * this->m[r0*n+c]; if (!this->m[r2*n+c].info(info_flags::numeric)) this->m[r2*n+c] = this->m[r2*n+c].normal(); } } // fill up left hand side with zeros - for (unsigned c=0; c<=r1; ++c) + for (unsigned c=r0; c<=c0; ++c) this->m[r2*n+c] = _ex0; } if (det) { @@ -1206,7 +1263,12 @@ int matrix::gauss_elimination(const bool det) ++r0; } } - + // clear remaining rows + for (unsigned r=r0+1; rm[r*n+c] = _ex0; + } + return sign; } @@ -1228,8 +1290,8 @@ int matrix::division_free_elimination(const bool det) int sign = 1; unsigned r0 = 0; - for (unsigned r1=0; (r10) sign = -sign; for (unsigned r2=r0+1; r2m[r2*n+c] = (this->m[r0*n+r1]*this->m[r2*n+c] - this->m[r2*n+r1]*this->m[r0*n+c]).expand(); + for (unsigned c=c0+1; cm[r2*n+c] = (this->m[r0*n+c0]*this->m[r2*n+c] - this->m[r2*n+c0]*this->m[r0*n+c]).expand(); // fill up left hand side with zeros - for (unsigned c=0; c<=r1; ++c) + for (unsigned c=r0; c<=c0; ++c) this->m[r2*n+c] = _ex0; } if (det) { @@ -1253,7 +1315,12 @@ int matrix::division_free_elimination(const bool det) ++r0; } } - + // clear remaining rows + for (unsigned r=r0+1; rm[r*n+c] = _ex0; + } + return sign; } @@ -1315,7 +1382,7 @@ int matrix::fraction_free_elimination(const bool det) // makes things more complicated than they need to be. matrix tmp_n(*this); matrix tmp_d(m,n); // for denominators, if needed - lst srl; // symbol replacement list + exmap srl; // symbol replacement list exvector::const_iterator cit = this->m.begin(), citend = this->m.end(); exvector::iterator tmp_n_it = tmp_n.m.begin(), tmp_d_it = tmp_d.m.begin(); while (cit != citend) { @@ -1326,8 +1393,8 @@ int matrix::fraction_free_elimination(const bool det) } unsigned r0 = 0; - for (unsigned r1=0; (r10) { sign = -sign; // tmp_n's rows r0 and indx were swapped, do the same in tmp_d: - for (unsigned c=r1; cm.begin(), itend = this->m.end(); tmp_n_it = tmp_n.m.begin(); @@ -1525,4 +1598,52 @@ ex symbolic_matrix(unsigned r, unsigned c, const std::string & base_name, const return M; } +ex reduced_matrix(const matrix& m, unsigned r, unsigned c) +{ + if (r+1>m.rows() || c+1>m.cols() || m.cols()<2 || m.rows()<2) + throw std::runtime_error("minor_matrix(): index out of bounds"); + + const unsigned rows = m.rows()-1; + const unsigned cols = m.cols()-1; + matrix &M = *new matrix(rows, cols); + M.setflag(status_flags::dynallocated | status_flags::evaluated); + + unsigned ro = 0; + unsigned ro2 = 0; + while (ro2m.rows() || c+nc>m.cols()) + throw std::runtime_error("sub_matrix(): index out of bounds"); + + matrix &M = *new matrix(nr, nc); + M.setflag(status_flags::dynallocated | status_flags::evaluated); + + for (unsigned ro=0; ro