* Restore speed lost in GiNaC-1.3.2 [V. Kisil].

diff --git a/check/exam_clifford.cpp b/check/exam_clifford.cpp
index 15daef1750bea9893149fe5ec56337dbe8cd56e4..c0650a478bc43f7bc119f8655ca88bf82a27e8da 100644 (file)
--- 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)
 {
 
 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);
                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;
                        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);
        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 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);
        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);
 
        // 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 64eba6d08cf1cccf6dbf455cbcce4c01356dd3e7..e9c4fee2e56e9bd09b9c36a8bf73cb4f79d2a820 100644 (file)
--- a/ginac/indexed.cpp
+++ b/ginac/indexed.cpp
@@ -1374,12 +1374,12 @@ exvector get_all_dummy_indices(const ex & e)
        return v;
 }
 
        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()) {
 {
        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());
        } 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;
                }
                        }
                        ++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);
 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)
 }
 
 ex expand_dummy_sum(const ex & e, bool subs_idx)

diff --git a/ginac/indexed.h b/ginac/indexed.h
index 1ce81945b432cdc40f0537a7699eb35441043eb8..8f2779ffbf34b59be5a4c73a249fcf3a10273175 100644 (file)
--- a/ginac/indexed.h
+++ b/ginac/indexed.h
@@ -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 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);
 
 /** 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.
 /** 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 90fd772ec2b16552c390e4697bfc30d3ca2871a8..f5c202252ba4f8f9a351107075680673ba0f6995 100644 (file)
--- 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);
 
                                // 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:
                                // 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();
                                        // 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:
                                                // 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 {
                                                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);
                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;
 
                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());
                        ex term = (new mul(factors, overall_coeff))->setflag(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 eeadf4e7feb4fa6acf7b48db84062b049e502526..c72b608ae3756f19119849098cba08f82e01564b 100644 (file)
--- 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;
        /* 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 {
                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;
                }
        }
 
        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))->
                }
 
                distrseq.push_back((new ncmul(term, true))->

diff --git a/ginac/power.cpp b/ginac/power.cpp
index 3dd0f48f487cc6c397ae13e70019f41e9fe2f7e9..f9dfafbfb80c55a5856195a4be9127a0f538c0fc 100644 (file)
--- 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;
        // 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++)
                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;
        }
 
                return result;
        }