From: Richard Kreckel Date: Fri, 18 Nov 2005 01:26:24 +0000 (+0000) Subject: * Restore speed lost in GiNaC-1.3.2 [V. Kisil]. X-Git-Tag: release_1-4-0~130 X-Git-Url: https://www.ginac.de/ginac.git//ginac.git?p=ginac.git;a=commitdiff_plain;h=2bb16a36176246de7cd5c1ab2044c43f5a35115b * Restore speed lost in GiNaC-1.3.2 [V. Kisil]. --- diff --git a/check/exam_clifford.cpp b/check/exam_clifford.cpp index 15daef17..c0650a47 100644 --- a/check/exam_clifford.cpp +++ b/check/exam_clifford.cpp @@ -48,9 +48,9 @@ static unsigned check_equal_simplify(const ex &e1, const ex &e2) 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 @@ static unsigned clifford_check6(const matrix & A) 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_unit(nu, A)).is_anticommuting(); int result = 0; @@ -412,7 +412,7 @@ static unsigned clifford_check6(const matrix & A) 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); diff --git a/ginac/indexed.cpp b/ginac/indexed.cpp index 64eba6d0..e9c4fee2 100644 --- a/ginac/indexed.cpp +++ b/ginac/indexed.cpp @@ -1374,12 +1374,12 @@ exvector get_all_dummy_indices(const ex & e) 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(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 @@ ex rename_dummy_indices_uniquely(const exvector & va, const exvector & vb, const } ++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(indices_subs.op(1)).begin(); i != ex_to(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(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) diff --git a/ginac/indexed.h b/ginac/indexed.h index 1ce81945..8f2779ff 100644 --- a/ginac/indexed.h +++ b/ginac/indexed.h @@ -254,12 +254,20 @@ template<> inline bool is_exactly_a(const basic & obj) /** 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. diff --git a/ginac/mul.cpp b/ginac/mul.cpp index 90fd772e..f5c20225 100644 --- a/ginac/mul.cpp +++ b/ginac/mul.cpp @@ -895,16 +895,33 @@ ex mul::expand(unsigned options) const // 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(rest)) { oc += ex_to(rest).mul(ex_to(i1->coeff).mul(ex_to(i2->coeff))); } else { @@ -932,10 +949,12 @@ ex mul::expand(unsigned options) const 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; isetflag(status_flags::dynallocated); if (can_be_further_expanded(term)) { distrseq.push_back(term.expand()); diff --git a/ginac/ncmul.cpp b/ginac/ncmul.cpp index eeadf4e7..c72b608a 100644 --- a/ginac/ncmul.cpp +++ b/ginac/ncmul.cpp @@ -171,20 +171,21 @@ ex ncmul::expand(unsigned options) const /* Rename indices in the static members of the product */ exvector expanded_seq_mod; size_t j = 0; - size_t i; - for (i=0; i diff --git a/ginac/power.cpp b/ginac/power.cpp index 3dd0f48f..f9dfafbf 100644 --- a/ginac/power.cpp +++ b/ginac/power.cpp @@ -864,8 +864,11 @@ ex power::expand_mul(const mul & m, const numeric & n, unsigned options, bool fr // 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; }