- 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

diff --git a/ginac/basic.cpp b/ginac/basic.cpp
index 2414c92971e99c07a8f4e91377d9afd9829d28ab..a7645ac5cba460d0622dcafc6df5b7223aaca22d 100644 (file)
--- 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 842417cad1ff6a5a780cc491b2a8658dfe65413e..ee2bef237d0aeb569efaa55f6362def40078d216 100644 (file)
--- 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 ea4266cba5ac8e0c5e93d26d58e998060bd61193..c3fe1d9a26e6d49ca311b0e793913d6050ef402f 100644 (file)
--- 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 a63275f54338e7688065a79176091e36902503a1..1a5d6953c3d994365d9d6b05a60ebf351119e7bd 100644 (file)
--- 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 7efbfed7cfc4c5f255f57f9ead86e889003d7f6e..8c24dcf9043a423c21521841928327eab22586a4 100644 (file)
--- 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 c1698ff79af62e63dbd30abd36c3499c7a0a8705..669a59f8de12d0a7cfb3d3eba694eaa9703d275d 100644 (file)
--- 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; i<e.nops(); i++)
                        c *= lcmcoeff(e.op(i), _num1());
                return lcm(c, l);
-       } else if (is_ex_exactly_of_type(e, power))
-               return pow(lcmcoeff(e.op(0), l), ex_to_numeric(e.op(1)));
+       } else if (is_ex_exactly_of_type(e, power)) {
+               if (is_ex_exactly_of_type(e.op(0), symbol))
+                       return l;
+               else
+                       return pow(lcmcoeff(e.op(0), l), ex_to_numeric(e.op(1)));
+       }
        return l;
 }
 
@@ -280,7 +284,10 @@ static ex multiply_lcm(const ex &e, const numeric &lcm)
                        c += multiply_lcm(e.op(i), lcm);
                return c;
        } else if (is_ex_exactly_of_type(e, power)) {
-               return pow(multiply_lcm(e.op(0), lcm.power(ex_to_numeric(e.op(1)).inverse())), e.op(1));
+               if (is_ex_exactly_of_type(e.op(0), symbol))
+                       return e * lcm;
+               else
+                       return pow(multiply_lcm(e.op(0), lcm.power(ex_to_numeric(e.op(1)).inverse())), e.op(1));
        } else
                return e * lcm;
 }

diff --git a/ginac/tensor.cpp b/ginac/tensor.cpp
index 28b520f1bc57242c716b1d8369f17e0b3b8e34e7..9d5feb63c4051dbe5c9d65ee3a0398fd5ada7c96 100644 (file)
--- a/ginac/tensor.cpp
+++ b/ginac/tensor.cpp
@@ -250,28 +250,28 @@ ex minkmetric::eval_indexed(const basic & i) const
 }
 
 /** Contraction of an indexed delta tensor with something else. */
-bool tensdelta::contract_with(ex & self, ex & other) const
+bool tensdelta::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), tensdelta));
+       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), 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; i<other.nops(); i++) {
-                       const idx &other_idx = ex_to_idx(other.op(i));
+               for (int i=1; i<other->nops(); 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; i<other.nops(); i++) {
-                       const idx &other_idx = ex_to_idx(other.op(i));
+               for (int i=1; i<other->nops(); 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 220b69bcb56e4a493103497160788373f65ada62..9c4d2c36137a8a6dd7b049bf3d055b639b7f4398 100644 (file)
--- 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;
 };