} else {
- c.s << "[[ ";
+ c.s << "[";
for (unsigned y=0; y<row-1; ++y) {
- c.s << "[[";
+ c.s << "[";
for (unsigned x=0; x<col-1; ++x) {
m[y*col+x].print(c);
c.s << ",";
}
m[col*(y+1)-1].print(c);
- c.s << "]], ";
+ c.s << "],";
}
- c.s << "[[";
+ c.s << "[";
for (unsigned x=0; x<col-1; ++x) {
m[(row-1)*col+x].print(c);
c.s << ",";
}
m[row*col-1].print(c);
- c.s << "]] ]]";
+ c.s << "]]";
}
}
return matrix(row, col, tmp);
}
-/** Search ocurrences. A matrix 'has' an expression if it is the expression
- * itself or one of the elements 'has' it. */
-bool matrix::has(const ex & other) const
-{
- GINAC_ASSERT(other.bp!=0);
-
- // tautology: it is the expression itself
- if (is_equal(*other.bp)) return true;
-
- // search all the elements
- for (exvector::const_iterator r=m.begin(); r!=m.end(); ++r)
- if ((*r).has(other)) return true;
-
- return false;
-}
-
/** Evaluate matrix entry by entry. */
ex matrix::eval(int level) const
{
return matrix(row, col, m2);
}
-ex matrix::subs(const lst & ls, const lst & lr) const
+ex matrix::subs(const lst & ls, const lst & lr, bool no_pattern) const
{
exvector m2(row * col);
for (unsigned r=0; r<row; ++r)
for (unsigned c=0; c<col; ++c)
- m2[r*col+c] = m[r*col+c].subs(ls, lr);
+ m2[r*col+c] = m[r*col+c].subs(ls, lr, no_pattern);
- return matrix(row, col, m2);
+ return ex(matrix(row, col, m2)).bp->basic::subs(ls, lr, no_pattern);
}
// protected
matrix matrix::add(const matrix & other) const
{
if (col != other.col || row != other.row)
- throw (std::logic_error("matrix::add(): incompatible matrices"));
+ throw std::logic_error("matrix::add(): incompatible matrices");
exvector sum(this->m);
exvector::iterator i;
matrix matrix::sub(const matrix & other) const
{
if (col != other.col || row != other.row)
- throw (std::logic_error("matrix::sub(): incompatible matrices"));
+ throw std::logic_error("matrix::sub(): incompatible matrices");
exvector dif(this->m);
exvector::iterator i;
matrix matrix::mul(const matrix & other) const
{
if (this->cols() != other.rows())
- throw (std::logic_error("matrix::mul(): incompatible matrices"));
+ throw std::logic_error("matrix::mul(): incompatible matrices");
exvector prod(this->rows()*other.cols());
}
+/** Product of matrix and scalar expression. */
+matrix matrix::mul_scalar(const ex & other) const
+{
+ if (other.return_type() != return_types::commutative)
+ throw std::runtime_error("matrix::mul_scalar(): non-commutative scalar");
+
+ exvector prod(row * col);
+
+ for (unsigned r=0; r<row; ++r)
+ for (unsigned c=0; c<col; ++c)
+ prod[r*col+c] = m[r*col+c] * other;
+
+ return matrix(row, col, prod);
+}
+
+
+/** Power of a matrix. Currently handles integer exponents only. */
+matrix matrix::pow(const ex & expn) const
+{
+ if (col!=row)
+ throw (std::logic_error("matrix::pow(): matrix not square"));
+
+ if (is_ex_exactly_of_type(expn, numeric)) {
+ // Integer cases are computed by successive multiplication, using the
+ // obvious shortcut of storing temporaries, like A^4 == (A*A)*(A*A).
+ if (expn.info(info_flags::integer)) {
+ numeric k;
+ matrix prod(row,col);
+ if (expn.info(info_flags::negative)) {
+ k = -ex_to_numeric(expn);
+ prod = this->inverse();
+ } else {
+ k = ex_to_numeric(expn);
+ prod = *this;
+ }
+ matrix result(row,col);
+ for (unsigned r=0; r<row; ++r)
+ result.set(r,r,_ex1());
+ numeric b(1);
+ // this loop computes the representation of k in base 2 and multiplies
+ // the factors whenever needed:
+ while (b.compare(k)<=0) {
+ b *= numeric(2);
+ numeric r(mod(k,b));
+ if (!r.is_zero()) {
+ k -= r;
+ result = result.mul(prod);
+ }
+ prod = prod.mul(prod);
+ }
+ return result;
+ }
+ }
+ throw (std::runtime_error("matrix::pow(): don't know how to handle exponent"));
+}
+
+
/** operator() to access elements.
*
* @param ro row of element
{
if (ro>=row || co>=col)
throw (std::range_error("matrix::set(): index out of range"));
+ if (value.return_type() != return_types::commutative)
+ throw std::runtime_error("matrix::set(): non-commutative argument");
ensure_if_modifiable();
m[ro*col+co] = value;
return matrix(this->cols(),this->rows(),trans);
}
-
/** Determinant of square matrix. This routine doesn't actually calculate the
* determinant, it only implements some heuristics about which algorithm to
* run. If all the elements of the matrix are elements of an integral domain
std::vector<unsigned> pre_sort;
for (std::vector<uintpair>::iterator i=c_zeros.begin(); i!=c_zeros.end(); ++i)
pre_sort.push_back(i->second);
- int sign = permutation_sign(pre_sort);
+ std::vector<unsigned> pre_sort_test(pre_sort); // permutation_sign() modifies the vector so we make a copy here
+ int sign = permutation_sign(pre_sort_test.begin(), pre_sort_test.end());
exvector result(row*col); // represents sorted matrix
unsigned c = 0;
for (std::vector<unsigned>::iterator i=pre_sort.begin();
if (row != col)
throw (std::logic_error("matrix::inverse(): matrix not square"));
- // NOTE: the Gauss-Jordan elimination used here can in principle be
- // replaced by two clever calls to gauss_elimination() and some to
- // transpose(). Wouldn't be more efficient (maybe less?), just more
- // orthogonal.
- matrix tmp(row,col);
- // set tmp to the unit matrix
- for (unsigned i=0; i<col; ++i)
- tmp.m[i*col+i] = _ex1();
+ // This routine actually doesn't do anything fancy at all. We compute the
+ // inverse of the matrix A by solving the system A * A^{-1} == Id.
- // create a copy of this matrix
- matrix cpy(*this);
- for (unsigned r1=0; r1<row; ++r1) {
- int indx = cpy.pivot(r1, r1);
- if (indx == -1) {
+ // First populate the identity matrix supposed to become the right hand side.
+ matrix identity(row,col);
+ for (unsigned i=0; i<row; ++i)
+ identity.set(i,i,_ex1());
+
+ // Populate a dummy matrix of variables, just because of compatibility with
+ // matrix::solve() which wants this (for compatibility with under-determined
+ // systems of equations).
+ matrix vars(row,col);
+ for (unsigned r=0; r<row; ++r)
+ for (unsigned c=0; c<col; ++c)
+ vars.set(r,c,symbol());
+
+ matrix sol(row,col);
+ try {
+ sol = this->solve(vars,identity);
+ } catch (const std::runtime_error & e) {
+ if (e.what()==std::string("matrix::solve(): inconsistent linear system"))
throw (std::runtime_error("matrix::inverse(): singular matrix"));
- }
- if (indx != 0) { // swap rows r and indx of matrix tmp
- for (unsigned i=0; i<col; ++i)
- tmp.m[r1*col+i].swap(tmp.m[indx*col+i]);
- }
- ex a1 = cpy.m[r1*col+r1];
- for (unsigned c=0; c<col; ++c) {
- cpy.m[r1*col+c] /= a1;
- tmp.m[r1*col+c] /= a1;
- }
- for (unsigned r2=0; r2<row; ++r2) {
- if (r2 != r1) {
- if (!cpy.m[r2*col+r1].is_zero()) {
- ex a2 = cpy.m[r2*col+r1];
- // yes, there is something to do in this column
- for (unsigned c=0; c<col; ++c) {
- cpy.m[r2*col+c] -= a2 * cpy.m[r1*col+c];
- if (!cpy.m[r2*col+c].info(info_flags::numeric))
- cpy.m[r2*col+c] = cpy.m[r2*col+c].normal();
- tmp.m[r2*col+c] -= a2 * tmp.m[r1*col+c];
- if (!tmp.m[r2*col+c].info(info_flags::numeric))
- tmp.m[r2*col+c] = tmp.m[r2*col+c].normal();
- }
- }
- }
- }
+ else
+ throw;
}
-
- return tmp;
+ return sol;
}
switch(algo) {
case solve_algo::gauss:
aug.gauss_elimination();
+ break;
case solve_algo::divfree:
aug.division_free_elimination();
+ break;
case solve_algo::bareiss:
default:
aug.fraction_free_elimination();
git://www.ginac.de / ginac.git / blobdiff