can carry indices.
@end example
The optional last argument to @code{subs()} is a combination of
-@code{subs_options} flags. There are two options available:
+@code{subs_options} flags. There are three options available:
@code{subs_options::no_pattern} disables pattern matching, which makes
large @code{subs()} operations significantly faster if you are not using
patterns. The second option, @code{subs_options::algebraic} enables
algebraic substitutions in products and powers.
@ref{Pattern Matching and Advanced Substitutions}, for more information
-about patterns and algebraic substitutions.
+about patterns and algebraic substitutions. The third option,
+@code{subs_options::no_index_renaming} disables the feature that dummy
+indices are renamed if the subsitution could give a result in which a
+dummy index occurs more than two times. This is sometimes necessary if
+you want to use @code{subs()} to rename your dummy indices.
@code{subs()} performs syntactic substitution of any complete algebraic
object; it does not try to match sub-expressions as is demonstrated by the
// non-virtual functions in this class
protected:
+ unsigned return_type() const { return return_types::commutative; }
void do_print(const print_context & c, unsigned level) const;
void do_print_latex(const print_latex & c, unsigned level) const;
};
// non-virtual functions in this class
protected:
+ unsigned return_type() const { return return_types::commutative; }
void do_print(const print_context & c, unsigned level) const;
void do_print_latex(const print_latex & c, unsigned level) const;
};
#include "archive.h"
#include "operators.h"
#include "utils.h"
+#include "indexed.h"
#if EXPAIRSEQ_USE_HASHTAB
#include <cmath>
construct_from_2_ex_via_exvector(lh,rh);
} else {
#endif // EXPAIRSEQ_USE_HASHTAB
- construct_from_2_expairseq(ex_to<expairseq>(lh),
- ex_to<expairseq>(rh));
+ if(is_a<mul>(lh))
+ {
+ ex newrh=rename_dummy_indices_uniquely(lh, rh);
+ construct_from_2_expairseq(ex_to<expairseq>(lh),
+ ex_to<expairseq>(newrh));
+ }
+ else
+ construct_from_2_expairseq(ex_to<expairseq>(lh),
+ ex_to<expairseq>(rh));
#if EXPAIRSEQ_USE_HASHTAB
}
#endif // EXPAIRSEQ_USE_HASHTAB
seq.reserve(v.size()+noperands-nexpairseqs);
// copy elements and split off numerical part
+ exvector dummy_indices;
cit = v.begin();
while (cit!=v.end()) {
if (ex_to<basic>(*cit).tinfo()==this->tinfo()) {
- const expairseq &subseqref = ex_to<expairseq>(*cit);
- combine_overall_coeff(subseqref.overall_coeff);
- epvector::const_iterator cit_s = subseqref.seq.begin();
- while (cit_s!=subseqref.seq.end()) {
+ const expairseq *subseqref;
+ ex newfactor;
+ if(is_a<mul>(*cit))
+ {
+ exvector dummies_of_factor = get_all_dummy_indices(*cit);
+ sort(dummies_of_factor.begin(), dummies_of_factor.end(), ex_is_less());
+ newfactor = rename_dummy_indices_uniquely(dummy_indices, dummies_of_factor, *cit);
+ subseqref = &(ex_to<expairseq>(newfactor));
+ exvector new_dummy_indices;
+ set_union(dummy_indices.begin(), dummy_indices.end(), dummies_of_factor.begin(), dummies_of_factor.end(), std::back_insert_iterator<exvector>(new_dummy_indices), ex_is_less());
+ dummy_indices.swap(new_dummy_indices);
+ }
+ else
+ subseqref = &ex_to<expairseq>(*cit);
+ combine_overall_coeff(subseqref->overall_coeff);
+ epvector::const_iterator cit_s = subseqref->seq.begin();
+ while (cit_s!=subseqref->seq.end()) {
seq.push_back(*cit_s);
++cit_s;
}
return std::auto_ptr<epvector>(0); // signalling nothing has changed
}
+class safe_inserter
+{
+ public:
+ safe_inserter(const ex&, const bool disable_renaming=false);
+ std::auto_ptr<epvector> getseq(){return epv;}
+ void insert_old_pair(const expair &p)
+ {
+ epv->push_back(p);
+ }
+ void insert_new_pair(const expair &p, const ex &orig_ex);
+ private:
+ std::auto_ptr<epvector> epv;
+ bool dodummies;
+ exvector dummy_indices;
+ void update_dummy_indices(const exvector&);
+};
+
+safe_inserter::safe_inserter(const ex&e, const bool disable_renaming)
+ :epv(new epvector)
+{
+ epv->reserve(e.nops());
+ dodummies=is_a<mul>(e);
+ if(disable_renaming)
+ dodummies=false;
+ if(dodummies) {
+ dummy_indices = get_all_dummy_indices(e);
+ sort(dummy_indices.begin(), dummy_indices.end(), ex_is_less());
+ }
+}
+
+void safe_inserter::update_dummy_indices(const exvector &v)
+{
+ exvector new_dummy_indices;
+ set_union(dummy_indices.begin(), dummy_indices.end(), v.begin(), v.end(),
+ std::back_insert_iterator<exvector>(new_dummy_indices), ex_is_less());
+ dummy_indices.swap(new_dummy_indices);
+}
+
+void safe_inserter::insert_new_pair(const expair &p, const ex &orig_ex)
+{
+ if(!dodummies) {
+ epv->push_back(p);
+ return;
+ }
+ exvector dummies_of_factor = get_all_dummy_indices(p.rest);
+ if(dummies_of_factor.size() == 0) {
+ epv->push_back(p);
+ return;
+ }
+ sort(dummies_of_factor.begin(), dummies_of_factor.end(), ex_is_less());
+ exvector dummies_of_orig_ex = get_all_dummy_indices(orig_ex);
+ sort(dummies_of_orig_ex.begin(), dummies_of_orig_ex.end(), ex_is_less());
+ exvector new_dummy_indices;
+ new_dummy_indices.reserve(dummy_indices.size());
+ set_difference(dummy_indices.begin(), dummy_indices.end(), dummies_of_orig_ex.begin(), dummies_of_orig_ex.end(),
+ std::back_insert_iterator<exvector>(new_dummy_indices), ex_is_less());
+ dummy_indices.swap(new_dummy_indices);
+ ex newfactor = rename_dummy_indices_uniquely(dummy_indices, dummies_of_factor, p.rest);
+ update_dummy_indices(dummies_of_factor);
+ epv -> push_back(expair(newfactor, p.coeff));
+}
/** Member-wise substitute in this sequence.
*
if (!are_ex_trivially_equal(orig_ex, subsed_ex)) {
// Something changed, copy seq, subs and return it
- std::auto_ptr<epvector> s(new epvector);
- s->reserve(seq.size());
+ safe_inserter s(*this, options & subs_options::no_index_renaming);
// Copy parts of seq which are known not to have changed
- s->insert(s->begin(), seq.begin(), cit);
+ for(epvector::const_iterator i=seq.begin(); i!=cit; ++i)
+ s.insert_old_pair(*i);
// Copy first changed element
- s->push_back(split_ex_to_pair(subsed_ex));
+ s.insert_new_pair(split_ex_to_pair(subsed_ex), orig_ex);
++cit;
// Copy rest
while (cit != last) {
- s->push_back(split_ex_to_pair(recombine_pair_to_ex(*cit).subs(m, options)));
+ ex orig_ex = recombine_pair_to_ex(*cit);
+ ex subsed_ex = orig_ex.subs(m, options);
+ if(are_ex_trivially_equal(orig_ex, subsed_ex))
+ s.insert_old_pair(*cit);
+ else
+ s.insert_new_pair(split_ex_to_pair(subsed_ex), orig_ex);
++cit;
}
- return s;
+ return s.getseq();
}
++cit;
if (!are_ex_trivially_equal(cit->rest, subsed_ex)) {
// Something changed, copy seq, subs and return it
- std::auto_ptr<epvector> s(new epvector);
- s->reserve(seq.size());
+ safe_inserter s(*this, options & subs_options::no_index_renaming);
// Copy parts of seq which are known not to have changed
- s->insert(s->begin(), seq.begin(), cit);
+ for(epvector::const_iterator i=seq.begin(); i!=cit; ++i)
+ s.insert_old_pair(*i);
// Copy first changed element
- s->push_back(combine_ex_with_coeff_to_pair(subsed_ex, cit->coeff));
+ s.insert_new_pair(combine_ex_with_coeff_to_pair(subsed_ex, cit->coeff), cit->rest);
++cit;
// Copy rest
while (cit != last) {
- s->push_back(combine_ex_with_coeff_to_pair(cit->rest.subs(m, options),
- cit->coeff));
+ const ex &orig_ex = cit->rest;
+ const ex &subsed_ex = cit->rest.subs(m, options);
+ if(are_ex_trivially_equal(orig_ex, subsed_ex))
+ s.insert_old_pair(*cit);
+ else
+ s.insert_new_pair(combine_ex_with_coeff_to_pair(subsed_ex, cit->coeff), orig_ex);
++cit;
}
- return s;
+ return s.getseq();
}
++cit;
algebraic = 0x0002, ///< enable algebraic substitutions
subs_algebraic = 0x0002, // for backwards compatibility
pattern_is_product = 0x0004, ///< used internally by expairseq::subschildren()
- pattern_is_not_product = 0x0008 ///< used internally by expairseq::subschildren()
+ pattern_is_not_product = 0x0008, ///< used internally by expairseq::subschildren()
+ no_index_renaming = 0x0010
};
};
return f * thiscontainer(v);
}
+ if(this->tinfo()==TINFO_indexed && seq.size()==1)
+ return base;
+
// Canonicalize indices according to the symmetry properties
if (seq.size() > 2) {
exvector v = seq;
return indexed(ex_to<symmetry>(symtree), vp);
}
+unsigned indexed::return_type() const
+{
+ if(is_a<matrix>(op(0)))
+ return return_types::commutative;
+ else
+ return op(0).return_type();
+}
+
ex indexed::expand(unsigned options) const
{
GINAC_ASSERT(seq.size() > 0);
return f.get_free_indices();
}
+template<class T> size_t number_of_type(const exvector&v)
+{
+ size_t number = 0;
+ for(exvector::const_iterator i=v.begin(); i!=v.end(); ++i)
+ if(is_exactly_a<T>(*i))
+ ++number;
+ return number;
+}
+
/** Rename dummy indices in an expression.
*
* @param e Expression to work on
* @param global_dummy_indices The set of dummy indices that have appeared
* before and which we would like to use in "e", too. This gets updated
* by the function */
-static ex rename_dummy_indices(const ex & e, exvector & global_dummy_indices, exvector & local_dummy_indices)
+template<class T> static ex rename_dummy_indices(const ex & e, exvector & global_dummy_indices, exvector & local_dummy_indices)
{
- size_t global_size = global_dummy_indices.size(),
- local_size = local_dummy_indices.size();
+ size_t global_size = number_of_type<T>(global_dummy_indices),
+ local_size = number_of_type<T>(local_dummy_indices);
// Any local dummy indices at all?
if (local_size == 0)
int remaining = local_size - global_size;
exvector::const_iterator it = local_dummy_indices.begin(), itend = local_dummy_indices.end();
while (it != itend && remaining > 0) {
- if (find_if(global_dummy_indices.begin(), global_dummy_indices.end(), bind2nd(op0_is_equal(), *it)) == global_dummy_indices.end()) {
+ if (is_exactly_a<T>(*it) && find_if(global_dummy_indices.begin(), global_dummy_indices.end(), bind2nd(idx_is_equal_ignore_dim(), *it)) == global_dummy_indices.end()) {
global_dummy_indices.push_back(*it);
global_size++;
remaining--;
exvector local_syms, global_syms;
local_syms.reserve(local_size);
global_syms.reserve(local_size);
- for (size_t i=0; i<local_size; i++)
- local_syms.push_back(local_dummy_indices[i].op(0));
+ for (size_t i=0; local_syms.size()!=local_size; i++)
+ if(is_exactly_a<T>(local_dummy_indices[i]))
+ local_syms.push_back(local_dummy_indices[i].op(0));
shaker_sort(local_syms.begin(), local_syms.end(), ex_is_less(), ex_swap());
- for (size_t i=0; i<local_size; i++) // don't use more global symbols than necessary
- global_syms.push_back(global_dummy_indices[i].op(0));
+ for (size_t i=0; global_syms.size()!=local_size; i++) // don't use more global symbols than necessary
+ if(is_exactly_a<T>(global_dummy_indices[i]))
+ global_syms.push_back(global_dummy_indices[i].op(0));
shaker_sort(global_syms.begin(), global_syms.end(), ex_is_less(), ex_swap());
// Remove common indices
}
}
+template<class T> ex idx_symmetrization(const ex& r,const exvector& local_dummy_indices)
+{ exvector dummy_syms;
+ dummy_syms.reserve(r.nops());
+ for (exvector::const_iterator it = local_dummy_indices.begin(); it != local_dummy_indices.end(); ++it)
+ if(is_exactly_a<T>(*it))
+ dummy_syms.push_back(it->op(0));
+ if(dummy_syms.size() < 2)
+ return r;
+ ex q=symmetrize(r, dummy_syms);
+ return q;
+}
+
/** Simplify product of indexed expressions (commutative, noncommutative and
* simple squares), return list of free indices. */
ex simplify_indexed_product(const ex & e, exvector & free_indices, exvector & dummy_indices, const scalar_products & sp)
// The result should be symmetric with respect to exchange of dummy
// indices, so if the symmetrization vanishes, the whole expression is
// zero. This detects things like eps.i.j.k * p.j * p.k = 0.
- if (local_dummy_indices.size() >= 2) {
- exvector dummy_syms;
- dummy_syms.reserve(local_dummy_indices.size());
- for (exvector::const_iterator it = local_dummy_indices.begin(); it != local_dummy_indices.end(); ++it)
- dummy_syms.push_back(it->op(0));
- if (symmetrize(r, dummy_syms).is_zero()) {
- free_indices.clear();
- return _ex0;
- }
+ ex q = idx_symmetrization<idx>(r, local_dummy_indices);
+ if (q.is_zero()) {
+ free_indices.clear();
+ return _ex0;
+ }
+ q = idx_symmetrization<varidx>(q, local_dummy_indices);
+ if (q.is_zero()) {
+ free_indices.clear();
+ return _ex0;
+ }
+ q = idx_symmetrization<spinidx>(q, local_dummy_indices);
+ if (q.is_zero()) {
+ free_indices.clear();
+ return _ex0;
}
// Dummy index renaming
- r = rename_dummy_indices(r, dummy_indices, local_dummy_indices);
+ r = rename_dummy_indices<idx>(r, dummy_indices, local_dummy_indices);
+ r = rename_dummy_indices<varidx>(r, dummy_indices, local_dummy_indices);
+ r = rename_dummy_indices<spinidx>(r, dummy_indices, local_dummy_indices);
// Product of indexed object with a scalar?
if (is_exactly_a<mul>(r) && r.nops() == 2
}
};
+bool hasindex(const ex &x, const ex &sym)
+{
+ if(is_a<idx>(x) && x.op(0)==sym)
+ return true;
+ else
+ for(size_t i=0; i<x.nops(); ++i)
+ if(hasindex(x.op(i), sym))
+ return true;
+ return false;
+}
+
/** Simplify indexed expression, return list of free indices. */
ex simplify_indexed(const ex & e, exvector & free_indices, exvector & dummy_indices, const scalar_products & sp)
{
}
// Rename the dummy indices
- return rename_dummy_indices(e_expanded, dummy_indices, local_dummy_indices);
+ e_expanded = rename_dummy_indices<idx>(e_expanded, dummy_indices, local_dummy_indices);
+ e_expanded = rename_dummy_indices<varidx>(e_expanded, dummy_indices, local_dummy_indices);
+ e_expanded = rename_dummy_indices<spinidx>(e_expanded, dummy_indices, local_dummy_indices);
+ return e_expanded;
}
// Simplification of sum = sum of simplifications, check consistency of
if (num_terms_orig < 2 || dummy_indices.size() < 2)
return sum;
- // Yes, construct vector of all dummy index symbols
- exvector dummy_syms;
- dummy_syms.reserve(dummy_indices.size());
- for (exvector::const_iterator it = dummy_indices.begin(); it != dummy_indices.end(); ++it)
- dummy_syms.push_back(it->op(0));
-
// Chop the sum into terms and symmetrize each one over the dummy
// indices
std::vector<terminfo> terms;
for (size_t i=0; i<sum.nops(); i++) {
const ex & term = sum.op(i);
- ex term_symm = symmetrize(term, dummy_syms);
+ exvector dummy_indices_of_term;
+ dummy_indices_of_term.reserve(dummy_indices.size());
+ for(exvector::iterator i=dummy_indices.begin(); i!=dummy_indices.end(); ++i)
+ if(hasindex(term,i->op(0)))
+ dummy_indices_of_term.push_back(*i);
+ ex term_symm = idx_symmetrization<idx>(term, dummy_indices_of_term);
+ term_symm = idx_symmetrization<varidx>(term_symm, dummy_indices_of_term);
+ term_symm = idx_symmetrization<spinidx>(term_symm, dummy_indices_of_term);
if (term_symm.is_zero())
continue;
terms.push_back(terminfo(term, term_symm));
return v;
}
-ex rename_dummy_indices_uniquely(const ex & a, const ex & b)
+ex rename_dummy_indices_uniquely(const exvector & va, const exvector & vb, const ex & b)
{
- exvector va = get_all_dummy_indices(a), vb = get_all_dummy_indices(b), common_indices;
+ exvector common_indices;
set_intersection(va.begin(), va.end(), vb.begin(), vb.end(), std::back_insert_iterator<exvector>(common_indices), ex_is_less());
if (common_indices.empty()) {
return b;
new_indices.reserve(2*common_indices.size());
exvector::const_iterator ip = common_indices.begin(), ipend = common_indices.end();
while (ip != ipend) {
- if (is_a<varidx>(*ip)) {
- varidx mu((new symbol)->setflag(status_flags::dynallocated), ex_to<varidx>(*ip).get_dim(), ex_to<varidx>(*ip).is_covariant());
- old_indices.push_back(*ip);
- new_indices.push_back(mu);
+ ex newsym=(new symbol)->setflag(status_flags::dynallocated);
+ ex newidx;
+ if(is_exactly_a<spinidx>(*ip))
+ newidx = (new spinidx(newsym, ex_to<spinidx>(*ip).get_dim(),
+ ex_to<spinidx>(*ip).is_covariant(),
+ ex_to<spinidx>(*ip).is_dotted()))
+ -> setflag(status_flags::dynallocated);
+ else if (is_exactly_a<varidx>(*ip))
+ newidx = (new varidx(newsym, ex_to<varidx>(*ip).get_dim(),
+ ex_to<varidx>(*ip).is_covariant()))
+ -> setflag(status_flags::dynallocated);
+ else
+ newidx = (new idx(newsym, ex_to<idx>(*ip).get_dim()))
+ -> setflag(status_flags::dynallocated);
+ old_indices.push_back(*ip);
+ new_indices.push_back(newidx);
+ if(is_a<varidx>(*ip)) {
old_indices.push_back(ex_to<varidx>(*ip).toggle_variance());
- new_indices.push_back(mu.toggle_variance());
- } else {
- old_indices.push_back(*ip);
- new_indices.push_back(idx((new symbol)->setflag(status_flags::dynallocated), ex_to<varidx>(*ip).get_dim()));
+ new_indices.push_back(ex_to<varidx>(newidx).toggle_variance());
}
++ip;
}
}
}
+ex rename_dummy_indices_uniquely(const ex & a, const ex & b)
+{
+ exvector va = get_all_dummy_indices(a);
+ exvector vb = get_all_dummy_indices(b);
+ sort(va.begin(), va.end(), ex_is_less());
+ sort(vb.begin(), vb.end(), ex_is_less());
+ return rename_dummy_indices_uniquely(va, vb, b);
+}
+
ex expand_dummy_sum(const ex & e, bool subs_idx)
{
ex e_expanded = e.expand();
ex derivative(const symbol & s) const;
ex thiscontainer(const exvector & v) const;
ex thiscontainer(std::auto_ptr<exvector> vp) const;
- unsigned return_type() const { return return_types::commutative; }
+ unsigned return_type() const;
+ unsigned return_type_tinfo() const { return op(0).return_type_tinfo(); }
ex expand(unsigned options = 0) const;
// new virtual functions which can be overridden by derived classes
/** Returns b with all dummy indices, which are common with a, renamed */
ex rename_dummy_indices_uniquely(const ex & a, const ex & b);
+/** Same as above, where va and vb contain the indices of a and b and are sorted */
+ex rename_dummy_indices_uniquely(const exvector & va, const exvector & vb, const ex & b);
+
/** This function returns the given expression with expanded sums
* for all dummy index summations, where the dimensionality of
* the dummy index is numeric.
lst new_lst;
for (unsigned i=0; i<num; i++)
new_lst.append(orig_lst.op(iv[i]));
- ex term = e.subs(orig_lst, new_lst, subs_options::no_pattern);
+ ex term = e.subs(orig_lst, new_lst, subs_options::no_pattern|subs_options::no_index_renaming);
if (asymmetric) {
memcpy(iv2, iv, num * sizeof(unsigned));
term *= permutation_sign(iv2, iv2 + num);
for (unsigned i=0; i<num-1; i++) {
ex perm = new_lst.op(0);
new_lst.remove_first().append(perm);
- sum += e.subs(orig_lst, new_lst, subs_options::no_pattern);
+ sum += e.subs(orig_lst, new_lst, subs_options::no_pattern|subs_options::no_index_renaming);
}
return sum / num;
}
// non-virtual functions in this class
protected:
+ unsigned return_type() const { return return_types::commutative; }
void do_print(const print_context & c, unsigned level) const;
void do_print_latex(const print_latex & c, unsigned level) const;
};
// non-virtual functions in this class
protected:
+ unsigned return_type() const { return return_types::commutative; }
void do_print(const print_context & c, unsigned level) const;
};
// non-virtual functions in this class
protected:
+ unsigned return_type() const { return return_types::commutative; }
void do_print(const print_context & c, unsigned level) const;
void do_print_latex(const print_latex & c, unsigned level) const;
// non-virtual functions in this class
protected:
+ unsigned return_type() const { return return_types::commutative; }
void do_print(const print_context & c, unsigned level) const;
void do_print_latex(const print_latex & c, unsigned level) const;