[GiNaC-devel] patch for simplify_indexed.

[Chris Dams]

Dear developers,

When simplify_indexed() is called two symmetrizations w.r.t dummy indices
are performed. In these two cases symmetrization w.r.t. all dummy indices
is done. However, it is also possible to do symmetrization with respect to
ordinary idx-es, varidx-es and spinidx-es seperately. The latter is
potentially much faster. E.g. Today I had an example where some terms
contained 5 ordinary dummy indices and 2 dummy indices with variance. In
this case a factor of 21 fewer terms are generated by the latter option.  
(I.e. 7!/5!/2!=21). Therefore symmetrizing seperately can be much faster.  
A patch is attached.

Best wishes,
Chris
-------------- next part --------------
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));