* Restore speed lost in GiNaC-1.3.2 [V. Kisil].
authorRichard Kreckel <Richard.Kreckel@uni-mainz.de>
Fri, 18 Nov 2005 01:26:24 +0000 (01:26 +0000)
committerRichard Kreckel <Richard.Kreckel@uni-mainz.de>
Fri, 18 Nov 2005 01:26:24 +0000 (01:26 +0000)
check/exam_clifford.cpp
ginac/indexed.cpp
ginac/indexed.h
ginac/mul.cpp
ginac/ncmul.cpp
ginac/power.cpp

index 15daef1..c0650a4 100644 (file)
@@ -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>(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);
index 64eba6d..e9c4fee 100644 (file)
@@ -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<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 @@ 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<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 1ce8194..8f2779f 100644 (file)
@@ -254,12 +254,20 @@ template<> inline bool is_exactly_a<indexed>(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.
index 90fd772..f5c2022 100644 (file)
@@ -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<numeric>(rest)) {
                                                        oc += ex_to<numeric>(rest).mul(ex_to<numeric>(i1->coeff).mul(ex_to<numeric>(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; 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 eeadf4e..c72b608 100644 (file)
@@ -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<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))->
index 3dd0f48..f9dfafb 100644 (file)
@@ -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;
        }