X-Git-Url: https://www.ginac.de/ginac.git//ginac.git?p=ginac.git;a=blobdiff_plain;f=ginac%2Findexed.cpp;h=0d74baf55ed92aa7f61e12ec03693fc9c82f0976;hp=d39692c4e859d2a8be1f381052caeba203551996;hb=5432feb4dae4c48f3b652d661aa323f07d20c00b;hpb=22d004fe5759658a8232736e40b49202987e9407 diff --git a/ginac/indexed.cpp b/ginac/indexed.cpp index d39692c4..0d74baf5 100644 --- a/ginac/indexed.cpp +++ b/ginac/indexed.cpp @@ -3,7 +3,7 @@ * Implementation of GiNaC's indexed expressions. */ /* - * GiNaC Copyright (C) 1999-2005 Johannes Gutenberg University Mainz, Germany + * GiNaC Copyright (C) 1999-2010 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,10 +20,6 @@ * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA */ -#include -#include -#include - #include "indexed.h" #include "idx.h" #include "add.h" @@ -39,6 +35,12 @@ #include "utils.h" #include "integral.h" #include "matrix.h" +#include "inifcns.h" + +#include +#include +#include +#include namespace GiNaC { @@ -53,7 +55,6 @@ GINAC_IMPLEMENT_REGISTERED_CLASS_OPT(indexed, exprseq, indexed::indexed() : symtree(not_symmetric()) { - tinfo_key = TINFO_indexed; } ////////// @@ -62,87 +63,75 @@ indexed::indexed() : symtree(not_symmetric()) indexed::indexed(const ex & b) : inherited(b), symtree(not_symmetric()) { - tinfo_key = TINFO_indexed; validate(); } indexed::indexed(const ex & b, const ex & i1) : inherited(b, i1), symtree(not_symmetric()) { - tinfo_key = TINFO_indexed; validate(); } indexed::indexed(const ex & b, const ex & i1, const ex & i2) : inherited(b, i1, i2), symtree(not_symmetric()) { - 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(not_symmetric()) { - 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(not_symmetric()) { - 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) { - 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) { - 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) { - tinfo_key = TINFO_indexed; validate(); } indexed::indexed(const ex & b, const exvector & v) : inherited(b), symtree(not_symmetric()) { seq.insert(seq.end(), v.begin(), v.end()); - tinfo_key = TINFO_indexed; validate(); } indexed::indexed(const ex & b, const symmetry & symm, const exvector & v) : inherited(b), symtree(symm) { seq.insert(seq.end(), v.begin(), v.end()); - tinfo_key = TINFO_indexed; validate(); } indexed::indexed(const symmetry & symm, const exprseq & es) : inherited(es), symtree(symm) { - tinfo_key = TINFO_indexed; } indexed::indexed(const symmetry & symm, const exvector & v, bool discardable) : inherited(v, discardable), symtree(symm) { - tinfo_key = TINFO_indexed; } indexed::indexed(const symmetry & symm, std::auto_ptr vp) : inherited(vp), symtree(symm) { - tinfo_key = TINFO_indexed; } ////////// // archiving ////////// -indexed::indexed(const archive_node &n, lst &sym_lst) : inherited(n, sym_lst) +void indexed::read_archive(const archive_node &n, lst &sym_lst) { + inherited::read_archive(n, sym_lst); if (!n.find_ex("symmetry", symtree, sym_lst)) { // GiNaC versions <= 0.9.0 had an unsigned "symmetry" property unsigned symm = 0; @@ -161,6 +150,7 @@ indexed::indexed(const archive_node &n, lst &sym_lst) : inherited(n, sym_lst) const_cast(ex_to(symtree)).validate(seq.size() - 1); } } +GINAC_BIND_UNARCHIVER(indexed); void indexed::archive(archive_node &n) const { @@ -168,8 +158,6 @@ void indexed::archive(archive_node &n) const n.add_ex("symmetry", symtree); } -DEFAULT_UNARCHIVE(indexed) - ////////// // functions overriding virtual functions from base classes ////////// @@ -297,7 +285,7 @@ ex indexed::eval(int level) const return f * thiscontainer(v); } - if(this->tinfo()==TINFO_indexed && seq.size()==1) + if((typeid(*this) == typeid(indexed)) && seq.size()==1) return base; // Canonicalize indices according to the symmetry properties @@ -305,7 +293,7 @@ ex indexed::eval(int level) const exvector v = seq; GINAC_ASSERT(is_exactly_a(symtree)); int sig = canonicalize(v.begin() + 1, ex_to(symtree)); - if (sig != INT_MAX) { + if (sig != std::numeric_limits::max()) { // Something has changed while sorting indices, more evaluations later if (sig == 0) return _ex0; @@ -317,6 +305,20 @@ ex indexed::eval(int level) const return ex_to(base).eval_indexed(*this); } +ex indexed::real_part() const +{ + if(op(0).info(info_flags::real)) + return *this; + return real_part_function(*this).hold(); +} + +ex indexed::imag_part() const +{ + if(op(0).info(info_flags::real)) + return 0; + return imag_part_function(*this).hold(); +} + ex indexed::thiscontainer(const exvector & v) const { return indexed(ex_to(symtree), v); @@ -624,10 +626,60 @@ bool reposition_dummy_indices(ex & e, exvector & variant_dummy_indices, exvector { bool something_changed = false; + // Find dummy symbols that occur twice in the same indexed object. + exvector local_var_dummies; + local_var_dummies.reserve(e.nops()/2); + for (size_t i=1; i(e.op(i))) + continue; + for (size_t j=i+1; jop(0)) { + variant_dummy_indices.erase(k); + break; + } + } + break; + } + } + } + + // In the case where a dummy symbol occurs twice in the same indexed object + // we try all posibilities of raising/lowering and keep the least one in + // the sense of ex_is_less. + ex optimal_e = e; + size_t numpossibs = 1 << local_var_dummies.size(); + for (size_t i=0; i(curr_idx).toggle_variance(); + m[curr_idx] = curr_toggle; + m[curr_toggle] = curr_idx; + } + try_e = e.subs(m, subs_options::no_pattern); + } + if(ex_is_less()(try_e, optimal_e)) + { optimal_e = try_e; + something_changed = true; + } + } + e = optimal_e; + + if (!is_a(e)) + return true; + + exvector seq = ex_to(e).seq; + // If a dummy index is encountered for the first time in the // product, pull it up, otherwise, pull it down - exvector::const_iterator it2, it2start, it2end; - for (it2start = ex_to(e).seq.begin(), it2end = ex_to(e).seq.end(), it2 = it2start + 1; it2 != it2end; ++it2) { + for (exvector::iterator it2 = seq.begin()+1, it2end = seq.end(); + it2 != it2end; ++it2) { if (!is_exactly_a(*it2)) continue; @@ -635,14 +687,20 @@ bool reposition_dummy_indices(ex & e, exvector & variant_dummy_indices, exvector for (vit = variant_dummy_indices.begin(), vitend = variant_dummy_indices.end(); vit != vitend; ++vit) { if (it2->op(0).is_equal(vit->op(0))) { if (ex_to(*it2).is_covariant()) { - e = e.subs(lst( - *it2 == ex_to(*it2).toggle_variance(), - ex_to(*it2).toggle_variance() == *it2 - ), subs_options::no_pattern); + /* + * N.B. we don't want to use + * + * e = e.subs(lst( + * *it2 == ex_to(*it2).toggle_variance(), + * ex_to(*it2).toggle_variance() == *it2 + * ), subs_options::no_pattern); + * + * since this can trigger non-trivial repositioning of indices, + * e.g. due to non-trivial symmetry properties of e, thus + * invalidating iterators + */ + *it2 = ex_to(*it2).toggle_variance(); something_changed = true; - it2 = ex_to(e).seq.begin() + (it2 - it2start); - it2start = ex_to(e).seq.begin(); - it2end = ex_to(e).seq.end(); } moved_indices.push_back(*vit); variant_dummy_indices.erase(vit); @@ -653,11 +711,8 @@ bool reposition_dummy_indices(ex & e, exvector & variant_dummy_indices, exvector for (vit = moved_indices.begin(), vitend = moved_indices.end(); vit != vitend; ++vit) { if (it2->op(0).is_equal(vit->op(0))) { if (ex_to(*it2).is_contravariant()) { - e = e.subs(*it2 == ex_to(*it2).toggle_variance(), subs_options::no_pattern); + *it2 = ex_to(*it2).toggle_variance(); something_changed = true; - it2 = ex_to(e).seq.begin() + (it2 - it2start); - it2start = ex_to(e).seq.begin(); - it2end = ex_to(e).seq.end(); } goto next_index; } @@ -666,6 +721,9 @@ bool reposition_dummy_indices(ex & e, exvector & variant_dummy_indices, exvector next_index: ; } + if (something_changed) + e = ex_to(e).thiscontainer(seq); + return something_changed; } @@ -1337,7 +1395,7 @@ exvector get_all_dummy_indices_safely(const ex & e) else if (is_a(e) || is_a(e)) { exvector dummies; exvector free_indices; - for (int i=0; i(e)) { exvector result; - for(int i=0; i fcn(expand_dummy_sum, subs_idx); if (is_a(e_expanded) || is_a(e_expanded) || is_a(e_expanded)) { return e_expanded.map(fcn); - } else if (is_a(e_expanded) || is_a(e_expanded) || is_a(e_expanded)) { - exvector v = get_all_dummy_indices(e_expanded); - exvector::const_iterator it = v.begin(), itend = v.end(); - while (it != itend) { - varidx nu = ex_to(*it); - if (nu.is_dim_numeric()) { - ex en = 0; - for (int i=0; i < ex_to(nu.get_dim()).to_int(); i++) { - if (is_a(nu) && !subs_idx) { - en += e_expanded.subs(lst(nu == varidx(i, nu.get_dim(), true), nu.toggle_variance() == varidx(i, nu.get_dim()))); - } else { - en += e_expanded.subs(lst(nu == idx(i, nu.get_dim()), nu.toggle_variance() == idx(i, nu.get_dim()))); - } - } - return expand_dummy_sum(en, subs_idx); - } - ++it; - } - return e; - } else if (is_a(e_expanded)) { - exvector v = ex_to(e_expanded).get_dummy_indices(); - exvector::const_iterator it = v.begin(), itend = v.end(); - while (it != itend) { - varidx nu = ex_to(*it); - if (nu.is_dim_numeric()) { + } else if (is_a(e_expanded) || is_a(e_expanded) || is_a(e_expanded) || is_a(e_expanded)) { + exvector v; + if (is_a(e_expanded)) + v = ex_to(e_expanded).get_dummy_indices(); + else + v = get_all_dummy_indices(e_expanded); + ex result = e_expanded; + for(exvector::const_iterator it=v.begin(); it!=v.end(); ++it) { + ex nu = *it; + if (ex_to(nu).get_dim().info(info_flags::nonnegint)) { + int idim = ex_to(ex_to(nu).get_dim()).to_int(); ex en = 0; - for (int i=0; i < ex_to(nu.get_dim()).to_int(); i++) { - if (is_a(nu) && !subs_idx) { - en += e_expanded.subs(lst(nu == varidx(i, nu.get_dim(), true), nu.toggle_variance() == varidx(i, nu.get_dim()))); + for (int i=0; i < idim; i++) { + if (subs_idx && is_a(nu)) { + ex other = ex_to(nu).toggle_variance(); + en += result.subs(lst( + nu == idx(i, idim), + other == idx(i, idim) + )); } else { - en += e_expanded.subs(lst(nu == idx(i, nu.get_dim()), nu.toggle_variance() == idx(i, nu.get_dim()))); + en += result.subs( nu.op(0) == i ); } } - return expand_dummy_sum(en, subs_idx); - } - ++it; + result = en; + } } - return e; + return result; } else { return e; }