Index: ginac/expairseq.cpp =================================================================== RCS file: /home/cvs/GiNaC/ginac/expairseq.cpp,v retrieving revision 1.77 diff -r1.77 expairseq.cpp 36a37 > #include "indexed.h" 760,761c761,769 < construct_from_2_expairseq(ex_to(lh), < ex_to(rh)); --- > if(is_a(lh)) > { > ex newrh=rename_dummy_indices_uniquely(lh, rh); > construct_from_2_expairseq(ex_to(lh), > ex_to(newrh)); > } > else > construct_from_2_expairseq(ex_to(lh), > ex_to(rh)); 1010a1019 > exvector dummy_indices; 1014,1017c1023,1038 < const expairseq &subseqref = ex_to(*cit); < combine_overall_coeff(subseqref.overall_coeff); < epvector::const_iterator cit_s = subseqref.seq.begin(); < while (cit_s!=subseqref.seq.end()) { --- > const expairseq *subseqref; > ex newfactor; > if(is_a(*cit)) > { > exvector dummies_of_factor = get_all_dummy_indices(*cit); > newfactor = rename_dummy_indices_uniquely(dummy_indices, dummies_of_factor, *cit); > subseqref = &(ex_to(newfactor)); > exvector new_dummy_indices; > set_union(dummy_indices.begin(), dummy_indices.end(), dummies_of_factor.begin(), dummies_of_factor.end(), std::back_insert_iterator(new_dummy_indices), ex_is_less()); > dummy_indices.swap(new_dummy_indices); > } > else > subseqref = &ex_to(*cit); > combine_overall_coeff(subseqref->overall_coeff); > epvector::const_iterator cit_s = subseqref->seq.begin(); > while (cit_s!=subseqref->seq.end()) { 1581a1603,1651 > class safe_inserter > { > public: > safe_inserter(const ex&); > std::auto_ptr getseq(){return epv;} > void insert_old_pair(const expair &p) > { > epv->push_back(p); > } > void insert_new_pair(const expair &p); > private: > std::auto_ptr epv; > bool dodummies; > exvector dummy_indices; > void update_dummy_indices(const exvector&); > }; > > safe_inserter::safe_inserter(const ex&e):epv(new epvector) > { > epv->reserve(e.nops()); > dodummies=is_a(e); > if(dodummies) > dummy_indices = get_all_dummy_indices(e); > } > > void safe_inserter::update_dummy_indices(const exvector &v) > { > exvector new_dummy_indices; > set_union(dummy_indices.begin(), dummy_indices.end(), v.begin(), v.end(), > std::back_insert_iterator(new_dummy_indices), ex_is_less()); > dummy_indices.swap(new_dummy_indices); > } > > void safe_inserter::insert_new_pair(const expair &p) > { > if(!dodummies) > { > epv->push_back(p); > return; > } > exvector dummies_of_factor = get_all_dummy_indices(p.rest); > if(dummies_of_factor.size() == 0) > { epv->push_back(p); > return; > } > ex newfactor = rename_dummy_indices_uniquely(dummy_indices, dummies_of_factor, p.rest); > update_dummy_indices(dummies_of_factor); > epv -> push_back(expair(newfactor, p.coeff)); > } 1617,1618c1687 < std::auto_ptr s(new epvector); < s->reserve(seq.size()); --- > safe_inserter s(*this); 1621c1690,1691 < s->insert(s->begin(), seq.begin(), cit); --- > for(epvector::const_iterator i=seq.begin(); i!=cit; ++i) > s.insert_old_pair(*i); 1624c1694 < s->push_back(split_ex_to_pair(subsed_ex)); --- > s.insert_new_pair(split_ex_to_pair(subsed_ex)); 1629c1699,1704 < s->push_back(split_ex_to_pair(recombine_pair_to_ex(*cit).subs(m, options))); --- > ex orig_ex = recombine_pair_to_ex(*cit); > ex subsed_ex = orig_ex.subs(m, options); > if(are_ex_trivially_equal(orig_ex, subsed_ex)) > s.insert_old_pair(*cit); > else > s.insert_new_pair(split_ex_to_pair(subsed_ex)); 1632c1707 < return s; --- > return s.getseq(); 1648,1649c1723 < std::auto_ptr s(new epvector); < s->reserve(seq.size()); --- > safe_inserter s(*this); 1652c1726,1727 < s->insert(s->begin(), seq.begin(), cit); --- > for(epvector::const_iterator i=seq.begin(); i!=cit; ++i) > s.insert_old_pair(*i); 1655c1730 < s->push_back(combine_ex_with_coeff_to_pair(subsed_ex, cit->coeff)); --- > s.insert_new_pair(combine_ex_with_coeff_to_pair(subsed_ex, cit->coeff)); 1660,1661c1735,1740 < s->push_back(combine_ex_with_coeff_to_pair(cit->rest.subs(m, options), < cit->coeff)); --- > ex orig_ex = cit->rest; > ex subsed_ex = orig_ex.subs(m, options); > if(are_ex_trivially_equal(orig_ex, subsed_ex)) > s.insert_old_pair(*cit); > else > s.insert_new_pair(combine_ex_with_coeff_to_pair(subsed_ex, cit->coeff)); 1664c1743 < return s; --- > return s.getseq(); Index: ginac/indexed.cpp =================================================================== RCS file: /home/cvs/GiNaC/ginac/indexed.cpp,v retrieving revision 1.96 diff -r1.96 indexed.cpp 534a535,543 > template 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(*i)) > ++number; > return number; > } > 543c552 < static ex rename_dummy_indices(const ex & e, exvector & global_dummy_indices, exvector & local_dummy_indices) --- > template static ex rename_dummy_indices(const ex & e, exvector & global_dummy_indices, exvector & local_dummy_indices) 545,546c554,555 < size_t global_size = global_dummy_indices.size(), < local_size = local_dummy_indices.size(); --- > size_t global_size = number_of_type(global_dummy_indices), > local_size = number_of_type(local_dummy_indices); 560c569 < if (find_if(global_dummy_indices.begin(), global_dummy_indices.end(), bind2nd(op0_is_equal(), *it)) == global_dummy_indices.end()) { --- > if (is_exactly_a(*it) && find_if(global_dummy_indices.begin(), global_dummy_indices.end(), bind2nd(idx_is_equal_ignore_dim(), *it)) == global_dummy_indices.end()) { 578,579c587,589 < for (size_t i=0; i for (size_t i=0; local_syms.size()!=local_size; i++) > if(is_exactly_a(local_dummy_indices[i])) > local_syms.push_back(local_dummy_indices[i].op(0)); 581,582c591,593 < for (size_t i=0; i for (size_t i=0; global_syms.size()!=local_size; i++) // don't use more global symbols than necessary > if(is_exactly_a(global_dummy_indices[i])) > global_syms.push_back(global_dummy_indices[i].op(0)); 706a718,729 > template 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(*it)) > dummy_syms.push_back(it->op(0)); > if(dummy_syms.size() < 2) > return r; > ex q=symmetrize(r, dummy_syms); > return q; > } > 867,875c890,903 < 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; < } --- > ex q = idx_symmetrization(r, local_dummy_indices); > if (q.is_zero()) { > free_indices.clear(); > return _ex0; > } > q = idx_symmetrization(q, local_dummy_indices); > if (q.is_zero()) { > free_indices.clear(); > return _ex0; > } > q = idx_symmetrization(q, local_dummy_indices); > if (q.is_zero()) { > free_indices.clear(); > return _ex0; 879c907,909 < r = rename_dummy_indices(r, dummy_indices, local_dummy_indices); --- > r = rename_dummy_indices(r, dummy_indices, local_dummy_indices); > r = rename_dummy_indices(r, dummy_indices, local_dummy_indices); > r = rename_dummy_indices(r, dummy_indices, local_dummy_indices); 945a976,986 > bool hasindex(const ex &x, const ex &sym) > { > if(is_a(x) && x.op(0)==sym) > return true; > else > for(size_t i=0; i if(hasindex(x.op(i), sym)) > return true; > return false; > } > 974c1015,1018 < return rename_dummy_indices(e_expanded, dummy_indices, local_dummy_indices); --- > e_expanded = rename_dummy_indices(e_expanded, dummy_indices, local_dummy_indices); > e_expanded = rename_dummy_indices(e_expanded, dummy_indices, local_dummy_indices); > e_expanded = rename_dummy_indices(e_expanded, dummy_indices, local_dummy_indices); > return e_expanded; 1018,1023d1061 < // 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)); < 1029c1067,1074 < ex term_symm = symmetrize(term, dummy_syms); --- > 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(term, dummy_indices_of_term); > term_symm = idx_symmetrization(term_symm, dummy_indices_of_term); > term_symm = idx_symmetrization(term_symm, dummy_indices_of_term); 1321c1366 < ex rename_dummy_indices_uniquely(const ex & a, const ex & b) --- > ex rename_dummy_indices_uniquely(const exvector & va, const exvector & vb, const ex & b) 1323c1368 < exvector va = get_all_dummy_indices(a), vb = get_all_dummy_indices(b), common_indices; --- > exvector common_indices; 1333,1336c1378,1394 < if (is_a(*ip)) { < varidx mu((new symbol)->setflag(status_flags::dynallocated), ex_to(*ip).get_dim(), ex_to(*ip).is_covariant()); < old_indices.push_back(*ip); < new_indices.push_back(mu); --- > ex newsym=(new symbol)->setflag(status_flags::dynallocated); > ex newidx; > if(is_exactly_a(*ip)) > newidx = (new spinidx(newsym, ex_to(*ip).get_dim(), > ex_to(*ip).is_covariant(), > ex_to(*ip).is_dotted())) > -> setflag(status_flags::dynallocated); > else if (is_exactly_a(*ip)) > newidx = (new varidx(newsym, ex_to(*ip).get_dim(), > ex_to(*ip).is_covariant())) > -> setflag(status_flags::dynallocated); > else > newidx = (new idx(newsym, ex_to(*ip).get_dim())) > -> setflag(status_flags::dynallocated); > old_indices.push_back(*ip); > new_indices.push_back(newidx); > if(is_a(*ip)) { 1338,1341c1396 < new_indices.push_back(mu.toggle_variance()); < } else { < old_indices.push_back(*ip); < new_indices.push_back(idx((new symbol)->setflag(status_flags::dynallocated), ex_to(*ip).get_dim())); --- > new_indices.push_back(ex_to(newidx).toggle_variance()); 1346a1402,1408 > } > > ex rename_dummy_indices_uniquely(const ex & a, const ex & b) > { > exvector va = get_all_dummy_indices(a); > exvector vb = get_all_dummy_indices(b); > return rename_dummy_indices_uniquely(va, vb, b); Index: ginac/indexed.h =================================================================== RCS file: /home/cvs/GiNaC/ginac/indexed.h,v retrieving revision 1.52 diff -r1.52 indexed.h 258a259,261 > /** Same as above, where va and vb contain the indices of a and b and are sorted */ > ex rename_dummy_indices_uniquely(const exvector & va, const exvector & vb, const ex & b); >