From 08d556dc3ac3fbf2b0ad3acd37016a1f925d7c02 Mon Sep 17 00:00:00 2001 From: Christian Bauer Date: Fri, 23 Mar 2001 17:49:54 +0000 Subject: [PATCH] - simplify_indexed() simplifies c*(M).i.j -> (M').i.j (c: numeric, M: matrix) with M'=c*M - new member function matrix::mul(const numeric &other) to multiply matrix with scalar - lst_to_matrix() takes a "lst &" instead of an "ex &" - added constructor of matrix from flat list - added function diag_matrix() to construct matrix from list of diagonal elements --- ginac/basic.cpp | 12 ++++++++ ginac/basic.h | 1 + ginac/indexed.cpp | 14 ++++++++-- ginac/input_parser.yy | 2 +- ginac/matrix.cpp | 65 +++++++++++++++++++++++++++++++++++++++---- ginac/matrix.h | 9 +++++- 6 files changed, 93 insertions(+), 10 deletions(-) diff --git a/ginac/basic.cpp b/ginac/basic.cpp index 6415f33f..1c6fe6fa 100644 --- a/ginac/basic.cpp +++ b/ginac/basic.cpp @@ -310,6 +310,18 @@ ex basic::add_indexed(const ex & self, const ex & other) const return self + other; } +/** Multiply an indexed expression with a scalar. This function is used + * internally by simplify_indexed(). + * + * @param self Indexed expression; it's base object is *this + * @param other Numeric value + * @return product of self and other + * @see ex::simplify_indexed() */ +ex basic::scalar_mul_indexed(const ex & self, const numeric & other) const +{ + return self * other; +} + /** Try to contract two indexed expressions that appear in the same product. * If a contraction exists, the function overwrites one or both of the * expressions and returns true. Otherwise it returns false. It is diff --git a/ginac/basic.h b/ginac/basic.h index 4277d949..2e1aefdf 100644 --- a/ginac/basic.h +++ b/ginac/basic.h @@ -131,6 +131,7 @@ public: // only const functions please (may break reference counting) virtual ex simplify_ncmul(const exvector & v) const; virtual ex eval_indexed(const basic & i) const; virtual ex add_indexed(const ex & self, const ex & other) const; + virtual ex scalar_mul_indexed(const ex & self, const numeric & other) const; virtual bool contract_with(exvector::iterator self, exvector::iterator other, exvector & v) const; protected: // non-const functions should be called from class ex only virtual ex derivative(const symbol & s) const; diff --git a/ginac/indexed.cpp b/ginac/indexed.cpp index e23709d5..86fe6115 100644 --- a/ginac/indexed.cpp +++ b/ginac/indexed.cpp @@ -679,13 +679,21 @@ try_again: } find_free_and_dummy(un.begin(), un.end(), free_indices, dummy_indices); + ex r; if (something_changed) { if (non_commutative) - return ncmul(v); + r = ncmul(v); else - return mul(v); + r = mul(v); } else - return e; + r = e; + + // Product of indexed object with a scalar? + if (is_ex_exactly_of_type(r, mul) && r.nops() == 2 + && is_ex_exactly_of_type(r.op(1), numeric) && is_ex_of_type(r.op(0), indexed)) + return r.op(0).op(0).bp->scalar_mul_indexed(r.op(0), ex_to_numeric(r.op(1))); + else + return r; } /** Simplify indexed expression, return list of free indices. */ diff --git a/ginac/input_parser.yy b/ginac/input_parser.yy index d1e87669..9c4ffd97 100644 --- a/ginac/input_parser.yy +++ b/ginac/input_parser.yy @@ -115,7 +115,7 @@ exp : T_NUMBER {$$ = $1;} | exp '!' {$$ = factorial($1);} | '(' exp ')' {$$ = $2;} | '[' list_or_empty ']' {$$ = $2;} - | T_MATRIX_BEGIN matrix T_MATRIX_END {$$ = lst_to_matrix($2);} + | T_MATRIX_BEGIN matrix T_MATRIX_END {$$ = lst_to_matrix(ex_to_lst($2));} ; exprseq : exp {$$ = exprseq($1);} diff --git a/ginac/matrix.cpp b/ginac/matrix.cpp index f34534af..650799e5 100644 --- a/ginac/matrix.cpp +++ b/ginac/matrix.cpp @@ -95,6 +95,25 @@ matrix::matrix(unsigned r, unsigned c, const exvector & m2) debugmsg("matrix ctor from unsigned,unsigned,exvector",LOGLEVEL_CONSTRUCT); } +/** 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) +{ + debugmsg("matrix ctor from unsigned,unsigned,lst",LOGLEVEL_CONSTRUCT); + m.resize(r*c, _ex0()); + + for (unsigned i=0; i= r) + break; // matrix smaller than list: throw away excessive elements + m[y*c+x] = l.op(i); + } +} + ////////// // archiving ////////// @@ -385,6 +404,7 @@ ex matrix::eval_indexed(const basic & i) const ex matrix::add_indexed(const ex & self, const ex & other) const { GINAC_ASSERT(is_ex_of_type(self, indexed)); + GINAC_ASSERT(is_ex_of_type(self.op(0), matrix)); GINAC_ASSERT(is_ex_of_type(other, indexed)); GINAC_ASSERT(self.nops() == 2 || self.nops() == 3); @@ -416,6 +436,21 @@ ex matrix::add_indexed(const ex & self, const ex & other) const return self + other; } +/** Product of an indexed matrix with a number. */ +ex matrix::scalar_mul_indexed(const ex & self, const numeric & other) const +{ + GINAC_ASSERT(is_ex_of_type(self, indexed)); + GINAC_ASSERT(is_ex_of_type(self.op(0), matrix)); + GINAC_ASSERT(self.nops() == 2 || self.nops() == 3); + + const matrix &self_matrix = ex_to_matrix(self.op(0)); + + if (self.nops() == 2) + return indexed(self_matrix.mul(other), self.op(1)); + else // self.nops() == 3 + return indexed(self_matrix.mul(other), self.op(1), self.op(2)); +} + /** Contraction of an indexed matrix with something else. */ bool matrix::contract_with(exvector::iterator self, exvector::iterator other, exvector & v) const { @@ -581,6 +616,19 @@ matrix matrix::mul(const matrix & other) const } +/** Product of matrix and scalar. */ +matrix matrix::mul(const numeric & other) const +{ + exvector prod(row * col); + + for (unsigned r=0; r(*e.bp); } -extern ex lst_to_matrix(const ex &l); +/** Convert list of lists to matrix. */ +extern ex lst_to_matrix(const lst & l); + +/** Convert list of diagonal elements to matrix. */ +extern ex diag_matrix(const lst & l); } // namespace GiNaC -- 2.44.0