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;
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;
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);
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());
}
++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)
/** 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.
// 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 {
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());
/* Rename indices in the static members of the product */
exvector expanded_seq_mod;
size_t j = 0;
- size_t i;
- for (i=0; i<expanded_seq.size(); i++) {
+ exvector va;
+
+ for (size_t 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]));
+ for (size_t i=0; i<number_of_adds; 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))->
// 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;
}