- Cleanups: My evil plot of making ex::bp private may finally be carried

[ginac.git] / ginac / indexed.cpp

diff --git a/ginac/indexed.cpp b/ginac/indexed.cpp
index 56fd0b6b4080ebaafa72845978bfea39d5198ae9..5759b24d4c93e6f21e64a07e9ed392145e62db70 100644 (file)
--- a/ginac/indexed.cpp
+++ b/ginac/indexed.cpp
@@ -174,7 +174,7 @@ indexed::indexed(const archive_node &n, const lst &sym_lst) : inherited(n, sym_l
                                symtree = sy_none();
                                break;
                }
-               ex_to_nonconst_symmetry(symtree).validate(seq.size() - 1);
+               const_cast<symmetry &>(ex_to<symmetry>(symtree)).validate(seq.size() - 1);
        }
 }
 
@@ -250,7 +250,7 @@ bool indexed::all_index_values_are(unsigned inf) const
 
 int indexed::compare_same_type(const basic & other) const
 {
-       GINAC_ASSERT(is_of_type(other, indexed));
+       GINAC_ASSERT(is_a<indexed>(other));
        return inherited::compare_same_type(other);
 }
 
@@ -277,7 +277,7 @@ ex indexed::eval(int level) const
        // Canonicalize indices according to the symmetry properties
        if (seq.size() > 2) {
                exvector v = seq;
-               GINAC_ASSERT(is_ex_exactly_of_type(symtree, symmetry));
+               GINAC_ASSERT(is_exactly_a<symmetry>(symtree));
                int sig = canonicalize(v.begin() + 1, ex_to<symmetry>(symtree));
                if (sig != INT_MAX) {
                        // Something has changed while sorting indices, more evaluations later
@@ -288,22 +288,22 @@ ex indexed::eval(int level) const
        }
 
        // Let the class of the base object perform additional evaluations
-       return base.bp->eval_indexed(*this);
+       return ex_to<basic>(base).eval_indexed(*this);
 }
 
 int indexed::degree(const ex & s) const
 {
-       return is_equal(*s.bp) ? 1 : 0;
+       return is_equal(ex_to<basic>(s)) ? 1 : 0;
 }
 
 int indexed::ldegree(const ex & s) const
 {
-       return is_equal(*s.bp) ? 1 : 0;
+       return is_equal(ex_to<basic>(s)) ? 1 : 0;
 }
 
 ex indexed::coeff(const ex & s, int n) const
 {
-       if (is_equal(*s.bp))
+       if (is_equal(ex_to<basic>(s)))
                return n==1 ? _ex1() : _ex0();
        else
                return n==0 ? ex(*this) : _ex0();
@@ -406,7 +406,7 @@ void indexed::validate(void) const
        if (!symtree.is_zero()) {
                if (!is_ex_exactly_of_type(symtree, symmetry))
                        throw(std::invalid_argument("symmetry of indexed object must be of type symmetry"));
-               ex_to_nonconst_symmetry(symtree).validate(seq.size() - 1);
+               const_cast<symmetry &>(ex_to<symmetry>(symtree)).validate(seq.size() - 1);
        }
 }
 
@@ -668,36 +668,13 @@ try_again:
                                }
                        }
 
-                       // Contraction of symmetric with antisymmetric object is zero
-                       if (num_dummies > 1
-                        && ex_to<symmetry>(ex_to<indexed>(*it1).symtree).has_symmetry()
-                        && ex_to<symmetry>(ex_to<indexed>(*it2).symtree).has_symmetry()) {
-
-                               // Check all pairs of dummy indices
-                               for (unsigned idx1=0; idx1<num_dummies-1; idx1++) {
-                                       for (unsigned idx2=idx1+1; idx2<num_dummies; idx2++) {
-
-                                               // Try and swap the index pair and check whether the
-                                               // relative sign changed
-                                               lst subs_lst(dummy[idx1].op(0), dummy[idx2].op(0)), repl_lst(dummy[idx2].op(0), dummy[idx1].op(0));
-                                               ex swapped1 = it1->subs(subs_lst, repl_lst);
-                                               ex swapped2 = it2->subs(subs_lst, repl_lst);
-                                               if (it1->is_equal(swapped1) && it2->is_equal(-swapped2)
-                                                || it1->is_equal(-swapped1) && it2->is_equal(swapped2)) {
-                                                       free_indices.clear();
-                                                       return _ex0();
-                                               }
-                                       }
-                               }
-                       }
-
                        // Try to contract the first one with the second one
-                       contracted = it1->op(0).bp->contract_with(it1, it2, v);
+                       contracted = ex_to<basic>(it1->op(0)).contract_with(it1, it2, v);
                        if (!contracted) {
 
                                // That didn't work; maybe the second object knows how to
                                // contract itself with the first one
-                               contracted = it2->op(0).bp->contract_with(it2, it1, v);
+                               contracted = ex_to<basic>(it2->op(0)).contract_with(it2, it1, v);
                        }
                        if (contracted) {
 contraction_done:
@@ -749,13 +726,26 @@ contraction_done:
        else
                r = e;
 
+       // The result should be symmetric with respect to exchange of dummy
+       // indices, so if the symmetrization vanishes, the whole expression is
+       // zero. This detects things like eps.i.j.k * p.j * p.k = 0.
+       if (local_dummy_indices.size() >= 2) {
+               lst dummy_syms;
+               for (int i=0; i<local_dummy_indices.size(); i++)
+                       dummy_syms.append(local_dummy_indices[i].op(0));
+               if (r.symmetrize(dummy_syms).is_zero()) {
+                       free_indices.clear();
+                       return _ex0();
+               }
+       }
+
        // Dummy index renaming
        r = rename_dummy_indices(r, dummy_indices, local_dummy_indices);
 
        // Product of indexed object with a scalar?
        if (is_ex_exactly_of_type(r, mul) && r.nops() == 2
         && is_ex_exactly_of_type(r.op(1), numeric) && is_ex_of_type(r.op(0), indexed))
-               return r.op(0).op(0).bp->scalar_mul_indexed(r.op(0), ex_to<numeric>(r.op(1)));
+               return ex_to<basic>(r.op(0).op(0)).scalar_mul_indexed(r.op(0), ex_to<numeric>(r.op(1)));
        else
                return r;
 }
@@ -794,7 +784,7 @@ ex simplify_indexed(const ex & e, exvector & free_indices, exvector & dummy_indi
                                        if (!indices_consistent(free_indices, free_indices_of_term))
                                                throw (std::runtime_error("simplify_indexed: inconsistent indices in sum"));
                                        if (is_ex_of_type(sum, indexed) && is_ex_of_type(term, indexed))
-                                               sum = sum.op(0).bp->add_indexed(sum, term);
+                                               sum = ex_to<basic>(sum.op(0)).add_indexed(sum, term);
                                        else
                                                sum += term;
                                }