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	24 Oct 2005 15:36:19 -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,1041 ----
  	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);
! 				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 ****
--- 1600,1654 ----
  	return std::auto_ptr<epvector>(0); // signalling nothing has changed
  }
  
+ class safe_inserter
+ {
+ 	public:
+ 		safe_inserter(const ex&);
+ 		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);
+ 	private:
+ 		std::auto_ptr<epvector> epv;
+ 		bool dodummies;
+ 		exvector dummy_indices;
+ 		void update_dummy_indices(const exvector&);
+ };
+ 
+ safe_inserter::safe_inserter(const ex&e):epv(new epvector)
+ {
+ 	epv->reserve(e.nops());
+ 	dodummies=is_a<mul>(e);
+ 	if(dodummies)
+ 		dummy_indices = get_all_dummy_indices(e);
+ }
+ 
+ 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)
+ {
+ 	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;
+ 	}
+ 	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;
--- 1684,1710 ----
  			if (!are_ex_trivially_equal(orig_ex, subsed_ex)) {
  
  				// Something changed, copy seq, subs and return it
! 				safe_inserter s(*this);
  
  				// 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));
  				++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));
  					++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;
--- 1720,1746 ----
  			if (!are_ex_trivially_equal(cit->rest, subsed_ex)) {
  			
  				// Something changed, copy seq, subs and return it
! 				safe_inserter s(*this);
  
  				// 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;
  
  				// 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));
  					++cit;
  				}
! 				return s.getseq();
  			}
  
  			++cit;
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	24 Oct 2005 15:36:19 -0000
***************
*** 297,302 ****
--- 297,305 ----
  		return f * thiscontainer(v);
  	}
  
+ 	if(seq.size()==1)
+ 		return base;
+ 
  	// Canonicalize indices according to the symmetry properties
  	if (seq.size() > 2) {
  		exvector v = seq;
***************
*** 532,537 ****
--- 535,549 ----
  	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)
--- 552,561 ----
   *  @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--;
--- 569,575 ----
  		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
--- 587,599 ----
  	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 ****
--- 718,735 ----
  	}
  }
  
+ 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
--- 890,915 ----
  	// 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 ****
--- 976,992 ----
  	}
  };
  
+ 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
--- 1015,1024 ----
  		}
  
  		// 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));
--- 1062,1080 ----
  		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;
--- 1366,1374 ----
  	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)
--- 1378,1414 ----
  		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);
+ 	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	24 Oct 2005 15:36:19 -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 { return op(0).return_type(); }
! 	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/operators.cpp
===================================================================
RCS file: /home/cvs/GiNaC/ginac/operators.cpp,v
retrieving revision 1.33
diff -c -r1.33 operators.cpp
*** ginac/operators.cpp	29 Jun 2005 14:09:24 -0000	1.33
--- ginac/operators.cpp	24 Oct 2005 15:36:19 -0000
***************
*** 47,52 ****
--- 47,55 ----
  	// Check if we are constructing a mul object or a ncmul object.  Due to
  	// ncmul::eval()'s rule to pull out commutative elements we need to check
  	// only one of the elements.
+ 	std::cout << "exmul: multiplying " << lh << " and " << rh << std::endl;
+ 	std::cout << "lhs is commutative " << int(lh.return_type()==return_types::commutative) << std::endl;
+ 	std::cout << "rhs is commutative " << int(rh.return_type()==return_types::commutative) << std::endl;
  	if (rh.return_type()==return_types::commutative ||
  	    lh.return_type()==return_types::commutative)
  		return (new mul(lh,rh))->setflag(status_flags::dynallocated);
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	24 Oct 2005 15:36:19 -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;