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	3 Oct 2005 09:33:21 -0000
***************
*** 532,537 ****
--- 532,546 ----
  	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)
--- 549,558 ----
   *  @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--;
--- 566,572 ----
  		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_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
--- 584,596 ----
  	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 ****
--- 715,732 ----
  	}
  }
  
+ 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
--- 887,912 ----
  	// 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 ****
--- 973,989 ----
  	}
  };
  
+ 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
--- 1012,1021 ----
  		}
  
  		// 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));
--- 1059,1077 ----
  		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));