[GiNaC-devel] more patch
Chris Dams
Chris.Dams at mi.infn.it
Thu Oct 27 14:28:29 CEST 2005
Hi!
This time the changes are:
(1) The implementation that is currently in CVS disregards that it is
necessary to sort a container before applying algorithms such as
set_union.
(2) Don't rename indices if they both occur in one factor and no new dummy
indices are subsed into that factor.
(3) Introduced and documented a subs_options that switchs off dummy index
renaming.
(4) Use aforementioned option when symmetrizing expressions.
Actually, I am not really entirely certain if this subs_option is the
right way to go. On the one hand, the fact that constructing powers of
expressions with summed indices automatically yields renamed results,
suggest that we really want to consider index renaming as something that
happens automatically without the user having to know anything about it.
On the other hand, this may break existing code that wants to rename dummy
indices. The problem is probably quite rare, though. One needs something
like an ncmul in the mul to run into it. It is, of course, also possible
to first just substitute and then check whether this leads to inconsistent
indices and then, if necessary, to rename indices. This sounds
inefficient, though.
I would appreciate some input on this issue.
Best,
Chris
-------------- next part --------------
Index: doc/tutorial/ginac.texi
===================================================================
RCS file: /home/cvs/GiNaC/doc/tutorial/ginac.texi,v
retrieving revision 1.171
diff -c -r1.171 ginac.texi
*** doc/tutorial/ginac.texi 23 Sep 2005 07:02:33 -0000 1.171
--- doc/tutorial/ginac.texi 27 Oct 2005 11:07:55 -0000
***************
*** 4249,4261 ****
@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::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.
@code{subs()} performs syntactic substitution of any complete algebraic
object; it does not try to match sub-expressions as is demonstrated by the
--- 4249,4265 ----
@end example
The optional last argument to @code{subs()} is a combination of
! @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. 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
Index: ginac/color.h
===================================================================
RCS file: /home/cvs/GiNaC/ginac/color.h,v
retrieving revision 1.47
diff -c -r1.47 color.h
*** ginac/color.h 1 May 2005 18:23:26 -0000 1.47
--- ginac/color.h 27 Oct 2005 11:07:55 -0000
***************
*** 109,114 ****
--- 109,115 ----
// 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;
};
***************
*** 125,130 ****
--- 126,132 ----
// 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;
};
Index: ginac/expairseq.cpp
===================================================================
RCS file: /home/cvs/GiNaC/ginac/expairseq.cpp,v
retrieving revision 1.77
diff -c -r1.77 expairseq.cpp
*** ginac/expairseq.cpp 1 May 2005 18:23:26 -0000 1.77
--- ginac/expairseq.cpp 27 Oct 2005 11:07:55 -0000
***************
*** 34,39 ****
--- 34,40 ----
#include "archive.h"
#include "operators.h"
#include "utils.h"
+ #include "indexed.h"
#if EXPAIRSEQ_USE_HASHTAB
#include <cmath>
***************
*** 757,764 ****
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 EXPAIRSEQ_USE_HASHTAB
}
#endif // EXPAIRSEQ_USE_HASHTAB
--- 758,772 ----
construct_from_2_ex_via_exvector(lh,rh);
} else {
#endif // EXPAIRSEQ_USE_HASHTAB
! 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
***************
*** 1008,1020 ****
seq.reserve(v.size()+noperands-nexpairseqs);
// copy elements and split off numerical part
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()) {
seq.push_back(*cit_s);
++cit_s;
}
--- 1016,1042 ----
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 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;
}
***************
*** 1579,1584 ****
--- 1601,1668 ----
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.
*
***************
*** 1614,1635 ****
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());
// Copy parts of seq which are known not to have changed
! s->insert(s->begin(), seq.begin(), cit);
// Copy first changed element
! s->push_back(split_ex_to_pair(subsed_ex));
++cit;
// Copy rest
while (cit != last) {
! s->push_back(split_ex_to_pair(recombine_pair_to_ex(*cit).subs(m, options)));
++cit;
}
! return s;
}
++cit;
--- 1698,1724 ----
if (!are_ex_trivially_equal(orig_ex, subsed_ex)) {
// Something changed, copy seq, subs and return it
! safe_inserter s(*this, options & subs_options::no_index_renaming);
// Copy parts of seq which are known not to have changed
! for(epvector::const_iterator i=seq.begin(); i!=cit; ++i)
! s.insert_old_pair(*i);
// Copy first changed element
! s.insert_new_pair(split_ex_to_pair(subsed_ex), orig_ex);
++cit;
// Copy rest
while (cit != last) {
! 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.getseq();
}
++cit;
***************
*** 1645,1667 ****
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());
// Copy parts of seq which are known not to have changed
! s->insert(s->begin(), seq.begin(), cit);
// Copy first changed element
! s->push_back(combine_ex_with_coeff_to_pair(subsed_ex, cit->coeff));
++cit;
// Copy rest
while (cit != last) {
! s->push_back(combine_ex_with_coeff_to_pair(cit->rest.subs(m, options),
! cit->coeff));
++cit;
}
! return s;
}
++cit;
--- 1734,1760 ----
if (!are_ex_trivially_equal(cit->rest, subsed_ex)) {
// Something changed, copy seq, subs and return it
! safe_inserter s(*this, options & subs_options::no_index_renaming);
// Copy parts of seq which are known not to have changed
! for(epvector::const_iterator i=seq.begin(); i!=cit; ++i)
! s.insert_old_pair(*i);
// Copy first changed element
! s.insert_new_pair(combine_ex_with_coeff_to_pair(subsed_ex, cit->coeff), cit->rest);
++cit;
// Copy rest
while (cit != last) {
! ex orig_ex = cit->rest;
! 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(combine_ex_with_coeff_to_pair(subsed_ex, cit->coeff), orig_ex);
++cit;
}
! return s.getseq();
}
++cit;
Index: ginac/flags.h
===================================================================
RCS file: /home/cvs/GiNaC/ginac/flags.h,v
retrieving revision 1.36
diff -c -r1.36 flags.h
*** ginac/flags.h 22 Sep 2005 20:20:49 -0000 1.36
--- ginac/flags.h 27 Oct 2005 11:07:55 -0000
***************
*** 43,49 ****
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()
};
};
--- 43,50 ----
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()
! no_index_renaming = 0x0010
};
};
Index: ginac/indexed.cpp
===================================================================
RCS file: /home/cvs/GiNaC/ginac/indexed.cpp,v
retrieving revision 1.96
diff -c -r1.96 indexed.cpp
*** ginac/indexed.cpp 19 May 2005 14:10:40 -0000 1.96
--- ginac/indexed.cpp 27 Oct 2005 11:07:55 -0000
***************
*** 297,302 ****
--- 297,305 ----
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;
***************
*** 324,329 ****
--- 327,340 ----
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);
***************
*** 532,537 ****
--- 543,557 ----
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
***************
*** 540,549 ****
* @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)
{
! size_t global_size = global_dummy_indices.size(),
! local_size = local_dummy_indices.size();
// Any local dummy indices at all?
if (local_size == 0)
--- 560,569 ----
* @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 */
! template<class T> static ex rename_dummy_indices(const ex & e, exvector & global_dummy_indices, exvector & local_dummy_indices)
{
! 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)
***************
*** 557,563 ****
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()) {
global_dummy_indices.push_back(*it);
global_size++;
remaining--;
--- 577,583 ----
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 (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--;
***************
*** 575,585 ****
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));
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));
shaker_sort(global_syms.begin(), global_syms.end(), ex_is_less(), ex_swap());
// Remove common indices
--- 595,607 ----
exvector local_syms, global_syms;
local_syms.reserve(local_size);
global_syms.reserve(local_size);
! 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; 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
***************
*** 704,709 ****
--- 726,743 ----
}
}
+ 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)
***************
*** 864,882 ****
// 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;
! }
}
// Dummy index renaming
! r = rename_dummy_indices(r, dummy_indices, local_dummy_indices);
// Product of indexed object with a scalar?
if (is_exactly_a<mul>(r) && r.nops() == 2
--- 898,923 ----
// 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.
! 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<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
***************
*** 943,948 ****
--- 984,1000 ----
}
};
+ 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)
{
***************
*** 971,977 ****
}
// Rename the dummy indices
! return rename_dummy_indices(e_expanded, dummy_indices, local_dummy_indices);
}
// Simplification of sum = sum of simplifications, check consistency of
--- 1023,1032 ----
}
// Rename the 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
***************
*** 1015,1032 ****
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);
if (term_symm.is_zero())
continue;
terms.push_back(terminfo(term, term_symm));
--- 1070,1088 ----
if (num_terms_orig < 2 || dummy_indices.size() < 2)
return sum;
// 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);
! 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));
***************
*** 1318,1326 ****
return v;
}
! ex rename_dummy_indices_uniquely(const ex & a, const ex & b)
{
! exvector va = get_all_dummy_indices(a), vb = get_all_dummy_indices(b), 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;
--- 1374,1382 ----
return v;
}
! ex rename_dummy_indices_uniquely(const exvector & va, const exvector & vb, const ex & b)
{
! 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;
***************
*** 1330,1349 ****
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);
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()));
}
++ip;
}
return b.subs(lst(old_indices.begin(), old_indices.end()), lst(new_indices.begin(), new_indices.end()), subs_options::no_pattern);
}
}
ex expand_dummy_sum(const ex & e, bool subs_idx)
--- 1386,1424 ----
new_indices.reserve(2*common_indices.size());
exvector::const_iterator ip = common_indices.begin(), ipend = common_indices.end();
while (ip != ipend) {
! 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(ex_to<varidx>(newidx).toggle_variance());
}
++ip;
}
return b.subs(lst(old_indices.begin(), old_indices.end()), lst(new_indices.begin(), new_indices.end()), subs_options::no_pattern);
}
+ }
+
+ 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)
Index: ginac/indexed.h
===================================================================
RCS file: /home/cvs/GiNaC/ginac/indexed.h,v
retrieving revision 1.52
diff -c -r1.52 indexed.h
*** ginac/indexed.h 19 May 2005 14:10:40 -0000 1.52
--- ginac/indexed.h 27 Oct 2005 11:07:55 -0000
***************
*** 153,159 ****
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; }
ex expand(unsigned options = 0) const;
// new virtual functions which can be overridden by derived classes
--- 153,160 ----
ex derivative(const symbol & s) const;
ex thiscontainer(const exvector & v) const;
ex thiscontainer(std::auto_ptr<exvector> vp) const;
! 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
***************
*** 255,260 ****
--- 256,264 ----
/** 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
Index: ginac/symmetry.cpp
===================================================================
RCS file: /home/cvs/GiNaC/ginac/symmetry.cpp,v
retrieving revision 1.24
diff -c -r1.24 symmetry.cpp
*** ginac/symmetry.cpp 1 May 2005 18:23:27 -0000 1.24
--- ginac/symmetry.cpp 27 Oct 2005 11:07:55 -0000
***************
*** 427,433 ****
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);
if (asymmetric) {
memcpy(iv2, iv, num * sizeof(unsigned));
term *= permutation_sign(iv2, iv2 + num);
--- 427,433 ----
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|subs_options::no_index_renaming);
if (asymmetric) {
memcpy(iv2, iv, num * sizeof(unsigned));
term *= permutation_sign(iv2, iv2 + num);
***************
*** 468,474 ****
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);
}
return sum / num;
}
--- 468,474 ----
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|subs_options::no_index_renaming);
}
return sum / num;
}
Index: ginac/tensor.h
===================================================================
RCS file: /home/cvs/GiNaC/ginac/tensor.h,v
retrieving revision 1.23
diff -c -r1.23 tensor.h
*** ginac/tensor.h 1 May 2005 18:23:27 -0000 1.23
--- ginac/tensor.h 27 Oct 2005 11:07:55 -0000
***************
*** 65,70 ****
--- 65,71 ----
// 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;
};
***************
*** 84,89 ****
--- 85,91 ----
// 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;
};
***************
*** 106,111 ****
--- 108,114 ----
// 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;
***************
*** 153,158 ****
--- 156,162 ----
// 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;
More information about the GiNaC-devel
mailing list