problem with simplify_indexed() plus patch

[Chris Dams]

Hello everybody,

I have a problem with simplify_indexed(). This function tries, among other
things, if some terms cancel after symmetrizing with respect to dummy
indices. It then checks if the symmetrization really lessened the number
of terms and if this is the case, returns the result of the
symmetrization. However, this may miss opportunities for simplification.
Imagine that we have a sum of four terms, say,

	T = t1 + t2 + t3 + t4,

and that after symmetrization t1 and t2 cancel. However, it turns out that
symmetrization yields six terms. We have

	T_symm = t3,1 + ... + t3,6 + t4,1 + ... + t4,6.

This contains twelve terms, which is more than the original four ones, and
therefore, is not accepted as a simplification. However, as a (more or
less) intelligent human being I do not find it very difficult to use the
just presented information to think of a simpler expression for T.

Furthermore, even if, in such a case, symmetrization yields a smaller
number of terms, the result may be somewhat awkward. What if I have
objects carrying a three-dimensional index for which I use symbols a, b,
c, etcetera and a four-dimensional index for which I use symbols mu, nu,
sig, etcetera. The symmetrized result may have latin letters in places
where I would normally expect greek ones and vice versa.

Both problems sketched above are solved if one uses the symmetrized terms
for determining what cancels but the original terms to get the, hopefully,
simplified result. A patch is attached that is supposed to make this
happen.

Greetings,
Chris Dams
-------------- next part --------------
*** indexed.cpp	Fri Nov  8 18:08:45 2002
--- indexed.cpp.modified	Fri Nov  8 18:07:43 2002
***************
*** 829,834 ****
--- 829,869 ----
  		return r;
  }
  
+ struct symminfo {
+ 	ex symmterm;
+ 	ex coeff;
+ 	ex orig;
+ 	symminfo (const ex & symmterm_, const ex & orig_) {
+ 		if ( is_a<mul>(orig_) )
+ 		{	ex tmp = orig_.op(orig_.nops()-1);
+ 			orig = orig_ / tmp;
+ 		} else 
+ 			orig = orig_;
+ 		if ( is_a<mul>(symmterm_) ) {
+ 			coeff = symmterm_.op(symmterm_.nops()-1);
+ 			symmterm = symmterm_ / coeff;
+ 		} else {
+ 			coeff = 1;
+ 			symmterm = symmterm_;
+ 		}
+ 	}
+ };
+ 
+ class symminfo_is_less {
+ public:
+ 	bool operator() ( const symminfo & si1, const symminfo & si2 ) {
+ 		int comp = si1.symmterm.compare(si2.symmterm);
+ 		if ( comp < 0 ) return true;
+ 		if ( comp > 0 ) return false;
+ 		comp = si1.orig.compare(si2.orig);
+ 		if ( comp < 0 ) return true;
+ 		if ( comp > 0 ) return false;
+ 		comp = si1.coeff.compare(si2.coeff);
+ 		if ( comp < 0 ) 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)
  {
***************
*** 897,910 ****
  			lst dummy_syms;
  			for (int i=0; i<dummy_indices.size(); i++)
  				dummy_syms.append(dummy_indices[i].op(0));
! 			ex sum_symm = sum.symmetrize(dummy_syms);
! 			if (sum_symm.is_zero()) {
! 				free_indices.clear();
! 				return _ex0;
  			}
! 			int num_terms = (is_a<add>(sum_symm) ? sum_symm.nops() : 1);
! 			if (num_terms < num_terms_orig)
! 				return sum_symm;
  		}
  
  		return sum;
--- 932,965 ----
  			lst dummy_syms;
  			for (int i=0; i<dummy_indices.size(); i++)
  				dummy_syms.append(dummy_indices[i].op(0));
! 			std::vector<symminfo>v;
! 			for (int i=0; i<sum.nops(); i++) {
! 				ex sum_symm = sum.op(i).symmetrize(dummy_syms);
! 				if(is_a<add>(sum_symm))
! 					for (int j=0; j<sum_symm.nops(); j++)
! 						v.push_back( 
! 							symminfo(sum_symm.op(j),
! 								sum.op(i)) );
! 				else
! 					v.push_back( symminfo(sum_symm,
! 								sum.op(i)) );
  			}
! 			exvector result;
! 			std::sort(v.begin(), v.end(), symminfo_is_less());
! 			for (std::vector<symminfo>::iterator i=v.begin();
! 								i!=v.end();) {
! 				std::vector<symminfo>::iterator j=i;
! 				for (j++; j!=v.end()
! 					&& i->symmterm == j->symmterm; j++);
! 				for (std::vector<symminfo>::iterator k=i; k!=j;
! 									k++)
! 					result.push_back((k->coeff)*(i->orig));
! 				i=j;
! 			}
! 			ex sum_symm=add(result);
! 			if (sum_symm.is_zero())
! 				free_indices.clear();
! 			return sum_symm;
  		}
  
  		return sum;