[GiNaC-devel] Patch for products of dummy summations
Vladimir Kisil
kisilv at maths.leeds.ac.uk
Wed Nov 16 22:55:19 CET 2005
> Was that patch for the 1.3 branch or for HEAD? In any case: When I try
> applying it to any of the two, some hunks fail.
I tried my patch on HEAD, but it may work with 1.3 branch as
well. I did syncronisation of CVS today and am attaching "diff -u" version
of this patch to this message. Hope it will work now.
Best wishes,
Vladimir
--
Vladimir Kisil email: kisilv at amsta.leeds.ac.uk
www: http://amsta.leeds.ac.uk/~kisilv
-------------- next part --------------
Index: ginac/mul.cpp
===================================================================
RCS file: /home/cvs/GiNaC/ginac/mul.cpp,v
retrieving revision 1.90
diff -u -r1.90 mul.cpp
--- ginac/mul.cpp 24 Jul 2005 21:02:39 -0000 1.90
+++ ginac/mul.cpp 16 Nov 2005 21:50:30 -0000
@@ -895,16 +895,33 @@
// Compute the new overall coefficient and put it together:
ex tmp_accu = (new add(distrseq, add1.overall_coeff*add2.overall_coeff))->setflag(status_flags::dynallocated);
+ exvector add1_dummy_indices, add2_dummy_indices, add_indices;
+
+ for (epvector::const_iterator i=add1begin; i!=add1end; ++i) {
+ add_indices = get_all_dummy_indices(i->rest);
+ add1_dummy_indices.insert(add1_dummy_indices.end(), add_indices.begin(), add_indices.end());
+ }
+ for (epvector::const_iterator i=add2begin; i!=add2end; ++i) {
+ add_indices = get_all_dummy_indices(i->rest);
+ add2_dummy_indices.insert(add2_dummy_indices.end(), add_indices.begin(), add_indices.end());
+ }
+
+ sort(add1_dummy_indices.begin(), add1_dummy_indices.end(), ex_is_less());
+ sort(add2_dummy_indices.begin(), add2_dummy_indices.end(), ex_is_less());
+ lst dummy_subs = rename_dummy_indices_uniquely(add1_dummy_indices, add2_dummy_indices);
+
// Multiply explicitly all non-numeric terms of add1 and add2:
- for (epvector::const_iterator i1=add1begin; i1!=add1end; ++i1) {
+ for (epvector::const_iterator i2=add2begin; i2!=add2end; ++i2) {
// We really have to combine terms here in order to compactify
// the result. Otherwise it would become waayy tooo bigg.
numeric oc;
distrseq.clear();
- for (epvector::const_iterator i2=add2begin; i2!=add2end; ++i2) {
+ ex i2_new = (dummy_subs.op(0).nops()>0?
+ i2->rest.subs((lst)dummy_subs.op(0), (lst)dummy_subs.op(1), subs_options::no_pattern) : i2->rest);
+ for (epvector::const_iterator i1=add1begin; i1!=add1end; ++i1) {
// Don't push_back expairs which might have a rest that evaluates to a numeric,
// since that would violate an invariant of expairseq:
- const ex rest = (new mul(i1->rest, rename_dummy_indices_uniquely(i1->rest, i2->rest)))->setflag(status_flags::dynallocated);
+ const ex rest = (new mul(i1->rest, i2_new))->setflag(status_flags::dynallocated);
if (is_exactly_a<numeric>(rest)) {
oc += ex_to<numeric>(rest).mul(ex_to<numeric>(i1->coeff).mul(ex_to<numeric>(i2->coeff)));
} else {
@@ -932,10 +949,12 @@
size_t n = last_expanded.nops();
exvector distrseq;
distrseq.reserve(n);
+ exvector va = get_all_dummy_indices(mul(non_adds));
+ sort(va.begin(), va.end(), ex_is_less());
for (size_t i=0; i<n; ++i) {
epvector factors = non_adds;
- factors.push_back(split_ex_to_pair(rename_dummy_indices_uniquely(mul(non_adds), last_expanded.op(i))));
+ factors.push_back(split_ex_to_pair(rename_dummy_indices_uniquely(va, last_expanded.op(i))));
ex term = (new mul(factors, overall_coeff))->setflag(status_flags::dynallocated);
if (can_be_further_expanded(term)) {
distrseq.push_back(term.expand());
Index: ginac/ncmul.cpp
===================================================================
RCS file: /home/cvs/GiNaC/ginac/ncmul.cpp,v
retrieving revision 1.57
diff -u -r1.57 ncmul.cpp
--- ginac/ncmul.cpp 19 May 2005 14:10:40 -0000 1.57
+++ ginac/ncmul.cpp 16 Nov 2005 21:50:30 -0000
@@ -170,21 +170,22 @@
/* Rename indices in the static members of the product */
exvector expanded_seq_mod;
- size_t j = 0;
- size_t i;
+ size_t j = 0, i;
+ exvector va;
+
for (i=0; i<expanded_seq.size(); i++) {
if (i == positions_of_adds[j]) {
expanded_seq_mod.push_back(_ex1);
j++;
} else {
- expanded_seq_mod.push_back(rename_dummy_indices_uniquely(ncmul(expanded_seq_mod), expanded_seq[i]));
+ expanded_seq_mod.push_back(rename_dummy_indices_uniquely(va, expanded_seq[i], true));
}
}
while (true) {
exvector term = expanded_seq_mod;
for (i=0; i<number_of_adds; i++) {
- term[positions_of_adds[i]] = rename_dummy_indices_uniquely(ncmul(term), expanded_seq[positions_of_adds[i]].op(k[i]));
+ term[positions_of_adds[i]] = rename_dummy_indices_uniquely(va, expanded_seq[positions_of_adds[i]].op(k[i]), true);
}
distrseq.push_back((new ncmul(term, true))->
Index: ginac/power.cpp
===================================================================
RCS file: /home/cvs/GiNaC/ginac/power.cpp,v
retrieving revision 1.102
diff -u -r1.102 power.cpp
--- ginac/power.cpp 9 Nov 2005 13:47:00 -0000 1.102
+++ ginac/power.cpp 16 Nov 2005 21:50:30 -0000
@@ -864,8 +864,11 @@
// Leave it to multiplication since dummy indices have to be renamed
if (get_all_dummy_indices(m).size() > 0 && n.is_positive()) {
ex result = m;
+ exvector va = get_all_dummy_indices(m);
+ sort(va.begin(), va.end(), ex_is_less());
+
for (int i=1; i < n.to_int(); i++)
- result *= rename_dummy_indices_uniquely(m,m);
+ result *= rename_dummy_indices_uniquely(va, m);
return result;
}
Index: ginac/indexed.h
===================================================================
RCS file: /home/cvs/GiNaC/ginac/indexed.h,v
retrieving revision 1.53
diff -u -r1.53 indexed.h
--- ginac/indexed.h 11 Nov 2005 12:05:24 -0000 1.53
+++ ginac/indexed.h 16 Nov 2005 21:50:30 -0000
@@ -254,12 +254,20 @@
/** Returns all dummy indices from the expression */
exvector get_all_dummy_indices(const ex & e);
+/** Returns b with all dummy indices, which are listed in va, renamed
+ * if modify_va is set to TRUE all dummy indices of b will be appended to va */
+ex rename_dummy_indices_uniquely(exvector & va, const ex & b, bool modify_va = false);
+
/** Returns b with all dummy indices, which are common with a, renamed */
ex rename_dummy_indices_uniquely(const ex & a, const ex & b);
/** 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);
+/** Similar to above, where va and vb are the same and the return value is a list of two lists
+ * for substitution in b */
+lst rename_dummy_indices_uniquely(const exvector & va, const exvector & vb);
+
/** This function returns the given expression with expanded sums
* for all dummy index summations, where the dimensionality of
* the dummy index is numeric.
Index: ginac/indexed.cpp
===================================================================
RCS file: /home/cvs/GiNaC/ginac/indexed.cpp,v
retrieving revision 1.97
diff -u -r1.97 indexed.cpp
--- ginac/indexed.cpp 11 Nov 2005 12:05:24 -0000 1.97
+++ ginac/indexed.cpp 16 Nov 2005 21:50:30 -0000
@@ -1374,12 +1374,12 @@
return v;
}
-ex rename_dummy_indices_uniquely(const exvector & va, const exvector & vb, const ex & b)
+lst rename_dummy_indices_uniquely(const exvector & va, const exvector & vb)
{
exvector common_indices;
set_intersection(va.begin(), va.end(), vb.begin(), vb.end(), std::back_insert_iterator<exvector>(common_indices), ex_is_less());
if (common_indices.empty()) {
- return b;
+ return lst(lst(), lst());
} else {
exvector new_indices, old_indices;
old_indices.reserve(2*common_indices.size());
@@ -1408,17 +1408,57 @@
}
++ip;
}
- return b.subs(lst(old_indices.begin(), old_indices.end()), lst(new_indices.begin(), new_indices.end()), subs_options::no_pattern);
+ return lst(lst(old_indices.begin(), old_indices.end()), lst(new_indices.begin(), new_indices.end()));
}
}
+ex rename_dummy_indices_uniquely(const exvector & va, const exvector & vb, const ex & b)
+{
+ lst indices_subs = rename_dummy_indices_uniquely(va, vb);
+ return (indices_subs.op(0).nops()>0 ? b.subs((lst)indices_subs.op(0), (lst)indices_subs.op(1), subs_options::no_pattern) : b);
+}
+
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);
- sort(va.begin(), va.end(), ex_is_less());
- sort(vb.begin(), vb.end(), ex_is_less());
- return rename_dummy_indices_uniquely(va, vb, b);
+ if (va.size() > 0) {
+ exvector vb = get_all_dummy_indices(b);
+ if (vb.size() > 0) {
+ sort(va.begin(), va.end(), ex_is_less());
+ sort(vb.begin(), vb.end(), ex_is_less());
+ lst indices_subs = rename_dummy_indices_uniquely(va, vb);
+ if (indices_subs.op(0).nops() > 0)
+ return b.subs((lst)indices_subs.op(0), (lst)indices_subs.op(1), subs_options::no_pattern);
+ }
+ }
+ return b;
+}
+
+ex rename_dummy_indices_uniquely(exvector & va, const ex & b, bool modify_va)
+{
+ if (va.size() > 0) {
+ exvector vb = get_all_dummy_indices(b);
+ if (vb.size() > 0) {
+ sort(vb.begin(), vb.end(), ex_is_less());
+ lst indices_subs = rename_dummy_indices_uniquely(va, vb);
+ if (indices_subs.op(0).nops() > 0) {
+ if (modify_va) {
+ for (lst::const_iterator i = ex_to<lst>(indices_subs.op(1)).begin(); i != ex_to<lst>(indices_subs.op(1)).end(); ++i)
+ va.push_back(*i);
+ exvector uncommon_indices;
+ set_difference(vb.begin(), vb.end(), indices_subs.op(0).begin(), indices_subs.op(0).end(), std::back_insert_iterator<exvector>(uncommon_indices), ex_is_less());
+ exvector::const_iterator ip = uncommon_indices.begin(), ipend = uncommon_indices.end();
+ while (ip != ipend) {
+ va.push_back(*ip);
+ ++ip;
+ }
+ sort(va.begin(), va.end(), ex_is_less());
+ }
+ return b.subs((lst)indices_subs.op(0), (lst)indices_subs.op(1), subs_options::no_pattern);
+ }
+ }
+ }
+ return b;
}
ex expand_dummy_sum(const ex & e, bool subs_idx)
Index: check/exam_clifford.cpp
===================================================================
RCS file: /home/cvs/GiNaC/check/exam_clifford.cpp,v
retrieving revision 1.27
diff -u -r1.27 exam_clifford.cpp
--- check/exam_clifford.cpp 12 Jul 2005 17:56:29 -0000 1.27
+++ check/exam_clifford.cpp 16 Nov 2005 21:50:30 -0000
@@ -48,9 +48,9 @@
static unsigned check_equal_lst(const ex & e1, const ex & e2)
{
- for(int i = 0; i++; i < e1.nops()) {
+ for (unsigned int i = 0; i < e1.nops(); i++) {
ex e = e1.op(i) - e2.op(i);
- if (!e.is_zero()) {
+ if (!e.normal().is_zero()) {
clog << "(" << e1 << ") - (" << e2 << ") erroneously returned "
<< e << " instead of 0 (in the entry " << i << ")" << endl;
return 1;
@@ -314,7 +314,7 @@
matrix A_symm(4,4), A2(4, 4);
A_symm = A.add(A.transpose()).mul(half);
A2 = A_symm.mul(A_symm);
-
+
ex e, e1;
bool anticommuting = ex_to<clifford>(clifford_unit(nu, A)).is_anticommuting();
int result = 0;
@@ -412,7 +412,7 @@
ex c = clifford_unit(nu, A, 1);
e = lst_to_clifford(lst(t, x, y, z), mu, A, 1) * lst_to_clifford(lst(1, 2, 3, 4), c);
e1 = clifford_inverse(e);
- result += check_equal_lst((e*e1).simplify_indexed(), dirac_ONE(1));
+ result += check_equal((e*e1).simplify_indexed(), dirac_ONE(1));
// Moebius map (both forms) checks for symmetric metrics only
matrix M1(2, 2), M2(2, 2);
More information about the GiNaC-devel
mailing list