From: Christian Bauer Date: Fri, 18 May 2001 00:53:07 +0000 (+0000) Subject: - simplify_indexed() renames dummy indices so, e.g., "a.i*a.i+a.j*a.j" gets X-Git-Tag: release_0-9-0~50 X-Git-Url: https://www.ginac.de/ginac.git//ginac.git?p=ginac.git;a=commitdiff_plain;h=c8f1702a09070f7d0d64118a0a3405521affb2cb - simplify_indexed() renames dummy indices so, e.g., "a.i*a.i+a.j*a.j" gets simplified to "2*a.i*a.i" (or 2*a.j*a.j, you can't know which one) - fixed possible crash when calling subs() on expressions with non-commutative products - added canonicalize_clifford() --- diff --git a/ginac/clifford.cpp b/ginac/clifford.cpp index 5c7f7c3b..a9038818 100644 --- a/ginac/clifford.cpp +++ b/ginac/clifford.cpp @@ -446,4 +446,62 @@ ex dirac_trace(const ex & e, unsigned char rl, const ex & trONE) return _ex0(); } +ex canonicalize_clifford(const ex & e) +{ + if (is_ex_exactly_of_type(e, add)) { + + ex sum = _ex0(); + for (unsigned i=0; iop(0), diracgamma5)) + itstart++; + while (next_to_last != itstart) { + exvector::iterator it = itstart; + while (it != next_to_last) { + if (it[0].op(1).compare(it[1].op(1)) > 0) { + ex save0 = it[0], save1 = it[1]; + it[0] = lorentz_g(it[0].op(1), it[1].op(1)); + it[1] = _ex2(); + ex sum = ncmul(v); + it[0] = save1; + it[1] = save0; + sum -= ncmul(v); + return canonicalize_clifford(sum); + } + it++; + } + next_to_last--; + } + return ncmul(v); + } + + return e; +} + } // namespace GiNaC diff --git a/ginac/clifford.h b/ginac/clifford.h index 42be538f..487c12d7 100644 --- a/ginac/clifford.h +++ b/ginac/clifford.h @@ -149,6 +149,10 @@ ex dirac_slash(const ex & e, const ex & dim, unsigned char rl = 0); * @param trONE Expression to be returned as the trace of the unit matrix */ ex dirac_trace(const ex & e, unsigned char rl = 0, const ex & trONE = 4); +/** Bring all products of clifford objects in an expression into a canonical + * order. This is not necessarily the most simple form but it will allow + * to checking two expressions for equality. */ +ex canonicalize_clifford(const ex & e); } // namespace GiNaC diff --git a/ginac/expairseq.cpp b/ginac/expairseq.cpp index 9f570617..b6a0f372 100644 --- a/ginac/expairseq.cpp +++ b/ginac/expairseq.cpp @@ -991,7 +991,6 @@ void expairseq::make_flat(const epvector &v) return; } - /** Brings this expairseq into a sorted (canonical) form. */ void expairseq::canonicalize(void) { diff --git a/ginac/indexed.cpp b/ginac/indexed.cpp index 3240c8dd..435d9527 100644 --- a/ginac/indexed.cpp +++ b/ginac/indexed.cpp @@ -21,6 +21,7 @@ */ #include +#include #include "indexed.h" #include "idx.h" @@ -253,31 +254,6 @@ int indexed::compare_same_type(const basic & other) const // reorder index pairs with known symmetry properties, while sort_index_vector() // always sorts the whole vector. -/** Bring a vector of indices into a canonic order (don't care about the - * symmetry of the objects carrying the indices). Dummy indices will lie - * next to each other after the sorting. - * - * @param v Index vector to be sorted */ -static void sort_index_vector(exvector &v) -{ - // Nothing to sort if less than 2 elements - if (v.size() < 2) - return; - - // Simple bubble sort algorithm should be sufficient for the small - // number of indices expected - exvector::iterator it1 = v.begin(), itend = v.end(), next_to_last_idx = itend - 1; - while (it1 != next_to_last_idx) { - exvector::iterator it2 = it1 + 1; - while (it2 != itend) { - if (it1->compare(*it2) > 0) - it1->swap(*it2); - it2++; - } - it1++; - } -} - /** Bring a vector of indices into a canonic order. This operation only makes * sense if the object carrying these indices is either symmetric or totally * antisymmetric with respect to the indices. @@ -570,9 +546,75 @@ exvector power::get_free_indices(void) const return basis.get_free_indices(); } +/* Function object for STL sort() */ +struct ex_is_less { + bool operator() (const ex &lh, const ex &rh) const + { + return lh.compare(rh) < 0; + } +}; + +/** Rename dummy indices in an expression. + * + * @param e Expression to be worked on + * @param local_dummy_indices The set of dummy indices that appear in the + * expression "e" + * @param global_dummy_indices The set of dummy indices that have appeared + * before and which we would like to use in "e", too. This gets updated + * by the function */ +static ex rename_dummy_indices(const ex & e, exvector & global_dummy_indices, exvector & local_dummy_indices) +{ + int global_size = global_dummy_indices.size(), + local_size = local_dummy_indices.size(); + + // Any local dummy indices at all? + if (local_size == 0) + return e; + + sort(local_dummy_indices.begin(), local_dummy_indices.end(), ex_is_less()); + + if (global_size < local_size) { + + // More local indices than we encountered before, add the new ones + // to the global set + int remaining = local_size - global_size; + exvector::const_iterator it = local_dummy_indices.begin(), itend = local_dummy_indices.end(); + while (it != itend && remaining > 0) { + exvector::const_iterator git = global_dummy_indices.begin(), gitend = global_dummy_indices.end(); + while (git != gitend) { + if (it->is_equal(*git)) + goto found; + git++; + } + global_dummy_indices.push_back(*it); + global_size++; + remaining--; +found: it++; + } + sort(global_dummy_indices.begin(), global_dummy_indices.end(), ex_is_less()); + } + + // Replace index symbols in expression + GINAC_ASSERT(local_size <= global_size); + bool all_equal = true; + lst local_syms, global_syms; + for (unsigned i=0; iget_free_indices(); un.insert(un.end(), free_indices_of_factor.begin(), free_indices_of_factor.end()); it1++; } - find_free_and_dummy(un, free_indices, dummy_indices); + find_free_and_dummy(un, free_indices, local_dummy_indices); ex r; if (something_changed) @@ -705,6 +747,9 @@ contraction_done: else r = e; + // 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)) @@ -714,17 +759,18 @@ contraction_done: } /** Simplify indexed expression, return list of free indices. */ -ex simplify_indexed(const ex & e, exvector & free_indices, const scalar_products & sp) +ex simplify_indexed(const ex & e, exvector & free_indices, exvector & dummy_indices, const scalar_products & sp) { // Expand the expression ex e_expanded = e.expand(); // Simplification of single indexed object: just find the free indices + // (and perform dummy index renaming if if (is_ex_of_type(e_expanded, indexed)) { const indexed &i = ex_to_indexed(e_expanded); - exvector dummy_indices; - find_free_and_dummy(i.seq.begin() + 1, i.seq.end(), free_indices, dummy_indices); - return e_expanded; + exvector local_dummy_indices; + find_free_and_dummy(i.seq.begin() + 1, i.seq.end(), free_indices, local_dummy_indices); + return rename_dummy_indices(e_expanded, dummy_indices, local_dummy_indices); } // Simplification of sum = sum of simplifications, check consistency of @@ -736,7 +782,7 @@ ex simplify_indexed(const ex & e, exvector & free_indices, const scalar_products for (unsigned i=0; isetflag(status_flags::dynallocated); diff --git a/ginac/ncmul.h b/ginac/ncmul.h index 05f99949..05b96da3 100644 --- a/ginac/ncmul.h +++ b/ginac/ncmul.h @@ -60,7 +60,6 @@ public: ex expand(unsigned options=0) const; ex coeff(const ex & s, int n=1) const; ex eval(int level=0) const; - ex subs(const lst & ls, const lst & lr) const; exvector get_free_indices(void) const; ex thisexprseq(const exvector & v) const; ex thisexprseq(exvector * vp) const; diff --git a/ginac/normal.cpp b/ginac/normal.cpp index 369e2c65..d3a407d5 100644 --- a/ginac/normal.cpp +++ b/ginac/normal.cpp @@ -36,7 +36,6 @@ #include "inifcns.h" #include "lst.h" #include "mul.h" -#include "ncmul.h" #include "numeric.h" #include "power.h" #include "relational.h"