X-Git-Url: https://ginac.de/ginac.git//ginac.git?a=blobdiff_plain;f=ginac%2Findexed.cpp;h=12561cf1d88f0557ef41c2f3e486a590361e80db;hb=09f37bdbd46f469b3a8a902a43d0f795c41a89bf;hp=c6a84e31ac70e6346091651bd0eb61f1a70fef5d;hpb=aa6281216091efd92dc5fcc3f96c7189114e80f1;p=ginac.git diff --git a/ginac/indexed.cpp b/ginac/indexed.cpp index c6a84e31..12561cf1 100644 --- a/ginac/indexed.cpp +++ b/ginac/indexed.cpp @@ -3,7 +3,7 @@ * Implementation of GiNaC's indexed expressions. */ /* - * GiNaC Copyright (C) 1999-2001 Johannes Gutenberg University Mainz, Germany + * GiNaC Copyright (C) 1999-2002 Johannes Gutenberg University Mainz, Germany * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by @@ -20,8 +20,8 @@ * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA */ +#include #include -#include #include "indexed.h" #include "idx.h" @@ -30,23 +30,22 @@ #include "ncmul.h" #include "power.h" #include "symmetry.h" +#include "operators.h" #include "lst.h" #include "print.h" #include "archive.h" #include "utils.h" -#include "debugmsg.h" namespace GiNaC { GINAC_IMPLEMENT_REGISTERED_CLASS(indexed, exprseq) ////////// -// default constructor, destructor, copy constructor assignment operator and helpers +// default ctor, dtor, copy ctor, assignment operator and helpers ////////// indexed::indexed() : symtree(sy_none()) { - debugmsg("indexed default constructor", LOGLEVEL_CONSTRUCT); tinfo_key = TINFO_indexed; } @@ -64,63 +63,54 @@ DEFAULT_DESTROY(indexed) indexed::indexed(const ex & b) : inherited(b), symtree(sy_none()) { - debugmsg("indexed constructor from ex", LOGLEVEL_CONSTRUCT); tinfo_key = TINFO_indexed; validate(); } indexed::indexed(const ex & b, const ex & i1) : inherited(b, i1), symtree(sy_none()) { - debugmsg("indexed constructor from ex,ex", LOGLEVEL_CONSTRUCT); tinfo_key = TINFO_indexed; validate(); } indexed::indexed(const ex & b, const ex & i1, const ex & i2) : inherited(b, i1, i2), symtree(sy_none()) { - debugmsg("indexed constructor from ex,ex,ex", LOGLEVEL_CONSTRUCT); tinfo_key = TINFO_indexed; validate(); } indexed::indexed(const ex & b, const ex & i1, const ex & i2, const ex & i3) : inherited(b, i1, i2, i3), symtree(sy_none()) { - debugmsg("indexed constructor from ex,ex,ex,ex", LOGLEVEL_CONSTRUCT); tinfo_key = TINFO_indexed; validate(); } indexed::indexed(const ex & b, const ex & i1, const ex & i2, const ex & i3, const ex & i4) : inherited(b, i1, i2, i3, i4), symtree(sy_none()) { - debugmsg("indexed constructor from ex,ex,ex,ex,ex", LOGLEVEL_CONSTRUCT); tinfo_key = TINFO_indexed; validate(); } indexed::indexed(const ex & b, const symmetry & symm, const ex & i1, const ex & i2) : inherited(b, i1, i2), symtree(symm) { - debugmsg("indexed constructor from ex,symmetry,ex,ex", LOGLEVEL_CONSTRUCT); tinfo_key = TINFO_indexed; validate(); } indexed::indexed(const ex & b, const symmetry & symm, const ex & i1, const ex & i2, const ex & i3) : inherited(b, i1, i2, i3), symtree(symm) { - debugmsg("indexed constructor from ex,symmetry,ex,ex,ex", LOGLEVEL_CONSTRUCT); tinfo_key = TINFO_indexed; validate(); } indexed::indexed(const ex & b, const symmetry & symm, const ex & i1, const ex & i2, const ex & i3, const ex & i4) : inherited(b, i1, i2, i3, i4), symtree(symm) { - debugmsg("indexed constructor from ex,symmetry,ex,ex,ex,ex", LOGLEVEL_CONSTRUCT); tinfo_key = TINFO_indexed; validate(); } indexed::indexed(const ex & b, const exvector & v) : inherited(b), symtree(sy_none()) { - debugmsg("indexed constructor from ex,exvector", LOGLEVEL_CONSTRUCT); seq.insert(seq.end(), v.begin(), v.end()); tinfo_key = TINFO_indexed; validate(); @@ -128,7 +118,6 @@ indexed::indexed(const ex & b, const exvector & v) : inherited(b), symtree(sy_no indexed::indexed(const ex & b, const symmetry & symm, const exvector & v) : inherited(b), symtree(symm) { - debugmsg("indexed constructor from ex,symmetry,exvector", LOGLEVEL_CONSTRUCT); seq.insert(seq.end(), v.begin(), v.end()); tinfo_key = TINFO_indexed; validate(); @@ -136,19 +125,16 @@ indexed::indexed(const ex & b, const symmetry & symm, const exvector & v) : inhe indexed::indexed(const symmetry & symm, const exprseq & es) : inherited(es), symtree(symm) { - debugmsg("indexed constructor from symmetry,exprseq", LOGLEVEL_CONSTRUCT); tinfo_key = TINFO_indexed; } indexed::indexed(const symmetry & symm, const exvector & v, bool discardable) : inherited(v, discardable), symtree(symm) { - debugmsg("indexed constructor from symmetry,exvector", LOGLEVEL_CONSTRUCT); tinfo_key = TINFO_indexed; } indexed::indexed(const symmetry & symm, exvector * vp) : inherited(vp), symtree(symm) { - debugmsg("indexed constructor from symmetry,exvector *", LOGLEVEL_CONSTRUCT); tinfo_key = TINFO_indexed; } @@ -158,7 +144,6 @@ indexed::indexed(const symmetry & symm, exvector * vp) : inherited(vp), symtree( indexed::indexed(const archive_node &n, const lst &sym_lst) : inherited(n, sym_lst) { - debugmsg("indexed constructor from archive_node", LOGLEVEL_CONSTRUCT); if (!n.find_ex("symmetry", symtree, sym_lst)) { // GiNaC versions <= 0.9.0 had an unsigned "symmetry" property unsigned symm = 0; @@ -174,7 +159,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(ex_to(symtree)).validate(seq.size() - 1); } } @@ -192,10 +177,9 @@ DEFAULT_UNARCHIVE(indexed) void indexed::print(const print_context & c, unsigned level) const { - debugmsg("indexed print", LOGLEVEL_PRINT); GINAC_ASSERT(seq.size() > 0); - if (is_of_type(c, print_tree)) { + if (is_a(c)) { c.s << std::string(level, ' ') << class_name() << std::hex << ", hash=0x" << hashvalue << ", flags=0x" << flags << std::dec @@ -207,21 +191,19 @@ void indexed::print(const print_context & c, unsigned level) const } else { - bool is_tex = is_of_type(c, print_latex); + bool is_tex = is_a(c); const ex & base = seq[0]; - bool need_parens = is_ex_exactly_of_type(base, add) || is_ex_exactly_of_type(base, mul) - || is_ex_exactly_of_type(base, ncmul) || is_ex_exactly_of_type(base, power) - || is_ex_of_type(base, indexed); + + if (precedence() <= level) + c.s << (is_tex ? "{(" : "("); if (is_tex) c.s << "{"; - if (need_parens) - c.s << "("; - base.print(c); - if (need_parens) - c.s << ")"; + base.print(c, precedence()); if (is_tex) c.s << "}"; printindices(c, level); + if (precedence() <= level) + c.s << (is_tex ? ")}" : ")"); } } @@ -250,7 +232,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(other)); return inherited::compare_same_type(other); } @@ -264,10 +246,10 @@ ex indexed::eval(int level) const // If the base object is 0, the whole object is 0 if (base.is_zero()) - return _ex0(); + return _ex0; // If the base object is a product, pull out the numeric factor - if (is_ex_exactly_of_type(base, mul) && is_ex_exactly_of_type(base.op(base.nops() - 1), numeric)) { + if (is_exactly_a(base) && is_exactly_a(base.op(base.nops() - 1))) { exvector v(seq); ex f = ex_to(base.op(base.nops() - 1)); v[0] = seq[0] / f; @@ -277,36 +259,18 @@ 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(symtree)); int sig = canonicalize(v.begin() + 1, ex_to(symtree)); if (sig != INT_MAX) { // Something has changed while sorting indices, more evaluations later if (sig == 0) - return _ex0(); + return _ex0; return ex(sig) * thisexprseq(v); } } // Let the class of the base object perform additional evaluations - return base.bp->eval_indexed(*this); -} - -int indexed::degree(const ex & s) const -{ - return is_equal(*s.bp) ? 1 : 0; -} - -int indexed::ldegree(const ex & s) const -{ - return is_equal(*s.bp) ? 1 : 0; -} - -ex indexed::coeff(const ex & s, int n) const -{ - if (is_equal(*s.bp)) - return n==1 ? _ex1() : _ex0(); - else - return n==0 ? ex(*this) : _ex0(); + return ex_to(base).eval_indexed(*this); } ex indexed::thisexprseq(const exvector & v) const @@ -323,11 +287,11 @@ ex indexed::expand(unsigned options) const { GINAC_ASSERT(seq.size() > 0); - if ((options & expand_options::expand_indexed) && is_ex_exactly_of_type(seq[0], add)) { + if ((options & expand_options::expand_indexed) && is_exactly_a(seq[0])) { // expand_indexed expands (a+b).i -> a.i + b.i const ex & base = seq[0]; - ex sum = _ex0(); + ex sum = _ex0; for (unsigned i=0; i(c)) { // TeX output: group by variance bool first = true; bool covariant = true; while (it != itend) { - bool cur_covariant = (is_ex_of_type(*it, varidx) ? ex_to(*it).is_covariant() : true); - if (first || cur_covariant != covariant) { + bool cur_covariant = (is_a(*it) ? ex_to(*it).is_covariant() : true); + if (first || cur_covariant != covariant) { // Variance changed + // The empty {} prevents indices from ending up on top of each other if (!first) - c.s << "}"; + c.s << "}{}"; covariant = cur_covariant; if (covariant) c.s << "_{"; @@ -398,15 +363,15 @@ void indexed::validate(void) const GINAC_ASSERT(seq.size() > 0); exvector::const_iterator it = seq.begin() + 1, itend = seq.end(); while (it != itend) { - if (!is_ex_of_type(*it, idx)) + if (!is_a(*it)) throw(std::invalid_argument("indices of indexed object must be of type idx")); it++; } if (!symtree.is_zero()) { - if (!is_ex_exactly_of_type(symtree, symmetry)) + if (!is_exactly_a(symtree)) throw(std::invalid_argument("symmetry of indexed object must be of type symmetry")); - ex_to_nonconst_symmetry(symtree).validate(seq.size() - 1); + const_cast(ex_to(symtree)).validate(seq.size() - 1); } } @@ -415,7 +380,7 @@ void indexed::validate(void) const * @see ex::diff */ ex indexed::derivative(const symbol & s) const { - return _ex0(); + return _ex0; } ////////// @@ -545,6 +510,7 @@ static ex rename_dummy_indices(const ex & e, exvector & global_dummy_indices, ex // More local indices than we encountered before, add the new ones // to the global set + int old_global_size = global_size; int remaining = local_size - global_size; exvector::const_iterator it = local_dummy_indices.begin(), itend = local_dummy_indices.end(); while (it != itend && remaining > 0) { @@ -555,26 +521,35 @@ static ex rename_dummy_indices(const ex & e, exvector & global_dummy_indices, ex } it++; } - } - // Replace index symbols in expression - GINAC_ASSERT(local_size <= global_size); - bool all_equal = true; - lst local_syms, global_syms; - for (unsigned i=0; i(local_dummy_indices[i]).get_dim().is_equal(ex_to(global_dummy_indices[i]).get_dim())) { - all_equal = false; - local_syms.append(loc_sym); - global_syms.append(glob_sym); - } + // If this is the first set of local indices, do nothing + if (old_global_size == 0) + return e; } - if (all_equal) + GINAC_ASSERT(local_size <= global_size); + + // Construct lists of index symbols + exlist local_syms, global_syms; + for (unsigned i=0; i(local_uniq), ex_is_less()); + set_difference(global_syms.begin(), global_syms.end(), local_syms.begin(), local_syms.end(), std::back_insert_iterator(global_uniq), ex_is_less()); + + // Replace remaining non-common local index symbols by global ones + if (local_uniq.empty()) return e; - else - return e.subs(local_syms, global_syms); + else { + while (global_uniq.size() > local_uniq.size()) + global_uniq.pop_back(); + return e.subs(lst(local_uniq), lst(global_uniq)); + } } /** Simplify product of indexed expressions (commutative, noncommutative and @@ -583,24 +558,24 @@ ex simplify_indexed_product(const ex & e, exvector & free_indices, exvector & du { // Remember whether the product was commutative or noncommutative // (because we chop it into factors and need to reassemble later) - bool non_commutative = is_ex_exactly_of_type(e, ncmul); + bool non_commutative = is_exactly_a(e); // Collect factors in an exvector, store squares twice exvector v; v.reserve(e.nops() * 2); - if (is_ex_exactly_of_type(e, power)) { + if (is_exactly_a(e)) { // We only get called for simple squares, split a^2 -> a*a - GINAC_ASSERT(e.op(1).is_equal(_ex2())); + GINAC_ASSERT(e.op(1).is_equal(_ex2)); v.push_back(e.op(0)); v.push_back(e.op(0)); } else { for (unsigned i=0; i(f) && f.op(1).is_equal(_ex2)) { v.push_back(f.op(0)); v.push_back(f.op(0)); - } else if (is_ex_exactly_of_type(f, ncmul)) { + } else if (is_exactly_a(f)) { // Noncommutative factor found, split it as well non_commutative = true; // everything becomes noncommutative, ncmul will sort out the commutative factors later for (unsigned j=0; j(*it1)) continue; bool first_noncommutative = (it1->return_type() != return_types::commutative); @@ -630,7 +605,7 @@ try_again: exvector::iterator it2; for (it2 = it1 + 1; it2 != itend; it2++) { - if (!is_ex_of_type(*it2, indexed)) + if (!is_a(*it2)) continue; bool second_noncommutative = (it2->return_type() != return_types::commutative); @@ -653,48 +628,25 @@ try_again: if (free.empty()) { if (sp.is_defined(*it1, *it2)) { *it1 = sp.evaluate(*it1, *it2); - *it2 = _ex1(); + *it2 = _ex1; goto contraction_done; } } - // Contraction of symmetric with antisymmetric object is zero - if (num_dummies > 1 - && ex_to(ex_to(*it1).symtree).has_symmetry() - && ex_to(ex_to(*it2).symtree).has_symmetry()) { - - // Check all pairs of dummy indices - for (unsigned idx1=0; idx1subs(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(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(it2->op(0)).contract_with(it2, it1, v); } if (contracted) { contraction_done: if (first_noncommutative || second_noncommutative - || is_ex_exactly_of_type(*it1, add) || is_ex_exactly_of_type(*it2, add) - || is_ex_exactly_of_type(*it1, mul) || is_ex_exactly_of_type(*it2, mul) - || is_ex_exactly_of_type(*it1, ncmul) || is_ex_exactly_of_type(*it2, ncmul)) { + || is_exactly_a(*it1) || is_exactly_a(*it2) + || is_exactly_a(*it1) || is_exactly_a(*it2) + || is_exactly_a(*it1) || is_exactly_a(*it2)) { // One of the factors became a sum or product: // re-expand expression and run again @@ -720,7 +672,7 @@ contraction_done: it1 = v.begin(); itend = v.end(); while (it1 != itend) { exvector free_indices_of_factor; - if (is_ex_of_type(*it1, indexed)) { + if (is_a(*it1)) { exvector dummy_indices_of_factor; find_free_and_dummy(ex_to(*it1).seq.begin() + 1, ex_to(*it1).seq.end(), free_indices_of_factor, dummy_indices_of_factor); individual_dummy_indices.insert(individual_dummy_indices.end(), dummy_indices_of_factor.begin(), dummy_indices_of_factor.end()); @@ -739,13 +691,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; iscalar_mul_indexed(r.op(0), ex_to(r.op(1))); + if (is_exactly_a(r) && r.nops() == 2 + && is_exactly_a(r.op(1)) && is_a(r.op(0))) + return ex_to(r.op(0).op(0)).scalar_mul_indexed(r.op(0), ex_to(r.op(1))); else return r; } @@ -758,7 +723,7 @@ ex simplify_indexed(const ex & e, exvector & free_indices, exvector & dummy_indi // Simplification of single indexed object: just find the free indices // and perform dummy index renaming - if (is_ex_of_type(e_expanded, indexed)) { + if (is_a(e_expanded)) { const indexed &i = ex_to(e_expanded); exvector local_dummy_indices; find_free_and_dummy(i.seq.begin() + 1, i.seq.end(), free_indices, local_dummy_indices); @@ -767,9 +732,9 @@ ex simplify_indexed(const ex & e, exvector & free_indices, exvector & dummy_indi // Simplification of sum = sum of simplifications, check consistency of // free indices in each term - if (is_ex_exactly_of_type(e_expanded, add)) { + if (is_exactly_a(e_expanded)) { bool first = true; - ex sum = _ex0(); + ex sum = _ex0; free_indices.clear(); for (unsigned i=0; iadd_indexed(sum, term); + if (is_a(sum) && is_a(term)) + sum = ex_to(sum.op(0)).add_indexed(sum, term); else sum += term; } @@ -795,9 +760,9 @@ ex simplify_indexed(const ex & e, exvector & free_indices, exvector & dummy_indi } // Simplification of products - if (is_ex_exactly_of_type(e_expanded, mul) - || is_ex_exactly_of_type(e_expanded, ncmul) - || (is_ex_exactly_of_type(e_expanded, power) && is_ex_of_type(e_expanded.op(0), indexed) && e_expanded.op(1).is_equal(_ex2()))) + if (is_exactly_a(e_expanded) + || is_exactly_a(e_expanded) + || (is_exactly_a(e_expanded) && is_a(e_expanded.op(0)) && e_expanded.op(1).is_equal(_ex2))) return simplify_indexed_product(e_expanded, free_indices, dummy_indices, sp); // Cannot do anything @@ -903,8 +868,8 @@ void scalar_products::debugprint(void) const spmapkey scalar_products::make_key(const ex & v1, const ex & v2) { // If indexed, extract base objects - ex s1 = is_ex_of_type(v1, indexed) ? v1.op(0) : v1; - ex s2 = is_ex_of_type(v2, indexed) ? v2.op(0) : v2; + ex s1 = is_a(v1) ? v1.op(0) : v1; + ex s2 = is_a(v2) ? v2.op(0) : v2; // Enforce canonical order in pair if (s1.compare(s2) > 0)