X-Git-Url: https://ginac.de/ginac.git//ginac.git?a=blobdiff_plain;f=ginac%2Findexed.cpp;h=255fed5322fa1ed77c319a5ab5a567574a7359c5;hb=7064a931c4b5bea8c4a924fc9e9edabe5584b57a;hp=f5c53a6bb956e6606b125c80a248c840840cb809;hpb=e53c6d21b67ce4f5adafa7a33e703e202aa17ed0;p=ginac.git diff --git a/ginac/indexed.cpp b/ginac/indexed.cpp index f5c53a6b..255fed53 100644 --- a/ginac/indexed.cpp +++ b/ginac/indexed.cpp @@ -52,11 +52,7 @@ void indexed::copy(const indexed & other) symmetry = other.symmetry; } -void indexed::destroy(bool call_parent) -{ - if (call_parent) - inherited::destroy(call_parent); -} +DEFAULT_DESTROY(indexed) ////////// // other constructors @@ -159,7 +155,6 @@ indexed::indexed(symmetry_type symm, exvector * vp) : inherited(vp), symmetry(sy // archiving ////////// -/** Construct object from archive_node. */ indexed::indexed(const archive_node &n, const lst &sym_lst) : inherited(n, sym_lst) { debugmsg("indexed constructor from archive_node", LOGLEVEL_CONSTRUCT); @@ -168,19 +163,14 @@ indexed::indexed(const archive_node &n, const lst &sym_lst) : inherited(n, sym_l throw (std::runtime_error("unknown indexed symmetry type in archive")); } -/** Unarchive the object. */ -ex indexed::unarchive(const archive_node &n, const lst &sym_lst) -{ - return (new indexed(n, sym_lst))->setflag(status_flags::dynallocated); -} - -/** Archive the object. */ void indexed::archive(archive_node &n) const { inherited::archive(n); n.add_unsigned("symmetry", symmetry); } +DEFAULT_UNARCHIVE(indexed) + ////////// // functions overriding virtual functions from bases classes ////////// @@ -330,10 +320,20 @@ ex indexed::eval(int level) const if (level > 1) return indexed(symmetry, evalchildren(level)); + const ex &base = seq[0]; + // If the base object is 0, the whole object is 0 - if (seq[0].is_zero()) + if (base.is_zero()) return _ex0(); + // If the base object is a product, pull out the numeric factor + if (is_ex_exactly_of_type(base, mul) && is_ex_exactly_of_type(base.op(base.nops() - 1), numeric)) { + exvector v = seq; + ex f = ex_to_numeric(base.op(base.nops() - 1)); + v[0] = seq[0] / f; + return f * thisexprseq(v); + } + // Canonicalize indices according to the symmetry properties if (seq.size() > 2 && (symmetry != unknown && symmetry != mixed)) { exvector v = seq; @@ -347,7 +347,7 @@ ex indexed::eval(int level) const } // Let the class of the base object perform additional evaluations - return seq[0].bp->eval_indexed(*this); + return base.bp->eval_indexed(*this); } ex indexed::thisexprseq(const exvector & v) const @@ -445,48 +445,6 @@ void indexed::assert_all_indices_of_type_idx(void) const // global functions ////////// -/** Given a vector of indices, split them into two vectors, one containing - * the free indices, the other containing the dummy indices. */ -static void find_free_and_dummy(exvector::const_iterator it, exvector::const_iterator itend, exvector & out_free, exvector & out_dummy) -{ - out_free.clear(); - out_dummy.clear(); - - // No indices? Then do nothing - if (it == itend) - return; - - // Only one index? Then it is a free one if it's not numeric - if (itend - it == 1) { - if (ex_to_idx(*it).is_symbolic()) - out_free.push_back(*it); - return; - } - - // Sort index vector. This will cause dummy indices come to lie next - // to each other (because the sort order is defined to guarantee this). - exvector v(it, itend); - sort_index_vector(v); - - // Find dummy pairs and free indices - it = v.begin(); itend = v.end(); - exvector::const_iterator last = it++; - while (it != itend) { - if (is_dummy_pair(*it, *last)) { - out_dummy.push_back(*last); - it++; - if (it == itend) - return; - } else { - if (!it->is_equal(*last) && ex_to_idx(*last).is_symbolic()) - out_free.push_back(*last); - } - last = it++; - } - if (ex_to_idx(*last).is_symbolic()) - out_free.push_back(*last); -} - /** Check whether two sorted index vectors are consistent (i.e. equal). */ static bool indices_consistent(const exvector & v1, const exvector & v2) { @@ -505,6 +463,12 @@ static bool indices_consistent(const exvector & v1, const exvector & v2) return true; } +exvector indexed::get_indices(void) const +{ + GINAC_ASSERT(seq.size() >= 1); + return exvector(seq.begin() + 1, seq.end()); +} + exvector indexed::get_dummy_indices(void) const { exvector free_indices, dummy_indices; @@ -512,6 +476,16 @@ exvector indexed::get_dummy_indices(void) const return dummy_indices; } +exvector indexed::get_dummy_indices(const indexed & other) const +{ + exvector indices = get_free_indices(); + exvector other_indices = other.get_free_indices(); + indices.insert(indices.end(), other_indices.begin(), other_indices.end()); + exvector dummy_indices; + find_dummy_indices(indices, dummy_indices); + return dummy_indices; +} + exvector indexed::get_free_indices(void) const { exvector free_indices, dummy_indices; @@ -545,7 +519,7 @@ exvector mul::get_free_indices(void) const // And remove the dummy indices exvector free_indices, dummy_indices; - find_free_and_dummy(un.begin(), un.end(), free_indices, dummy_indices); + find_free_and_dummy(un, free_indices, dummy_indices); return free_indices; } @@ -560,7 +534,7 @@ exvector ncmul::get_free_indices(void) const // And remove the dummy indices exvector free_indices, dummy_indices; - find_free_and_dummy(un.begin(), un.end(), free_indices, dummy_indices); + find_free_and_dummy(un, free_indices, dummy_indices); return free_indices; } @@ -596,7 +570,7 @@ ex simplify_indexed_product(const ex & e, exvector & free_indices, const scalar_ } else if (is_ex_exactly_of_type(f, ncmul)) { // Noncommutative factor found, split it as well non_commutative = true; // everything becomes noncommutative, ncmul will sort out the commutative factors later - for (int j=0; j 1) { + free_indices.clear(); + return _ex0(); + } + // Try to contract the first one with the second one bool contracted = it1->op(0).bp->contract_with(it1, it2, v); if (!contracted) { @@ -667,15 +659,23 @@ try_again: } it1++; } - find_free_and_dummy(un.begin(), un.end(), free_indices, dummy_indices); + find_free_and_dummy(un, 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. */ @@ -695,15 +695,27 @@ ex simplify_indexed(const ex & e, exvector & free_indices, const scalar_products // Simplification of sum = sum of simplifications, check consistency of // free indices in each term if (is_ex_exactly_of_type(e_expanded, add)) { + bool first = true; ex sum = _ex0(); + free_indices.clear(); for (unsigned i=0; iadd_indexed(sum, term); + else + sum += term; + } + } } return sum;