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	30 Sep 2005 14:43:59 -0000
***************
*** 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--;
--- 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(idx_is_equal_ignore_dim(), *it)) == global_dummy_indices.end()) {
  				global_dummy_indices.push_back(*it);
  				global_size++;
  				remaining--;
***************
*** 704,709 ****
--- 704,721 ----
  	}
  }
  
+ 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,878 ****
  	// 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
--- 876,895 ----
  	// 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
***************
*** 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));
--- 1032,1045 ----
  		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);
! 			ex term_symm = idx_symmetrization<idx>(term, dummy_indices);
! 			term_symm = idx_symmetrization<varidx>(term_symm, dummy_indices);
! 			term_symm = idx_symmetrization<spinidx>(term_symm, dummy_indices);
  			if (term_symm.is_zero())
  				continue;
  			terms.push_back(terminfo(term, term_symm));