[GiNaC-list] Yet another version of the patch

Chris Dams Chris.Dams at mi.infn.it
Wed Oct 5 18:51:16 CEST 2005


Hi again!

I found yet another error. This time because 
combine_ex_with_coeff_to_pair should be used to insert something into an 
expairseq. The latest version of my patch is attached.

Best wishes,
Chris
-------------- next part --------------
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<expairseq>(lh),
< 				                           ex_to<expairseq>(rh));
---
> 				if(is_a<mul>(lh))
> 				{	
> 					ex newrh=rename_dummy_indices_uniquely(lh, rh);
> 					construct_from_2_expairseq(ex_to<expairseq>(lh),
> 					                           ex_to<expairseq>(newrh));
> 				}
> 				else
> 					construct_from_2_expairseq(ex_to<expairseq>(lh),
> 					                           ex_to<expairseq>(rh));
1010a1019
> 	exvector dummy_indices;
1014,1017c1023,1038
< 			const expairseq &subseqref = ex_to<expairseq>(*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<mul>(*cit))
> 			{
> 				exvector dummies_of_factor = get_all_dummy_indices(*cit);
> 				newfactor = rename_dummy_indices_uniquely(dummy_indices, dummies_of_factor, *cit);
> 				subseqref = &(ex_to<expairseq>(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<exvector>(new_dummy_indices), ex_is_less());
> 				dummy_indices.swap(new_dummy_indices);
> 			}
> 			else
> 				subseqref = &ex_to<expairseq>(*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<epvector> 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<epvector> 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());
> 	bool dodummies=is_a<mul>(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<exvector>(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<epvector> 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<epvector> 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<class T> 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<T>(*i))
> 			++number;
> 	return number;
> }
> 
543c552
< static ex rename_dummy_indices(const ex & e, exvector & global_dummy_indices, exvector & local_dummy_indices)
---
> template<class T> 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<T>(global_dummy_indices),
> 	       local_size = number_of_type<T>(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<T>(*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<local_size; i++)
< 		local_syms.push_back(local_dummy_indices[i].op(0));
---
> 	for (size_t i=0; local_syms.size()!=local_size; i++)
> 		if(is_exactly_a<T>(local_dummy_indices[i]))
> 			local_syms.push_back(local_dummy_indices[i].op(0));
581,582c591,593
< 	for (size_t i=0; i<local_size; i++) // don't use more global symbols than necessary
< 		global_syms.push_back(global_dummy_indices[i].op(0));
---
> 	for (size_t i=0; global_syms.size()!=local_size; i++) // don't use more global symbols than necessary
> 		if(is_exactly_a<T>(global_dummy_indices[i]))
> 			global_syms.push_back(global_dummy_indices[i].op(0));
706a718,729
> 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;
> }
> 
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<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;
879c907,909
< 	r = rename_dummy_indices(r, dummy_indices, local_dummy_indices);
---
> 	r = rename_dummy_indices<idx>(r, dummy_indices, local_dummy_indices);
> 	r = rename_dummy_indices<varidx>(r, dummy_indices, local_dummy_indices);
> 	r = rename_dummy_indices<spinidx>(r, dummy_indices, local_dummy_indices);
945a976,986
> bool hasindex(const ex &x, const ex &sym)
> {	
> 	if(is_a<idx>(x) && x.op(0)==sym)
> 		return true;
> 	else
> 		for(size_t i=0; i<x.nops(); ++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<idx>(e_expanded, dummy_indices, local_dummy_indices);
> 		e_expanded = rename_dummy_indices<varidx>(e_expanded, dummy_indices, local_dummy_indices);
> 		e_expanded = rename_dummy_indices<spinidx>(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<idx>(term, dummy_indices_of_term);
> 			term_symm = idx_symmetrization<varidx>(term_symm, dummy_indices_of_term);
> 			term_symm = idx_symmetrization<spinidx>(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<varidx>(*ip)) {
< 				varidx mu((new symbol)->setflag(status_flags::dynallocated), ex_to<varidx>(*ip).get_dim(), ex_to<varidx>(*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<spinidx>(*ip))
> 				newidx = (new spinidx(newsym, ex_to<spinidx>(*ip).get_dim(),
> 						ex_to<spinidx>(*ip).is_covariant(),
> 						ex_to<spinidx>(*ip).is_dotted()))
> 					-> setflag(status_flags::dynallocated);
> 			else if (is_exactly_a<varidx>(*ip))
> 				newidx = (new varidx(newsym, ex_to<varidx>(*ip).get_dim(),
> 						ex_to<varidx>(*ip).is_covariant()))
> 					-> setflag(status_flags::dynallocated);
> 			else
> 				newidx = (new idx(newsym, ex_to<idx>(*ip).get_dim()))
> 					-> setflag(status_flags::dynallocated);
> 			old_indices.push_back(*ip);
> 			new_indices.push_back(newidx);
> 			if(is_a<varidx>(*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<varidx>(*ip).get_dim()));
---
> 				new_indices.push_back(ex_to<varidx>(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);
> 


More information about the GiNaC-list mailing list