*** 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;