From 539de73dc25de63f2888f8696d73ac3fc83db45f Mon Sep 17 00:00:00 2001 From: Christian Bauer Date: Fri, 9 Mar 2001 22:21:35 +0000 Subject: [PATCH] - lcm_of_coefficients_denominators(1/2+x^10) returned 1024 instead of 2 and multiply_lcm() on that expression barfed - contract_with() now takes two exvector::iterators and and exvector; this makes it possible to look for more contractions in the implementation --- ginac/basic.cpp | 17 +++++----- ginac/basic.h | 2 +- ginac/indexed.cpp | 4 +-- ginac/matrix.cpp | 79 ++++++++++++++++++++++------------------------- ginac/matrix.h | 2 +- ginac/normal.cpp | 13 ++++++-- ginac/tensor.cpp | 54 ++++++++++++++++---------------- ginac/tensor.h | 4 +-- 8 files changed, 89 insertions(+), 86 deletions(-) diff --git a/ginac/basic.cpp b/ginac/basic.cpp index 2414c929..a7645ac5 100644 --- a/ginac/basic.cpp +++ b/ginac/basic.cpp @@ -297,16 +297,17 @@ ex basic::eval_indexed(const basic & i) const return i.hold(); } -/** Try to contract two indexed expressions. If a contraction exists, the - * function overwrites one or both arguments and returns true. Otherwise it - * returns false. It is guaranteed that both expressions are of class - * indexed (or a subclass) and that at least one dummy index has been - * found. +/** 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 + * guaranteed that both expressions are of class indexed (or a subclass) + * and that at least one dummy index has been found. * - * @param self The first indexed expression; it's base object is *this - * @param other The second indexed expression + * @param self Pointer to first indexed expression; it's base object is *this + * @param other Pointer to second indexed expression + * @param v The complete vector of factors * @return true if the contraction was successful, false otherwise */ -bool basic::contract_with(ex & self, ex & other) const +bool basic::contract_with(exvector::iterator self, exvector::iterator other, exvector & v) const { // Do nothing return false; diff --git a/ginac/basic.h b/ginac/basic.h index 842417ca..ee2bef23 100644 --- a/ginac/basic.h +++ b/ginac/basic.h @@ -130,7 +130,7 @@ public: // only const functions please (may break reference counting) virtual exvector get_free_indices(void) const; virtual ex simplify_ncmul(const exvector & v) const; virtual ex eval_indexed(const basic & i) const; - virtual bool contract_with(ex & self, ex & 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; virtual int compare_same_type(const basic & other) const; diff --git a/ginac/indexed.cpp b/ginac/indexed.cpp index ea4266cb..c3fe1d9a 100644 --- a/ginac/indexed.cpp +++ b/ginac/indexed.cpp @@ -629,12 +629,12 @@ try_again: } // Try to contract the first one with the second one - bool contracted = it1->op(0).bp->contract_with(*it1, *it2); + bool contracted = it1->op(0).bp->contract_with(it1, it2, v); if (!contracted) { // That didn't work; maybe the second object knows how to // contract itself with the first one - contracted = it2->op(0).bp->contract_with(*it2, *it1); + contracted = it2->op(0).bp->contract_with(it2, it1, v); } if (contracted) { something_changed = true; diff --git a/ginac/matrix.cpp b/ginac/matrix.cpp index a63275f5..1a5d6953 100644 --- a/ginac/matrix.cpp +++ b/ginac/matrix.cpp @@ -382,103 +382,98 @@ ex matrix::eval_indexed(const basic & i) const } /** Contraction of an indexed matrix with something else. */ -bool matrix::contract_with(ex & self, ex & other) const +bool matrix::contract_with(exvector::iterator self, exvector::iterator other, exvector & v) const { - GINAC_ASSERT(is_ex_of_type(self, indexed)); - GINAC_ASSERT(is_ex_of_type(other, indexed)); - GINAC_ASSERT(self.nops() == 2 || self.nops() == 3); - GINAC_ASSERT(is_ex_of_type(self.op(0), matrix)); + GINAC_ASSERT(is_ex_of_type(*self, indexed)); + GINAC_ASSERT(is_ex_of_type(*other, indexed)); + GINAC_ASSERT(self->nops() == 2 || self->nops() == 3); + GINAC_ASSERT(is_ex_of_type(self->op(0), matrix)); // Only contract with other matrices - if (!is_ex_of_type(other.op(0), matrix)) + if (!is_ex_of_type(other->op(0), matrix)) return false; - GINAC_ASSERT(other.nops() == 2 || other.nops() == 3); + GINAC_ASSERT(other->nops() == 2 || other->nops() == 3); - const matrix &self_matrix = ex_to_matrix(self.op(0)); - const matrix &other_matrix = ex_to_matrix(other.op(0)); + const matrix &self_matrix = ex_to_matrix(self->op(0)); + const matrix &other_matrix = ex_to_matrix(other->op(0)); - if (self.nops() == 2) { + if (self->nops() == 2) { unsigned self_dim = (self_matrix.col == 1) ? self_matrix.row : self_matrix.col; - if (other.nops() == 2) { // vector * vector (scalar product) + if (other->nops() == 2) { // vector * vector (scalar product) unsigned other_dim = (other_matrix.col == 1) ? other_matrix.row : other_matrix.col; if (self_matrix.col == 1) { if (other_matrix.col == 1) { // Column vector * column vector, transpose first vector - self = self_matrix.transpose().mul(other_matrix)(0, 0); + *self = self_matrix.transpose().mul(other_matrix)(0, 0); } else { // Column vector * row vector, swap factors - self = other_matrix.mul(self_matrix)(0, 0); + *self = other_matrix.mul(self_matrix)(0, 0); } } else { if (other_matrix.col == 1) { // Row vector * column vector, perfect - self = self_matrix.mul(other_matrix)(0, 0); + *self = self_matrix.mul(other_matrix)(0, 0); } else { // Row vector * row vector, transpose second vector - self = self_matrix.mul(other_matrix.transpose())(0, 0); + *self = self_matrix.mul(other_matrix.transpose())(0, 0); } } - other = _ex1(); + *other = _ex1(); return true; } else { // vector * matrix - GINAC_ASSERT(other.nops() == 3); - // B_i * A_ij = (B*A)_j (B is row vector) - if (is_dummy_pair(self.op(1), other.op(1))) { + if (is_dummy_pair(self->op(1), other->op(1))) { if (self_matrix.row == 1) - self = indexed(self_matrix.mul(other_matrix), other.op(2)); + *self = indexed(self_matrix.mul(other_matrix), other->op(2)); else - self = indexed(self_matrix.transpose().mul(other_matrix), other.op(2)); - other = _ex1(); + *self = indexed(self_matrix.transpose().mul(other_matrix), other->op(2)); + *other = _ex1(); return true; } // B_j * A_ij = (A*B)_i (B is column vector) - if (is_dummy_pair(self.op(1), other.op(2))) { + if (is_dummy_pair(self->op(1), other->op(2))) { if (self_matrix.col == 1) - self = indexed(other_matrix.mul(self_matrix), other.op(1)); + *self = indexed(other_matrix.mul(self_matrix), other->op(1)); else - self = indexed(other_matrix.mul(self_matrix.transpose()), other.op(1)); - other = _ex1(); + *self = indexed(other_matrix.mul(self_matrix.transpose()), other->op(1)); + *other = _ex1(); return true; } } - } else if (other.nops() == 3) { // matrix * matrix - - GINAC_ASSERT(self.nops() == 3); - GINAC_ASSERT(other.nops() == 3); + } else if (other->nops() == 3) { // matrix * matrix // A_ij * B_jk = (A*B)_ik - if (is_dummy_pair(self.op(2), other.op(1))) { - self = indexed(self_matrix.mul(other_matrix), self.op(1), other.op(2)); - other = _ex1(); + if (is_dummy_pair(self->op(2), other->op(1))) { + *self = indexed(self_matrix.mul(other_matrix), self->op(1), other->op(2)); + *other = _ex1(); return true; } // A_ij * B_kj = (A*Btrans)_ik - if (is_dummy_pair(self.op(2), other.op(2))) { - self = indexed(self_matrix.mul(other_matrix.transpose()), self.op(1), other.op(1)); - other = _ex1(); + if (is_dummy_pair(self->op(2), other->op(2))) { + *self = indexed(self_matrix.mul(other_matrix.transpose()), self->op(1), other->op(1)); + *other = _ex1(); return true; } // A_ji * B_jk = (Atrans*B)_ik - if (is_dummy_pair(self.op(1), other.op(1))) { - self = indexed(self_matrix.transpose().mul(other_matrix), self.op(2), other.op(2)); - other = _ex1(); + if (is_dummy_pair(self->op(1), other->op(1))) { + *self = indexed(self_matrix.transpose().mul(other_matrix), self->op(2), other->op(2)); + *other = _ex1(); return true; } // A_ji * B_kj = (B*A)_ki - if (is_dummy_pair(self.op(1), other.op(2))) { - self = indexed(other_matrix.mul(self_matrix), other.op(1), self.op(2)); - other = _ex1(); + if (is_dummy_pair(self->op(1), other->op(2))) { + *self = indexed(other_matrix.mul(self_matrix), other->op(1), self->op(2)); + *other = _ex1(); return true; } } diff --git a/ginac/matrix.h b/ginac/matrix.h index 7efbfed7..8c24dcf9 100644 --- a/ginac/matrix.h +++ b/ginac/matrix.h @@ -52,7 +52,7 @@ public: ex evalf(int level=0) const; ex subs(const lst & ls, const lst & lr) const; ex eval_indexed(const basic & i) const; - bool contract_with(ex & self, ex & other) const; + bool contract_with(exvector::iterator self, exvector::iterator other, exvector & v) const; protected: unsigned return_type(void) const { return return_types::noncommutative; }; // new virtual functions which can be overridden by derived classes diff --git a/ginac/normal.cpp b/ginac/normal.cpp index c1698ff7..669a59f8 100644 --- a/ginac/normal.cpp +++ b/ginac/normal.cpp @@ -240,8 +240,12 @@ static numeric lcmcoeff(const ex &e, const numeric &l) for (unsigned i=0; inops() == 3); + GINAC_ASSERT(is_ex_of_type(self->op(0), tensdelta)); // Try to contract first index - const idx *self_idx = &ex_to_idx(self.op(1)); - const idx *free_idx = &ex_to_idx(self.op(2)); + const idx *self_idx = &ex_to_idx(self->op(1)); + const idx *free_idx = &ex_to_idx(self->op(2)); bool first_index_tried = false; again: if (self_idx->is_symbolic()) { - for (int i=1; inops(); i++) { + const idx &other_idx = ex_to_idx(other->op(i)); if (is_dummy_pair(*self_idx, other_idx)) { // Contraction found, remove delta tensor and substitute // index in second object - self = _ex1(); - other = other.subs(other_idx == *free_idx); + *self = _ex1(); + *other = other->subs(other_idx == *free_idx); return true; } } @@ -280,8 +280,8 @@ again: if (!first_index_tried) { // No contraction with first index found, try second index - self_idx = &ex_to_idx(self.op(2)); - free_idx = &ex_to_idx(self.op(1)); + self_idx = &ex_to_idx(self->op(2)); + free_idx = &ex_to_idx(self->op(1)); first_index_tried = true; goto again; } @@ -290,33 +290,33 @@ again: } /** Contraction of an indexed metric tensor with something else. */ -bool tensmetric::contract_with(ex & self, ex & other) const +bool tensmetric::contract_with(exvector::iterator self, exvector::iterator other, exvector & v) const { - GINAC_ASSERT(is_ex_of_type(self, indexed)); - GINAC_ASSERT(is_ex_of_type(other, indexed)); - GINAC_ASSERT(self.nops() == 3); - GINAC_ASSERT(is_ex_of_type(self.op(0), tensmetric)); + GINAC_ASSERT(is_ex_of_type(*self, indexed)); + GINAC_ASSERT(is_ex_of_type(*other, indexed)); + GINAC_ASSERT(self->nops() == 3); + GINAC_ASSERT(is_ex_of_type(self->op(0), tensmetric)); // If contracting with the delta tensor, let the delta do it // (don't raise/lower delta indices) - if (is_ex_exactly_of_type(other.op(0), tensdelta)) + if (is_ex_exactly_of_type(other->op(0), tensdelta)) return false; // Try to contract first index - const idx *self_idx = &ex_to_idx(self.op(1)); - const idx *free_idx = &ex_to_idx(self.op(2)); + const idx *self_idx = &ex_to_idx(self->op(1)); + const idx *free_idx = &ex_to_idx(self->op(2)); bool first_index_tried = false; again: if (self_idx->is_symbolic()) { - for (int i=1; inops(); i++) { + const idx &other_idx = ex_to_idx(other->op(i)); if (is_dummy_pair(*self_idx, other_idx)) { // Contraction found, remove metric tensor and substitute // index in second object - self = _ex1(); - other = other.subs(other_idx == *free_idx); + *self = _ex1(); + *other = other->subs(other_idx == *free_idx); return true; } } @@ -325,8 +325,8 @@ again: if (!first_index_tried) { // No contraction with first index found, try second index - self_idx = &ex_to_idx(self.op(2)); - free_idx = &ex_to_idx(self.op(1)); + self_idx = &ex_to_idx(self->op(2)); + free_idx = &ex_to_idx(self->op(1)); first_index_tried = true; goto again; } diff --git a/ginac/tensor.h b/ginac/tensor.h index 220b69bc..9c4d2c36 100644 --- a/ginac/tensor.h +++ b/ginac/tensor.h @@ -55,7 +55,7 @@ class tensdelta : public tensor public: void print(std::ostream & os, unsigned upper_precedence=0) const; ex eval_indexed(const basic & i) const; - bool contract_with(ex & self, ex & other) const; + bool contract_with(exvector::iterator self, exvector::iterator other, exvector & v) const; }; @@ -70,7 +70,7 @@ class tensmetric : public tensor public: void print(std::ostream & os, unsigned upper_precedence=0) const; ex eval_indexed(const basic & i) const; - bool contract_with(ex & self, ex & other) const; + bool contract_with(exvector::iterator self, exvector::iterator other, exvector & v) const; }; -- 2.44.0