- Instead of just totally symmetric or antisymmetric, complex symmetries
authorChristian Bauer <Christian.Bauer@uni-mainz.de>
Mon, 11 Jun 2001 23:48:43 +0000 (23:48 +0000)
committerChristian Bauer <Christian.Bauer@uni-mainz.de>
Mon, 11 Jun 2001 23:48:43 +0000 (23:48 +0000)
  can now be defined for indexed objects. Symmetries are described by a
  tree of "symmetry" objects that is constructed with the sy_none(),
  sy_symm(), sy_anti() and sy_cycl() functions. The symmetry of a function
  with respect to its arguments can also be defined (this is currently
  only used for the Beta function).
- color_trace() and dirac_trace() can be applied to a more general class
  of expressions, e.g. using it on a relation will take the trace on both
  sides etc.
- Generalized map() to take a function object instead of a function pointer.
  This allows passing an arbitrary number of additional state to the
  function being called.
- The unarchiving functinos find_bool(), find_unsigned() and find_string()
  can take an additional "index" argument.

29 files changed:
ginac/Makefile.am
ginac/archive.cpp
ginac/archive.h
ginac/basic.cpp
ginac/basic.h
ginac/clifford.cpp
ginac/clifford.h
ginac/color.cpp
ginac/container.pl
ginac/ex.h
ginac/expairseq.cpp
ginac/expairseq.h
ginac/function.pl
ginac/ginac.h
ginac/idx.cpp
ginac/indexed.cpp
ginac/indexed.h
ginac/inifcns.cpp
ginac/inifcns.h
ginac/inifcns_gamma.cpp
ginac/power.cpp
ginac/power.h
ginac/pseries.cpp
ginac/relational.cpp
ginac/relational.h
ginac/tensor.cpp
ginac/tinfos.h
ginac/utils.h
ginac/wildcard.h

index 595e1f5..913457c 100644 (file)
@@ -8,7 +8,7 @@ libginac_la_SOURCES = add.cpp archive.cpp basic.cpp constant.cpp ex.cpp \
   symbol.cpp pseries.cpp utils.cpp ncmul.cpp structure.cpp exprseq_suppl.cpp \
   lst.cpp lst_suppl.cpp input_parser.yy input_lexer.ll input_lexer.h \
   remember.h remember.cpp debugmsg.h utils.h idx.cpp indexed.cpp tensor.cpp \
-  color.cpp clifford.cpp wildcard.cpp
+  color.cpp clifford.cpp wildcard.cpp symmetry.cpp
 libginac_la_LDFLAGS = -version-info $(LT_CURRENT):$(LT_REVISION):$(LT_AGE) \
   -release $(LT_RELEASE)
 ginacincludedir = $(includedir)/ginac
@@ -16,7 +16,8 @@ ginacinclude_HEADERS = ginac.h add.h archive.h basic.h constant.h ex.h \
   expair.h expairseq.h exprseq.h fail.h flags.h function.h inifcns.h \
   lst.h matrix.h mul.h ncmul.h normal.h numeric.h operators.h power.h \
   registrar.h relational.h pseries.h structure.h symbol.h tinfos.h assertion.h \
-  version.h idx.h indexed.h tensor.h color.h clifford.h wildcard.h print.h
+  version.h idx.h indexed.h tensor.h color.h clifford.h wildcard.h print.h \
+  symmetry.h
 LFLAGS = -Pginac_yy -olex.yy.c
 YFLAGS = -p ginac_yy -d
 EXTRA_DIST = container.pl function.pl structure.pl input_parser.h version.h.in
index cdd7619..2292228 100644 (file)
@@ -92,7 +92,7 @@ found:
        return nodes[i->root].unarchive(sym_lst);
 }
 
-ex archive::unarchive_ex(const lst &sym_lst, unsigned int index) const
+ex archive::unarchive_ex(const lst &sym_lst, unsigned index) const
 {
        if (index >= exprs.size())
                throw (std::range_error("index of archived expression out of range"));
@@ -101,7 +101,7 @@ ex archive::unarchive_ex(const lst &sym_lst, unsigned int index) const
        return nodes[exprs[index].root].unarchive(sym_lst);
 }
 
-ex archive::unarchive_ex(const lst &sym_lst, std::string &name, unsigned int index) const
+ex archive::unarchive_ex(const lst &sym_lst, std::string &name, unsigned index) const
 {
        if (index >= exprs.size())
                throw (std::range_error("index of archived expression out of range"));
@@ -113,12 +113,12 @@ ex archive::unarchive_ex(const lst &sym_lst, std::string &name, unsigned int ind
        return nodes[exprs[index].root].unarchive(sym_lst);
 }
 
-unsigned int archive::num_expressions(void) const
+unsigned archive::num_expressions(void) const
 {
        return exprs.size();
 }
 
-const archive_node &archive::get_top_node(unsigned int index) const
+const archive_node &archive::get_top_node(unsigned index) const
 {
        if (index >= exprs.size())
                throw (std::range_error("index of archived expression out of range"));
@@ -167,7 +167,7 @@ const archive_node &archive::get_top_node(unsigned int index) const
  */
 
 /** Write unsigned integer quantity to stream. */
-static void write_unsigned(std::ostream &os, unsigned int val)
+static void write_unsigned(std::ostream &os, unsigned val)
 {
        while (val >= 0x80) {
                os.put((val & 0x7f) | 0x80);
@@ -177,11 +177,11 @@ static void write_unsigned(std::ostream &os, unsigned int val)
 }
 
 /** Read unsigned integer quantity from stream. */
-static unsigned int read_unsigned(std::istream &is)
+static unsigned read_unsigned(std::istream &is)
 {
        unsigned char b;
-       unsigned int ret = 0;
-       unsigned int shift = 0;
+       unsigned ret = 0;
+       unsigned shift = 0;
        do {
         char b2;
                is.get(b2);
@@ -196,9 +196,9 @@ static unsigned int read_unsigned(std::istream &is)
 std::ostream &operator<<(std::ostream &os, const archive_node &n)
 {
        // Write properties
-       unsigned int num_props = n.props.size();
+       unsigned num_props = n.props.size();
        write_unsigned(os, num_props);
-       for (unsigned int i=0; i<num_props; i++) {
+       for (unsigned i=0; i<num_props; i++) {
                write_unsigned(os, n.props[i].type | (n.props[i].name << 3));
                write_unsigned(os, n.props[i].value);
        }
@@ -216,23 +216,23 @@ std::ostream &operator<<(std::ostream &os, const archive &ar)
        write_unsigned(os, ARCHIVE_VERSION);
 
        // Write atoms
-       unsigned int num_atoms = ar.atoms.size();
+       unsigned num_atoms = ar.atoms.size();
        write_unsigned(os, num_atoms);
-       for (unsigned int i=0; i<num_atoms; i++)
+       for (unsigned i=0; i<num_atoms; i++)
                os << ar.atoms[i] << std::ends;
 
        // Write expressions
-       unsigned int num_exprs = ar.exprs.size();
+       unsigned num_exprs = ar.exprs.size();
        write_unsigned(os, num_exprs);
-       for (unsigned int i=0; i<num_exprs; i++) {
+       for (unsigned i=0; i<num_exprs; i++) {
                write_unsigned(os, ar.exprs[i].name);
                write_unsigned(os, ar.exprs[i].root);
        }
 
        // Write nodes
-       unsigned int num_nodes = ar.nodes.size();
+       unsigned num_nodes = ar.nodes.size();
        write_unsigned(os, num_nodes);
-       for (unsigned int i=0; i<num_nodes; i++)
+       for (unsigned i=0; i<num_nodes; i++)
                os << ar.nodes[i];
        return os;
 }
@@ -241,10 +241,10 @@ std::ostream &operator<<(std::ostream &os, const archive &ar)
 std::istream &operator>>(std::istream &is, archive_node &n)
 {
        // Read properties
-       unsigned int num_props = read_unsigned(is);
+       unsigned num_props = read_unsigned(is);
        n.props.resize(num_props);
-       for (unsigned int i=0; i<num_props; i++) {
-               unsigned int name_type = read_unsigned(is);
+       for (unsigned i=0; i<num_props; i++) {
+               unsigned name_type = read_unsigned(is);
                n.props[i].type = (archive_node::property_type)(name_type & 7);
                n.props[i].name = name_type >> 3;
                n.props[i].value = read_unsigned(is);
@@ -260,29 +260,29 @@ std::istream &operator>>(std::istream &is, archive &ar)
        is.get(c1); is.get(c2); is.get(c3); is.get(c4);
        if (c1 != 'G' || c2 != 'A' || c3 != 'R' || c4 != 'C')
                throw (std::runtime_error("not a GiNaC archive (signature not found)"));
-       unsigned int version = read_unsigned(is);
+       unsigned version = read_unsigned(is);
        if (version > ARCHIVE_VERSION || version < ARCHIVE_VERSION - ARCHIVE_AGE)
                throw (std::runtime_error("archive version " + ToString(version) + " cannot be read by this GiNaC library (which supports versions " + ToString(ARCHIVE_VERSION-ARCHIVE_AGE) + " thru " + ToString(ARCHIVE_VERSION)));
 
        // Read atoms
-       unsigned int num_atoms = read_unsigned(is);
+       unsigned num_atoms = read_unsigned(is);
        ar.atoms.resize(num_atoms);
-       for (unsigned int i=0; i<num_atoms; i++)
+       for (unsigned i=0; i<num_atoms; i++)
                getline(is, ar.atoms[i], '\0');
 
        // Read expressions
-       unsigned int num_exprs = read_unsigned(is);
+       unsigned num_exprs = read_unsigned(is);
        ar.exprs.resize(num_exprs);
-       for (unsigned int i=0; i<num_exprs; i++) {
+       for (unsigned i=0; i<num_exprs; i++) {
                archive_atom name = read_unsigned(is);
                archive_node_id root = read_unsigned(is);
                ar.exprs[i] = archive::archived_ex(name, root);
        }
 
        // Read nodes
-       unsigned int num_nodes = read_unsigned(is);
+       unsigned num_nodes = read_unsigned(is);
        ar.nodes.resize(num_nodes, ar);
-       for (unsigned int i=0; i<num_nodes; i++)
+       for (unsigned i=0; i<num_nodes; i++)
                is >> ar.nodes[i];
        return is;
 }
@@ -360,7 +360,7 @@ void archive_node::add_bool(const std::string &name, bool value)
        props.push_back(property(a.atomize(name), PTYPE_BOOL, value));
 }
 
-void archive_node::add_unsigned(const std::string &name, unsigned int value)
+void archive_node::add_unsigned(const std::string &name, unsigned value)
 {
        props.push_back(property(a.atomize(name), PTYPE_UNSIGNED, value));
 }
@@ -378,73 +378,83 @@ void archive_node::add_ex(const std::string &name, const ex &value)
 }
 
 
-bool archive_node::find_bool(const std::string &name, bool &ret) const
+bool archive_node::find_bool(const std::string &name, bool &ret, unsigned index) const
 {
        archive_atom name_atom = a.atomize(name);
        std::vector<property>::const_iterator i = props.begin(), iend = props.end();
+       unsigned found_index = 0;
        while (i != iend) {
                if (i->type == PTYPE_BOOL && i->name == name_atom) {
-                       ret = i->value;
-                       return true;
+                       if (found_index == index) {
+                               ret = i->value;
+                               return true;
+                       }
+                       found_index++;
                }
                i++;
        }
        return false;
 }
 
-bool archive_node::find_unsigned(const std::string &name, unsigned int &ret) const
+bool archive_node::find_unsigned(const std::string &name, unsigned &ret, unsigned index) const
 {
        archive_atom name_atom = a.atomize(name);
        std::vector<property>::const_iterator i = props.begin(), iend = props.end();
+       unsigned found_index = 0;
        while (i != iend) {
                if (i->type == PTYPE_UNSIGNED && i->name == name_atom) {
-                       ret = i->value;
-                       return true;
+                       if (found_index == index) {
+                               ret = i->value;
+                               return true;
+                       }
+                       found_index++;
                }
                i++;
        }
        return false;
 }
 
-bool archive_node::find_string(const std::string &name, std::string &ret) const
+bool archive_node::find_string(const std::string &name, std::string &ret, unsigned index) const
 {
        archive_atom name_atom = a.atomize(name);
        std::vector<property>::const_iterator i = props.begin(), iend = props.end();
+       unsigned found_index = 0;
        while (i != iend) {
                if (i->type == PTYPE_STRING && i->name == name_atom) {
-                       ret = a.unatomize(i->value);
-                       return true;
+                       if (found_index == index) {
+                               ret = a.unatomize(i->value);
+                               return true;
+                       }
+                       found_index++;
                }
                i++;
        }
        return false;
 }
 
-bool archive_node::find_ex(const std::string &name, ex &ret, const lst &sym_lst, unsigned int index) const
+bool archive_node::find_ex(const std::string &name, ex &ret, const lst &sym_lst, unsigned index) const
 {
        archive_atom name_atom = a.atomize(name);
        std::vector<property>::const_iterator i = props.begin(), iend = props.end();
-       unsigned int found_index = 0;
+       unsigned found_index = 0;
        while (i != iend) {
                if (i->type == PTYPE_NODE && i->name == name_atom) {
-                       if (found_index == index)
-                               goto found;
+                       if (found_index == index) {
+                               ret = a.get_node(i->value).unarchive(sym_lst);
+                               return true;
+                       }
                        found_index++;
                }
                i++;
        }
        return false;
-
-found:
-       ret = a.get_node(i->value).unarchive(sym_lst);
-       return true;
 }
 
-const archive_node &archive_node::find_ex_node(const std::string &name, unsigned int index) const
+const archive_node &archive_node::find_ex_node(const std::string &name, unsigned index) const
 {
        archive_atom name_atom = a.atomize(name);
        std::vector<property>::const_iterator i = props.begin(), iend = props.end();
-       unsigned int found_index = 0;
+       unsigned found_index = 0;
        while (i != iend) {
                if (i->type == PTYPE_NODE && i->name == name_atom) {
                        if (found_index == index)
@@ -570,7 +580,7 @@ void archive::printraw(std::ostream &os) const
        os << "Expressions:\n";
        {
                std::vector<archived_ex>::const_iterator i = exprs.begin(), iend = exprs.end();
-               unsigned int index = 0;
+               unsigned index = 0;
                while (i != iend) {
                        os << " " << index << " \"" << unatomize(i->name) << "\" root node " << i->root << std::endl;
                        i++; index++;
index b2b403b..012086b 100644 (file)
@@ -36,10 +36,10 @@ class archive;
 
 
 /** Numerical ID value to refer to an archive_node. */
-typedef unsigned int archive_node_id;
+typedef unsigned archive_node_id;
 
 /** Numerical ID value to refer to a string. */
-typedef unsigned int archive_atom;
+typedef unsigned archive_atom;
 
 
 /** This class stores all properties needed to record/retrieve the state
@@ -88,7 +88,7 @@ public:
        void add_bool(const std::string &name, bool value);
 
        /** Add property of type "unsigned int" to node. */
-       void add_unsigned(const std::string &name, unsigned int value);
+       void add_unsigned(const std::string &name, unsigned value);
 
        /** Add property of type "string" to node. */
        void add_string(const std::string &name, const std::string &value);
@@ -98,23 +98,23 @@ public:
 
        /** Retrieve property of type "bool" from node.
         *  @return "true" if property was found, "false" otherwise */
-       bool find_bool(const std::string &name, bool &ret) const;
+       bool find_bool(const std::string &name, bool &ret, unsigned index = 0) const;
 
        /** Retrieve property of type "unsigned" from node.
         *  @return "true" if property was found, "false" otherwise */
-       bool find_unsigned(const std::string &name, unsigned int &ret) const;
+       bool find_unsigned(const std::string &name, unsigned &ret, unsigned index = 0) const;
 
        /** Retrieve property of type "string" from node.
         *  @return "true" if property was found, "false" otherwise */
-       bool find_string(const std::string &name, std::string &ret) const;
+       bool find_string(const std::string &name, std::string &ret, unsigned index = 0) const;
 
        /** Retrieve property of type "ex" from node.
         *  @return "true" if property was found, "false" otherwise */
-       bool find_ex(const std::string &name, ex &ret, const lst &sym_lst, unsigned int index = 0) const;
+       bool find_ex(const std::string &name, ex &ret, const lst &sym_lst, unsigned index = 0) const;
 
        /** Retrieve property of type "ex" from node, returning the node of
         *  the sub-expression. */
-       const archive_node &find_ex_node(const std::string &name, unsigned int index = 0) const;
+       const archive_node &find_ex_node(const std::string &name, unsigned index = 0) const;
 
        /** Return vector of properties stored in node. */
        void get_properties(propinfovector &v) const;
@@ -131,7 +131,7 @@ private:
        /** Archived property (data type, name and associated data) */
        struct property {
                property() {}
-               property(archive_atom n, property_type t, unsigned int v) : type(t), name(n), value(v) {}
+               property(archive_atom n, property_type t, unsigned v) : type(t), name(n), value(v) {}
                ~property() {}
 
                property(const property &other) : type(other.type), name(other.name), value(other.value) {}
@@ -139,7 +139,7 @@ private:
 
                property_type type; /**< Data type of property. */
                archive_atom name;  /**< Name of property. */
-               unsigned int value; /**< Stored value. */
+               unsigned value; /**< Stored value. */
        };
 
        /** Reference to the archive to which this node belongs. */
@@ -190,19 +190,19 @@ public:
        /** Retrieve expression from archive by index.
         *  @param sym_lst list of pre-defined symbols
      *  @see count_expressions */
-       ex unarchive_ex(const lst &sym_lst, unsigned int index = 0) const;
+       ex unarchive_ex(const lst &sym_lst, unsigned index = 0) const;
 
        /** Retrieve expression and its name from archive by index.
         *  @param sym_lst list of pre-defined symbols
         *  @param name receives the name of the expression
      *  @see count_expressions */
-       ex unarchive_ex(const lst &sym_lst, std::string &name, unsigned int index = 0) const;
+       ex unarchive_ex(const lst &sym_lst, std::string &name, unsigned index = 0) const;
 
        /** Return number of archived expressions. */
-       unsigned int num_expressions(void) const;
+       unsigned num_expressions(void) const;
 
        /** Return reference to top node of an expression specified by index. */
-       const archive_node &get_top_node(unsigned int index = 0) const;
+       const archive_node &get_top_node(unsigned index = 0) const;
 
        /** Clear all archived expressions. */
        void clear(void);
index a2e10f7..eb2c85d 100644 (file)
@@ -232,7 +232,7 @@ bool basic::has(const ex & other) const
 
 /** Construct new expression by applying the specified function to all
  *  sub-expressions (one level only, not recursively). */
-ex basic::map(map_func f) const
+ex basic::map(map_function & f) const
 {
        unsigned num = nops();
        if (num == 0)
@@ -368,8 +368,12 @@ ex basic::evalm(void) const
 {
        if (nops() == 0)
                return *this;
-       else
-               return map(GiNaC::evalm);
+       else {
+               struct evalm_map_function : public map_function {
+                       ex operator()(const ex & e) { return GiNaC::evalm(e); }
+               } fcn;
+               return map(fcn);
+       }
 }
 
 /** Perform automatic symbolic evaluations on indexed expression that
index 9199e34..4ace75f 100644 (file)
@@ -47,7 +47,14 @@ class print_context;
 
 typedef std::vector<ex> exvector;
 
-typedef ex (*map_func)(const ex & e);
+
+/** Function object for map(). */
+struct map_function {
+       typedef const ex & argument_type;
+       typedef ex result_type;
+       virtual ex operator()(const ex & e) = 0;
+};
+
 
 /** This class is the ABC (abstract base class) of GiNaC's class hierarchy.
  *  It is responsible for the reference counting. */
@@ -106,7 +113,7 @@ public: // only const functions please (may break reference counting)
        virtual ex operator[](const ex & index) const;
        virtual ex operator[](int i) const;
        virtual bool has(const ex & other) const;
-       virtual ex map(map_func f) const;
+       virtual ex map(map_function & f) const;
        virtual int degree(const ex & s) const;
        virtual int ldegree(const ex & s) const;
        virtual ex coeff(const ex & s, int n = 1) const;
index e62def5..c985bd8 100644 (file)
@@ -26,6 +26,7 @@
 #include "ncmul.h"
 #include "symbol.h"
 #include "numeric.h" // for I
+#include "symmetry.h"
 #include "lst.h"
 #include "relational.h"
 #include "print.h"
@@ -86,13 +87,13 @@ clifford::clifford(const ex & b, const ex & mu, unsigned char rl) : inherited(b,
        tinfo_key = TINFO_clifford;
 }
 
-clifford::clifford(unsigned char rl, const exvector & v, bool discardable) : inherited(indexed::unknown, v, discardable), representation_label(rl)
+clifford::clifford(unsigned char rl, const exvector & v, bool discardable) : inherited(sy_none(), v, discardable), representation_label(rl)
 {
        debugmsg("clifford constructor from unsigned char,exvector", LOGLEVEL_CONSTRUCT);
        tinfo_key = TINFO_clifford;
 }
 
-clifford::clifford(unsigned char rl, exvector * vp) : inherited(indexed::unknown, vp), representation_label(rl)
+clifford::clifford(unsigned char rl, exvector * vp) : inherited(sy_none(), vp), representation_label(rl)
 {
        debugmsg("clifford constructor from unsigned char,exvector *", LOGLEVEL_CONSTRUCT);
        tinfo_key = TINFO_clifford;
@@ -391,14 +392,6 @@ ex dirac_trace(const ex & e, unsigned char rl, const ex & trONE)
                else
                        return _ex0();
 
-       } else if (is_ex_exactly_of_type(e, add)) {
-
-               // Trace of sum = sum of traces
-               ex sum = _ex0();
-               for (unsigned i=0; i<e.nops(); i++)
-                       sum += dirac_trace(e.op(i), rl, trONE);
-               return sum;
-
        } else if (is_ex_exactly_of_type(e, mul)) {
 
                // Trace of product: pull out non-clifford factors
@@ -491,9 +484,15 @@ ex dirac_trace(const ex & e, unsigned char rl, const ex & trONE)
 
                        return trONE * trace_string(iv.begin(), num);
                }
-       }
 
-       return _ex0();
+       } else if (e.nops() > 0) {
+
+               // Trace maps to all other container classes (this includes sums)
+               pointer_to_map_function_2args<unsigned char, const ex &> fcn(dirac_trace, rl, trONE);
+               return e.map(fcn);
+
+       } else
+               return _ex0();
 }
 
 ex canonicalize_clifford(const ex & e)
index 487c12d..2877b88 100644 (file)
@@ -151,7 +151,7 @@ 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. */
+ *  to check two expressions for equality. */
 ex canonicalize_clifford(const ex & e);
 
 } // namespace GiNaC
index adcd874..a3876f3 100644 (file)
@@ -26,6 +26,7 @@
 #include "color.h"
 #include "idx.h"
 #include "ncmul.h"
+#include "symmetry.h"
 #include "numeric.h"
 #include "power.h" // for sqrt()
 #include "symbol.h"
@@ -86,13 +87,13 @@ color::color(const ex & b, const ex & i1, unsigned char rl) : inherited(b, i1),
        tinfo_key = TINFO_color;
 }
 
-color::color(unsigned char rl, const exvector & v, bool discardable) : inherited(indexed::unknown, v, discardable), representation_label(rl)
+color::color(unsigned char rl, const exvector & v, bool discardable) : inherited(sy_none(), v, discardable), representation_label(rl)
 {
        debugmsg("color constructor from unsigned char,exvector", LOGLEVEL_CONSTRUCT);
        tinfo_key = TINFO_color;
 }
 
-color::color(unsigned char rl, exvector * vp) : inherited(indexed::unknown, vp), representation_label(rl)
+color::color(unsigned char rl, exvector * vp) : inherited(sy_none(), vp), representation_label(rl)
 {
        debugmsg("color constructor from unsigned char,exvector *", LOGLEVEL_CONSTRUCT);
        tinfo_key = TINFO_color;
@@ -482,7 +483,7 @@ ex color_f(const ex & a, const ex & b, const ex & c)
        if (!ex_to_idx(a).get_dim().is_equal(8) || !ex_to_idx(b).get_dim().is_equal(8) || !ex_to_idx(c).get_dim().is_equal(8))
                throw(std::invalid_argument("index dimension for color_f must be 8"));
 
-       return indexed(su3f(), indexed::antisymmetric, a, b, c);
+       return indexed(su3f(), sy_anti(), a, b, c);
 }
 
 ex color_d(const ex & a, const ex & b, const ex & c)
@@ -492,7 +493,7 @@ ex color_d(const ex & a, const ex & b, const ex & c)
        if (!ex_to_idx(a).get_dim().is_equal(8) || !ex_to_idx(b).get_dim().is_equal(8) || !ex_to_idx(c).get_dim().is_equal(8))
                throw(std::invalid_argument("index dimension for color_d must be 8"));
 
-       return indexed(su3d(), indexed::symmetric, a, b, c);
+       return indexed(su3d(), sy_symm(), a, b, c);
 }
 
 ex color_h(const ex & a, const ex & b, const ex & c)
@@ -517,14 +518,6 @@ ex color_trace(const ex & e, unsigned char rl)
                else
                        return _ex0();
 
-       } else if (is_ex_exactly_of_type(e, add)) {
-
-               // Trace of sum = sum of traces
-               ex sum = _ex0();
-               for (unsigned i=0; i<e.nops(); i++)
-                       sum += color_trace(e.op(i), rl);
-               return sum;
-
        } else if (is_ex_exactly_of_type(e, mul)) {
 
                // Trace of product: pull out non-color factors
@@ -581,9 +574,15 @@ ex color_trace(const ex & e, unsigned char rl)
                        return delta_tensor(next_to_last_index, last_index) * color_trace(ncmul(v1), rl) / 6
                               + color_h(next_to_last_index, last_index, summation_index) * color_trace(ncmul(v2), rl) / 2;
                }
-       }
 
-       return _ex0();
+       } else if (e.nops() > 0) {
+
+               // Trace maps to all other container classes (this includes sums)
+               pointer_to_map_function_1arg<unsigned char> fcn(color_trace, rl);
+               return e.map(fcn);
+
+       } else
+               return _ex0();
 }
 
 } // namespace GiNaC
index 2754615..5b8b59e 100755 (executable)
@@ -210,7 +210,7 @@ public:
        bool info(unsigned inf) const;
        unsigned nops() const;
        ex & let_op(int i);
-       ex map(map_func f) const;
+       ex map(map_function & f) const;
        ex expand(unsigned options=0) const;
        ex eval(int level=0) const;
        ex evalf(int level=0) const;
@@ -432,7 +432,7 @@ unsigned ${CONTAINER}::nops() const
 
 ${LET_OP_IMPLEMENTATION}
 
-ex ${CONTAINER}::map(map_func f) const
+ex ${CONTAINER}::map(map_function & f) const
 {
        ${STLT} s;
        RESERVE(s,seq.size());
index 26183ef..b1b639b 100644 (file)
@@ -26,6 +26,8 @@
 #include "basic.h"
 #include "operators.h"
 
+#include <functional>
+
 namespace GiNaC {
 
 // Sorry, this is the only constant to pollute the global scope, the other ones
@@ -37,6 +39,7 @@ class symbol;
 class lst;
 class scalar_products;
 
+
 /** Lightweight wrapper for GiNaC's symbolic objects.  Basically all it does is
  *  to hold a pointer to the other objects, manage the reference counting and
  *  provide methods for manipulation of these objects.  (Some people call such
@@ -84,7 +87,8 @@ public:
        unsigned nops() const { return bp->nops(); }
        ex expand(unsigned options=0) const;
        bool has(const ex & other) const { return bp->has(other); }
-       ex map(map_func f) const { return bp->map(f); }
+       ex map(map_function & f) const { return bp->map(f); }
+       ex map(ex (*f)(const ex & e)) const;
        int degree(const ex & s) const { return bp->degree(s); }
        int ldegree(const ex & s) const { return bp->ldegree(s); }
        ex coeff(const ex & s, int n = 1) const { return bp->coeff(s, n); }
@@ -435,6 +439,70 @@ inline bool is_zero(const ex & thisex)
 inline void swap(ex & e1, ex & e2)
 { e1.swap(e2); }
 
+
+/* Function objects for STL sort() etc. */
+struct ex_is_less : public std::binary_function<ex, ex, bool> {
+       bool operator() (const ex &lh, const ex &rh) const { return lh.compare(rh) < 0; }
+};
+
+struct ex_is_equal : public std::binary_function<ex, ex, bool> {
+       bool operator() (const ex &lh, const ex &rh) const { return lh.is_equal(rh); }
+};
+
+struct ex_swap : public std::binary_function<ex, ex, void> {
+       void operator() (ex &lh, ex &rh) const { lh.swap(rh); }
+};
+
+
+/* Convert function pointer to function object suitable for map(). */
+class pointer_to_map_function : public map_function {
+protected:
+       ex (*ptr)(const ex &);
+public:
+       explicit pointer_to_map_function(ex (*x)(const ex &)) : ptr(x) {}
+       ex operator()(const ex & e) { return ptr(e); }
+};
+
+template<class T1>
+class pointer_to_map_function_1arg : public map_function {
+protected:
+       ex (*ptr)(const ex &, T1);
+       T1 arg1;
+public:
+       explicit pointer_to_map_function_1arg(ex (*x)(const ex &, T1), T1 a1) : ptr(x), arg1(a1) {}
+       ex operator()(const ex & e) { return ptr(e, arg1); }
+};
+
+template<class T1, class T2>
+class pointer_to_map_function_2args : public map_function {
+protected:
+       ex (*ptr)(const ex &, T1, T2);
+       T1 arg1;
+       T2 arg2;
+public:
+       explicit pointer_to_map_function_2args(ex (*x)(const ex &, T1, T2), T1 a1, T2 a2) : ptr(x), arg1(a1), arg2(a2) {}
+       ex operator()(const ex & e) { return ptr(e, arg1, arg2); }
+};
+
+template<class T1, class T2, class T3>
+class pointer_to_map_function_3args : public map_function {
+protected:
+       ex (*ptr)(const ex &, T1, T2, T3);
+       T1 arg1;
+       T2 arg2;
+       T3 arg3;
+public:
+       explicit pointer_to_map_function_3args(ex (*x)(const ex &, T1, T2, T3), T1 a1, T2 a2, T3 a3) : ptr(x), arg1(a1), arg2(a2), arg3(a3) {}
+       ex operator()(const ex & e) { return ptr(e, arg1, arg2, arg3); }
+};
+
+inline ex ex::map(ex (*f)(const ex & e)) const
+{
+       pointer_to_map_function fcn(f);
+       return bp->map(fcn);
+}
+
+
 } // namespace GiNaC
 
 #endif // ndef __GINAC_EX_H__
index 22669e5..6370fc2 100644 (file)
@@ -299,7 +299,7 @@ ex &expairseq::let_op(int i)
        throw(std::logic_error("let_op not defined for expairseq and derived classes (add,mul,...)"));
 }
 
-ex expairseq::map(map_func f) const
+ex expairseq::map(map_function & f) const
 {
        epvector *v = new epvector;
        v->reserve(seq.size());
index be5e014..0bfd2ed 100644 (file)
@@ -93,7 +93,7 @@ public:
        unsigned nops() const;
        ex op(int i) const;
        ex & let_op(int i);
-       ex map(map_func f) const;
+       ex map(map_function & f) const;
        ex eval(int level=0) const;
        ex evalf(int level=0) const;
        ex normal(lst &sym_lst, lst &repl_lst, int level=0) const;
index e228709..7d45d6d 100755 (executable)
@@ -242,6 +242,7 @@ if (!automatic_typecheck) { \\
 namespace GiNaC {
 
 class function;
+class symmetry;
 
 typedef ex (* eval_funcp)();
 typedef ex (* evalf_funcp)();
@@ -276,6 +277,7 @@ $series_func_interface
        function_options & remember(unsigned size, unsigned assoc_size=0,
                                    unsigned strategy=remember_strategies::delete_never);
        function_options & overloaded(unsigned o);
+       function_options & set_symmetry(const symmetry & s);
        void test_and_set_nparams(unsigned n);
        std::string get_name(void) const { return name; }
        unsigned get_nparams(void) const { return nparams; }
@@ -303,6 +305,8 @@ protected:
        unsigned remember_strategy;
 
        unsigned functions_with_same_name;
+
+       ex symtree;
 };
 
 /** The class function is used to implement builtin functions like sin, cos...
@@ -423,6 +427,7 @@ $implementation=<<END_OF_IMPLEMENTATION;
 #include "function.h"
 #include "ex.h"
 #include "lst.h"
+#include "symmetry.h"
 #include "print.h"
 #include "archive.h"
 #include "inifcns.h"
@@ -461,6 +466,7 @@ void function_options::initialize(void)
        use_return_type=false;
        use_remember=false;
        functions_with_same_name=1;
+       symtree = 0;
 }
 
 function_options & function_options::set_name(std::string const & n,
@@ -518,6 +524,12 @@ function_options & function_options::overloaded(unsigned o)
        functions_with_same_name=o;
        return *this;
 }
+
+function_options & function_options::set_symmetry(const symmetry & s)
+{
+       symtree = s;
+       return *this;
+}
        
 void function_options::test_and_set_nparams(unsigned n)
 {
@@ -719,17 +731,32 @@ ex function::eval(int level) const
                return function(serial,evalchildren(level));
        }
 
-       if (registered_functions()[serial].eval_f==0) {
+       const function_options &opt = registered_functions()[serial];
+
+       // Canonicalize argument order according to the symmetry properties
+       if (seq.size() > 1 && !(opt.symtree.is_zero())) {
+               exvector v = seq;
+               GINAC_ASSERT(is_ex_exactly_of_type(opt.symtree, symmetry));
+               int sig = canonicalize(v.begin(), ex_to_symmetry(opt.symtree));
+               if (sig != INT_MAX) {
+                       // Something has changed while sorting arguments, more evaluations later
+                       if (sig == 0)
+                               return _ex0();
+                       return ex(sig) * thisexprseq(v);
+               }
+       }
+
+       if (opt.eval_f==0) {
                return this->hold();
        }
 
-       bool use_remember=registered_functions()[serial].use_remember;
+       bool use_remember=opt.use_remember;
        ex eval_result;
        if (use_remember && lookup_remember_table(eval_result)) {
                return eval_result;
        }
 
-       switch (registered_functions()[serial].nparams) {
+       switch (opt.nparams) {
                // the following lines have been generated for max. ${maxargs} parameters
 ${eval_switch_statement}
                // end of generated lines
index e7684e5..3c0e9e8 100644 (file)
 #include "relational.h"
 #include "structure.h"
 #include "symbol.h"
+#include "pseries.h"
 #include "wildcard.h"
+#include "symmetry.h"
 
 #include "expair.h"
 #include "expairseq.h"
 #include "add.h"
 #include "mul.h"
-#include "pseries.h"
 
 #include "exprseq.h"
 #include "function.h"
index c72e478..e2b76be 100644 (file)
@@ -491,7 +491,7 @@ void find_free_and_dummy(exvector::const_iterator it, exvector::const_iterator i
        // Sort index vector. This will cause dummy indices come to lie next
        // to each other (because the sort order is defined to guarantee this).
        exvector v(it, itend);
-       shaker_sort(v.begin(), v.end(), ex_is_less());
+       shaker_sort(v.begin(), v.end(), ex_is_less(), ex_swap());
 
        // Find dummy pairs and free indices
        it = v.begin(); itend = v.end();
index 3f2cb71..fd96eeb 100644 (file)
@@ -29,8 +29,8 @@
 #include "mul.h"
 #include "ncmul.h"
 #include "power.h"
+#include "symmetry.h"
 #include "lst.h"
-#include "inifcns.h" // for symmetrize()
 #include "print.h"
 #include "archive.h"
 #include "utils.h"
@@ -44,7 +44,7 @@ GINAC_IMPLEMENT_REGISTERED_CLASS(indexed, exprseq)
 // default constructor, destructor, copy constructor assignment operator and helpers
 //////////
 
-indexed::indexed() : symmetry(unknown)
+indexed::indexed() : symtree(sy_none())
 {
        debugmsg("indexed default constructor", LOGLEVEL_CONSTRUCT);
        tinfo_key = TINFO_indexed;
@@ -53,7 +53,7 @@ indexed::indexed() : symmetry(unknown)
 void indexed::copy(const indexed & other)
 {
        inherited::copy(other);
-       symmetry = other.symmetry;
+       symtree = other.symtree;
 }
 
 DEFAULT_DESTROY(indexed)
@@ -62,97 +62,94 @@ DEFAULT_DESTROY(indexed)
 // other constructors
 //////////
 
-indexed::indexed(const ex & b) : inherited(b), symmetry(unknown)
+indexed::indexed(const ex & b) : inherited(b), symtree(sy_none())
 {
        debugmsg("indexed constructor from ex", LOGLEVEL_CONSTRUCT);
        tinfo_key = TINFO_indexed;
-       assert_all_indices_of_type_idx();
+       validate();
 }
 
-indexed::indexed(const ex & b, const ex & i1) : inherited(b, i1), symmetry(unknown)
+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;
-       assert_all_indices_of_type_idx();
+       validate();
 }
 
-indexed::indexed(const ex & b, const ex & i1, const ex & i2) : inherited(b, i1, i2), symmetry(unknown)
+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;
-       assert_all_indices_of_type_idx();
+       validate();
 }
 
-indexed::indexed(const ex & b, const ex & i1, const ex & i2, const ex & i3) : inherited(b, i1, i2, i3), symmetry(unknown)
+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;
-       assert_all_indices_of_type_idx();
+       validate();
 }
 
-indexed::indexed(const ex & b, const ex & i1, const ex & i2, const ex & i3, const ex & i4) : inherited(b, i1, i2, i3, i4), symmetry(unknown)
+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;
-       assert_all_indices_of_type_idx();
+       validate();
 }
 
-indexed::indexed(const ex & b, symmetry_type symm, const ex & i1, const ex & i2) : inherited(b, i1, i2), symmetry(symm)
+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;
-       assert_all_indices_of_type_idx();
+       validate();
 }
 
-indexed::indexed(const ex & b, symmetry_type symm, const ex & i1, const ex & i2, const ex & i3) : inherited(b, i1, i2, i3), symmetry(symm)
+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;
-       assert_all_indices_of_type_idx();
+       validate();
 }
 
-indexed::indexed(const ex & b, symmetry_type symm, const ex & i1, const ex & i2, const ex & i3, const ex & i4) : inherited(b, i1, i2, i3, i4), symmetry(symm)
+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;
-       assert_all_indices_of_type_idx();
+       validate();
 }
 
-indexed::indexed(const ex & b, const exvector & v) : inherited(b), symmetry(unknown)
+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;
-       assert_all_indices_of_type_idx();
+       validate();
 }
 
-indexed::indexed(const ex & b, symmetry_type symm, const exvector & v) : inherited(b), symmetry(symm)
+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;
-       assert_all_indices_of_type_idx();
+       validate();
 }
 
-indexed::indexed(symmetry_type symm, const exprseq & es) : inherited(es), symmetry(symm)
+indexed::indexed(const symmetry & symm, const exprseq & es) : inherited(es), symtree(symm)
 {
        debugmsg("indexed constructor from symmetry,exprseq", LOGLEVEL_CONSTRUCT);
        tinfo_key = TINFO_indexed;
-       assert_all_indices_of_type_idx();
 }
 
-indexed::indexed(symmetry_type symm, const exvector & v, bool discardable) : inherited(v, discardable), symmetry(symm)
+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;
-       assert_all_indices_of_type_idx();
 }
 
-indexed::indexed(symmetry_type symm, exvector * vp) : inherited(vp), symmetry(symm)
+indexed::indexed(const symmetry & symm, exvector * vp) : inherited(vp), symtree(symm)
 {
        debugmsg("indexed constructor from symmetry,exvector *", LOGLEVEL_CONSTRUCT);
        tinfo_key = TINFO_indexed;
-       assert_all_indices_of_type_idx();
 }
 
 //////////
@@ -162,15 +159,29 @@ indexed::indexed(symmetry_type symm, exvector * vp) : inherited(vp), symmetry(sy
 indexed::indexed(const archive_node &n, const lst &sym_lst) : inherited(n, sym_lst)
 {
        debugmsg("indexed constructor from archive_node", LOGLEVEL_CONSTRUCT);
-       unsigned int symm;
-       if (!(n.find_unsigned("symmetry", symm)))
-               throw (std::runtime_error("unknown indexed symmetry type in archive"));
+       if (!n.find_ex("symmetry", symtree, sym_lst)) {
+               // GiNaC versions <= 0.9.0 had an unsigned "symmetry" property
+               unsigned symm = 0;
+               n.find_unsigned("symmetry", symm);
+               switch (symm) {
+                       case 1:
+                               symtree = sy_symm();
+                               break;
+                       case 2:
+                               symtree = sy_anti();
+                               break;
+                       default:
+                               symtree = sy_none();
+                               break;
+               }
+               ex_to_nonconst_symmetry(symtree).validate(seq.size() - 1);
+       }
 }
 
 void indexed::archive(archive_node &n) const
 {
        inherited::archive(n);
-       n.add_unsigned("symmetry", symmetry);
+       n.add_ex("symmetry", symtree);
 }
 
 DEFAULT_UNARCHIVE(indexed)
@@ -188,12 +199,8 @@ void indexed::print(const print_context & c, unsigned level) const
 
                c.s << std::string(level, ' ') << class_name()
                    << std::hex << ", hash=0x" << hashvalue << ", flags=0x" << flags << std::dec
-                   << ", " << seq.size()-1 << " indices";
-               switch (symmetry) {
-                       case symmetric: c.s << ", symmetric"; break;
-                       case antisymmetric: c.s << ", antisymmetric"; break;
-                       default: break;
-               }
+                   << ", " << seq.size()-1 << " indices"
+                   << ", symmetry=" << symtree << std::endl;
                c.s << std::endl;
                unsigned delta_indent = static_cast<const print_tree &>(c).delta_indent;
                seq[0].print(c, level + delta_indent);
@@ -248,57 +255,11 @@ int indexed::compare_same_type(const basic & other) const
        return inherited::compare_same_type(other);
 }
 
-// The main difference between sort_index_vector() and canonicalize_indices()
-// is that the latter takes the symmetry of the object into account. Once we
-// implement mixed symmetries, canonicalize_indices() will only be able to
-// 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. This operation only makes
- *  sense if the object carrying these indices is either symmetric or totally
- *  antisymmetric with respect to the indices.
- *
- *  @param itbegin Start of index vector
- *  @param itend End of index vector
- *  @param antisymm Whether the object is antisymmetric
- *  @return the sign introduced by the reordering of the indices if the object
- *          is antisymmetric (or 0 if two equal indices are encountered). For
- *          symmetric objects, this is always +1. If the index vector was
- *          already in a canonic order this function returns INT_MAX. */
-static int canonicalize_indices(exvector::iterator itbegin, exvector::iterator itend, bool antisymm)
-{
-       bool something_changed = false;
-       int sig = 1;
-
-       // Simple bubble sort algorithm should be sufficient for the small
-       // number of indices expected
-       exvector::iterator it1 = itbegin, next_to_last_idx = itend - 1;
-       while (it1 != next_to_last_idx) {
-               exvector::iterator it2 = it1 + 1;
-               while (it2 != itend) {
-                       int cmpval = it1->compare(*it2);
-                       if (cmpval == 1) {
-                               it1->swap(*it2);
-                               something_changed = true;
-                               if (antisymm)
-                                       sig = -sig;
-                       } else if (cmpval == 0 && antisymm) {
-                               something_changed = true;
-                               sig = 0;
-                       }
-                       it2++;
-               }
-               it1++;
-       }
-
-       return something_changed ? sig : INT_MAX;
-}
-
 ex indexed::eval(int level) const
 {
        // First evaluate children, then we will end up here again
        if (level > 1)
-               return indexed(symmetry, evalchildren(level));
+               return indexed(ex_to_symmetry(symtree), evalchildren(level));
 
        const ex &base = seq[0];
 
@@ -315,9 +276,10 @@ ex indexed::eval(int level) const
        }
 
        // Canonicalize indices according to the symmetry properties
-       if (seq.size() > 2 && (symmetry == symmetric || symmetry == antisymmetric)) {
-               exvector v(seq);
-               int sig = canonicalize_indices(v.begin() + 1, v.end(), symmetry == antisymmetric);
+       if (seq.size() > 2) {
+               exvector v = seq;
+               GINAC_ASSERT(is_ex_exactly_of_type(symtree, symmetry));
+               int sig = canonicalize(v.begin() + 1, ex_to_symmetry(symtree));
                if (sig != INT_MAX) {
                        // Something has changed while sorting indices, more evaluations later
                        if (sig == 0)
@@ -350,12 +312,12 @@ ex indexed::coeff(const ex & s, int n) const
 
 ex indexed::thisexprseq(const exvector & v) const
 {
-       return indexed(symmetry, v);
+       return indexed(ex_to_symmetry(symtree), v);
 }
 
 ex indexed::thisexprseq(exvector * vp) const
 {
-       return indexed(symmetry, vp);
+       return indexed(ex_to_symmetry(symtree), vp);
 }
 
 ex indexed::expand(unsigned options) const
@@ -429,10 +391,10 @@ void indexed::printindices(const print_context & c, unsigned level) const
        }
 }
 
-/** Check whether all indices are of class idx. This function is used
- *  internally to make sure that all constructed indexed objects really
- *  carry indices and not some other classes. */
-void indexed::assert_all_indices_of_type_idx(void) const
+/** Check whether all indices are of class idx and validate the symmetry
+ *  tree. This function is used internally to make sure that all constructed
+ *  indexed objects really carry indices and not some other classes. */
+void indexed::validate(void) const
 {
        GINAC_ASSERT(seq.size() > 0);
        exvector::const_iterator it = seq.begin() + 1, itend = seq.end();
@@ -441,6 +403,12 @@ void indexed::assert_all_indices_of_type_idx(void) const
                        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))
+                       throw(std::invalid_argument("symmetry of indexed object must be of type symmetry"));
+               ex_to_nonconst_symmetry(symtree).validate(seq.size() - 1);
+       }
 }
 
 //////////
@@ -682,13 +650,26 @@ try_again:
                        }
 
                        // Contraction of symmetric with antisymmetric object is zero
-                       if ((ex_to_indexed(*it1).symmetry == indexed::symmetric &&
-                            ex_to_indexed(*it2).symmetry == indexed::antisymmetric
-                         || ex_to_indexed(*it1).symmetry == indexed::antisymmetric &&
-                            ex_to_indexed(*it2).symmetry == indexed::symmetric)
-                        && dummy.size() > 1) {
-                               free_indices.clear();
-                               return _ex0();
+                       if (dummy.size() > 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<dummy.size()-1; idx1++) {
+                                       for (unsigned idx2=idx1+1; idx2<dummy.size(); 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
index babb9c2..2221070 100644 (file)
@@ -31,6 +31,7 @@ namespace GiNaC {
 
 
 class scalar_products;
+class symmetry;
 
 /** This class holds an indexed expression. It consists of a 'base' expression
  *  (the expression being indexed) which can be accessed as op(0), and n (n >= 0)
@@ -42,16 +43,6 @@ class indexed : public exprseq
        friend ex simplify_indexed(const ex & e, exvector & free_indices, exvector & dummy_indices, const scalar_products & sp);
        friend ex simplify_indexed_product(const ex & e, exvector & free_indices, exvector & dummy_indices, const scalar_products & sp);
 
-       // types
-public:
-       /** Type of symmetry of the object with respect to commutation of its indices. */
-       typedef enum {
-               unknown,       /**< symmetry properties unknown */
-               symmetric,     /**< totally symmetric */
-               antisymmetric, /**< totally antisymmetric */
-               mixed          /**< mixed symmetry (unimplemented) */
-       } symmetry_type;
-
        // other constructors
 public:
        /** Construct indexed object with no index.
@@ -102,7 +93,7 @@ public:
         *  @param i1 First index
         *  @param i2 Second index
         *  @return newly constructed indexed object */
-       indexed(const ex & b, symmetry_type symm, const ex & i1, const ex & i2);
+       indexed(const ex & b, const symmetry & symm, const ex & i1, const ex & i2);
 
        /** Construct indexed object with three indices and a specified symmetry.
         *  The indices must be of class idx.
@@ -113,7 +104,7 @@ public:
         *  @param i2 Second index
         *  @param i3 Third index
         *  @return newly constructed indexed object */
-       indexed(const ex & b, symmetry_type symm, const ex & i1, const ex & i2, const ex & i3);
+       indexed(const ex & b, const symmetry & symm, const ex & i1, const ex & i2, const ex & i3);
 
        /** Construct indexed object with four indices and a specified symmetry. The
         *  indices must be of class idx.
@@ -125,7 +116,7 @@ public:
         *  @param i3 Third index
         *  @param i4 Fourth index
         *  @return newly constructed indexed object */
-       indexed(const ex & b, symmetry_type symm, const ex & i1, const ex & i2, const ex & i3, const ex & i4);
+       indexed(const ex & b, const symmetry & symm, const ex & i1, const ex & i2, const ex & i3, const ex & i4);
 
        /** Construct indexed object with a specified vector of indices. The indices
         *  must be of class idx.
@@ -142,12 +133,12 @@ public:
         *  @param symm Symmetry of indices
         *  @param iv Vector of indices
         *  @return newly constructed indexed object */
-       indexed(const ex & b, symmetry_type symm, const exvector & iv);
+       indexed(const ex & b, const symmetry & symm, const exvector & iv);
 
        // internal constructors
-       indexed(symmetry_type symm, const exprseq & es);
-       indexed(symmetry_type symm, const exvector & v, bool discardable = false);
-       indexed(symmetry_type symm, exvector * vp); // vp will be deleted
+       indexed(const symmetry & symm, const exprseq & es);
+       indexed(const symmetry & symm, const exvector & v, bool discardable = false);
+       indexed(const symmetry & symm, exvector * vp); // vp will be deleted
 
        // functions overriding virtual functions from base classes
 public:
@@ -188,13 +179,16 @@ public:
         *  with a given index. */
        bool has_dummy_index_for(const ex & i) const;
 
+       /** Return symmetry properties. */
+       ex get_symmetry(void) const {return symtree;}
+
 protected:
        void printindices(const print_context & c, unsigned level) const;
-       void assert_all_indices_of_type_idx(void) const;
+       void validate(void) const;
 
        // member variables
 protected:
-       symmetry_type symmetry; /**< Index symmetry */
+       ex symtree; /**< Index symmetry (tree of symmetry objects) */
 };
 
 
index 0ee9b27..e2714cf 100644 (file)
@@ -29,8 +29,6 @@
 #include "lst.h"
 #include "matrix.h"
 #include "mul.h"
-#include "ncmul.h"
-#include "numeric.h"
 #include "power.h"
 #include "relational.h"
 #include "pseries.h"
@@ -509,111 +507,6 @@ ex lsolve(const ex &eqns, const ex &symbols)
        return sollist;
 }
 
-// Symmetrize/antisymmetrize over a vector of objects
-static ex symm(const ex & e, exvector::const_iterator first, exvector::const_iterator last, bool asymmetric)
-{
-       // Need at least 2 objects for this operation
-       int num = last - first;
-       if (num < 2)
-               return e;
-
-       // Transform object vector to a list
-       exlist iv_lst;
-       iv_lst.insert(iv_lst.begin(), first, last);
-       lst orig_lst(iv_lst, true);
-
-       // Create index vectors for permutation
-       unsigned *iv = new unsigned[num], *iv2;
-       for (unsigned i=0; i<num; i++)
-               iv[i] = i;
-       iv2 = (asymmetric ? new unsigned[num] : NULL);
-
-       // Loop over all permutations (the first permutation, which is the
-       // identity, is unrolled)
-       ex sum = e;
-       while (std::next_permutation(iv, iv + num)) {
-               lst new_lst;
-               for (unsigned i=0; i<num; i++)
-                       new_lst.append(orig_lst.op(iv[i]));
-               ex term = e.subs(orig_lst, new_lst);
-               if (asymmetric) {
-                       memcpy(iv2, iv, num * sizeof(unsigned));
-                       term *= permutation_sign(iv2, iv2 + num);
-               }
-               sum += term;
-       }
-
-       delete[] iv;
-       delete[] iv2;
-
-       return sum / factorial(numeric(num));
-}
-
-ex symmetrize(const ex & e, exvector::const_iterator first, exvector::const_iterator last)
-{
-       return symm(e, first, last, false);
-}
-
-ex antisymmetrize(const ex & e, exvector::const_iterator first, exvector::const_iterator last)
-{
-       return symm(e, first, last, true);
-}
-
-ex symmetrize_cyclic(const ex & e, exvector::const_iterator first, exvector::const_iterator last)
-{
-       // Need at least 2 objects for this operation
-       int num = last - first;
-       if (num < 2)
-               return e;
-
-       // Transform object vector to a list
-       exlist iv_lst;
-       iv_lst.insert(iv_lst.begin(), first, last);
-       lst orig_lst(iv_lst, true);
-       lst new_lst = orig_lst;
-
-       // Loop over all cyclic permutations (the first permutation, which is
-       // the identity, is unrolled)
-       ex sum = e;
-       for (unsigned i=0; i<num-1; i++) {
-               ex perm = new_lst.op(0);
-               new_lst.remove_first().append(perm);
-               sum += e.subs(orig_lst, new_lst);
-       }
-       return sum / num;
-}
-
-/** Symmetrize expression over a list of objects (symbols, indices). */
-ex ex::symmetrize(const lst & l) const
-{
-       exvector v;
-       v.reserve(l.nops());
-       for (unsigned i=0; i<l.nops(); i++)
-               v.push_back(l.op(i));
-       return symm(*this, v.begin(), v.end(), false);
-}
-
-/** Antisymmetrize expression over a list of objects (symbols, indices). */
-ex ex::antisymmetrize(const lst & l) const
-{
-       exvector v;
-       v.reserve(l.nops());
-       for (unsigned i=0; i<l.nops(); i++)
-               v.push_back(l.op(i));
-       return symm(*this, v.begin(), v.end(), true);
-}
-
-/** Symmetrize expression by cyclic permutation over a list of objects
- *  (symbols, indices). */
-ex ex::symmetrize_cyclic(const lst & l) const
-{
-       exvector v;
-       v.reserve(l.nops());
-       for (unsigned i=0; i<l.nops(); i++)
-               v.push_back(l.op(i));
-       return GiNaC::symmetrize_cyclic(*this, v.begin(), v.end());
-}
-
 /** Force inclusion of functions from initcns_gamma and inifcns_zeta
  *  for static lib (so ginsh will see them). */
 unsigned force_include_tgamma = function_index_tgamma;
index 314228f..fd0d9e4 100644 (file)
@@ -133,35 +133,6 @@ DECLARE_FUNCTION_2P(Derivative)
 
 ex lsolve(const ex &eqns, const ex &symbols);
 
-/** Symmetrize expression over a set of objects (symbols, indices). */
-ex symmetrize(const ex & e, exvector::const_iterator first, exvector::const_iterator last);
-
-/** Symmetrize expression over a set of objects (symbols, indices). */
-inline ex symmetrize(const ex & e, const exvector & v)
-{
-       return symmetrize(e, v.begin(), v.end());
-}
-
-/** Antisymmetrize expression over a set of objects (symbols, indices). */
-ex antisymmetrize(const ex & e, exvector::const_iterator first, exvector::const_iterator last);
-
-/** Antisymmetrize expression over a set of objects (symbols, indices). */
-inline ex antisymmetrize(const ex & e, const exvector & v)
-{
-       return antisymmetrize(e, v.begin(), v.end());
-}
-
-/** Symmetrize expression by cyclic permuation over a set of objects
- *  (symbols, indices). */
-ex symmetrize_cyclic(const ex & e, exvector::const_iterator first, exvector::const_iterator last);
-
-/** Symmetrize expression by cyclic permutation over a set of objects
- *  (symbols, indices). */
-inline ex symmetrize_cyclic(const ex & e, const exvector & v)
-{
-       return symmetrize(e, v.begin(), v.end());
-}
-
 /** Check whether a function is the Order (O(n)) function. */
 inline bool is_order_function(const ex & e)
 {
index 8fa1998..dfc5ee1 100644 (file)
@@ -31,6 +31,7 @@
 #include "power.h"
 #include "relational.h"
 #include "symbol.h"
+#include "symmetry.h"
 #include "utils.h"
 
 namespace GiNaC {
@@ -319,7 +320,8 @@ REGISTER_FUNCTION(beta, eval_func(beta_eval).
                         evalf_func(beta_evalf).
                         derivative_func(beta_deriv).
                         series_func(beta_series).
-                        latex_name("\\mbox{B}"));
+                        latex_name("\\mbox{B}").
+                                               set_symmetry(sy_symm(0, 1)));
 
 
 //////////
index eefd438..8996625 100644 (file)
@@ -30,7 +30,7 @@
 #include "mul.h"
 #include "ncmul.h"
 #include "numeric.h"
-#include "inifcns.h"
+#include "inifcns.h" // for log() in power::derivative()
 #include "matrix.h"
 #include "symbol.h"
 #include "print.h"
@@ -242,7 +242,7 @@ ex & power::let_op(int i)
        return i==0 ? basis : exponent;
 }
 
-ex power::map(map_func f) const
+ex power::map(map_function & f) const
 {
        return (new power(f(basis), f(exponent)))->setflag(status_flags::dynallocated);
 }
index 0a15654..d13463a 100644 (file)
@@ -53,7 +53,7 @@ public:
        bool info(unsigned inf) const;
        unsigned nops() const;
        ex & let_op(int i);
-       ex map(map_func f) const;
+       ex map(map_function & f) const;
        int degree(const ex & s) const;
        int ldegree(const ex & s) const;
        ex coeff(const ex & s, int n = 1) const;
index 1477449..88f94ce 100644 (file)
@@ -25,7 +25,7 @@
 
 #include "pseries.h"
 #include "add.h"
-#include "inifcns.h"
+#include "inifcns.h" // for Order function
 #include "lst.h"
 #include "mul.h"
 #include "power.h"
index c4b3ec0..d539b70 100644 (file)
@@ -172,11 +172,6 @@ ex & relational::let_op(int i)
        return i==0 ? lh : rh;
 }
 
-ex relational::map(map_func f) const
-{
-       return (new relational(f(lh), f(rh), o))->setflag(status_flags::dynallocated);
-}
-       
 ex relational::eval(int level) const
 {
        if (level==1)
index 8f4efed..c3b223b 100644 (file)
@@ -56,7 +56,6 @@ public:
        bool info(unsigned inf) const;
        unsigned nops() const;
        ex & let_op(int i);
-       ex map(map_func f) const;
        ex eval(int level=0) const;
        ex evalf(int level=0) const;
        ex normal(lst &sym_lst, lst &repl_lst, int level=0) const;
index f4e4945..4a10f54 100644 (file)
@@ -26,6 +26,7 @@
 #include "tensor.h"
 #include "idx.h"
 #include "indexed.h"
+#include "symmetry.h"
 #include "relational.h"
 #include "lst.h"
 #include "numeric.h"
@@ -501,7 +502,7 @@ ex delta_tensor(const ex & i1, const ex & i2)
        if (!is_ex_of_type(i1, idx) || !is_ex_of_type(i2, idx))
                throw(std::invalid_argument("indices of delta tensor must be of type idx"));
 
-       return indexed(tensdelta(), indexed::symmetric, i1, i2);
+       return indexed(tensdelta(), sy_symm(), i1, i2);
 }
 
 ex metric_tensor(const ex & i1, const ex & i2)
@@ -509,7 +510,7 @@ ex metric_tensor(const ex & i1, const ex & i2)
        if (!is_ex_of_type(i1, varidx) || !is_ex_of_type(i2, varidx))
                throw(std::invalid_argument("indices of metric tensor must be of type varidx"));
 
-       return indexed(tensmetric(), indexed::symmetric, i1, i2);
+       return indexed(tensmetric(), sy_symm(), i1, i2);
 }
 
 ex lorentz_g(const ex & i1, const ex & i2, bool pos_sig)
@@ -517,7 +518,7 @@ ex lorentz_g(const ex & i1, const ex & i2, bool pos_sig)
        if (!is_ex_of_type(i1, varidx) || !is_ex_of_type(i2, varidx))
                throw(std::invalid_argument("indices of metric tensor must be of type varidx"));
 
-       return indexed(minkmetric(pos_sig), indexed::symmetric, i1, i2);
+       return indexed(minkmetric(pos_sig), sy_symm(), i1, i2);
 }
 
 ex spinor_metric(const ex & i1, const ex & i2)
@@ -527,7 +528,7 @@ ex spinor_metric(const ex & i1, const ex & i2)
        if (!ex_to_idx(i1).get_dim().is_equal(2) || !ex_to_idx(i2).get_dim().is_equal(2))
                throw(std::runtime_error("index dimension for spinor metric must be 2"));
 
-       return indexed(spinmetric(), indexed::antisymmetric, i1, i2);
+       return indexed(spinmetric(), sy_anti(), i1, i2);
 }
 
 ex epsilon_tensor(const ex & i1, const ex & i2)
@@ -541,7 +542,7 @@ ex epsilon_tensor(const ex & i1, const ex & i2)
        if (!ex_to_idx(i1).get_dim().is_equal(_ex2()))
                throw(std::runtime_error("index dimension of epsilon tensor must match number of indices"));
 
-       return indexed(tensepsilon(), indexed::antisymmetric, i1, i2);
+       return indexed(tensepsilon(), sy_anti(), i1, i2);
 }
 
 ex epsilon_tensor(const ex & i1, const ex & i2, const ex & i3)
@@ -555,7 +556,7 @@ ex epsilon_tensor(const ex & i1, const ex & i2, const ex & i3)
        if (!ex_to_idx(i1).get_dim().is_equal(_ex3()))
                throw(std::runtime_error("index dimension of epsilon tensor must match number of indices"));
 
-       return indexed(tensepsilon(), indexed::antisymmetric, i1, i2, i3);
+       return indexed(tensepsilon(), sy_anti(), i1, i2, i3);
 }
 
 ex lorentz_eps(const ex & i1, const ex & i2, const ex & i3, const ex & i4, bool pos_sig)
@@ -569,7 +570,7 @@ ex lorentz_eps(const ex & i1, const ex & i2, const ex & i3, const ex & i4, bool
        if (!ex_to_idx(i1).get_dim().is_equal(_ex4()))
                throw(std::runtime_error("index dimension of epsilon tensor must match number of indices"));
 
-       return indexed(tensepsilon(true, pos_sig), indexed::antisymmetric, i1, i2, i3, i4);
+       return indexed(tensepsilon(true, pos_sig), sy_anti(), i1, i2, i3, i4);
 }
 
 ex eps0123(const ex & i1, const ex & i2, const ex & i3, const ex & i4, bool pos_sig)
@@ -581,7 +582,7 @@ ex eps0123(const ex & i1, const ex & i2, const ex & i3, const ex & i4, bool pos_
        if (dim.is_equal(4))
                return lorentz_eps(i1, i2, i3, i4, pos_sig);
        else
-               return indexed(tensepsilon(true, pos_sig), indexed::antisymmetric, i1, i2, i3, i4);
+               return indexed(tensepsilon(true, pos_sig), sy_anti(), i1, i2, i3, i4);
 }
 
 } // namespace GiNaC
index 3faea2f..775cfc7 100644 (file)
@@ -84,6 +84,8 @@ const unsigned TINFO_diracgamma5   = 0x000e100eU;
 
 const unsigned TINFO_wildcard      = 0x000f0001U;
 
+const unsigned TINFO_symmetry      = 0x00100001U;
+
 } // namespace GiNaC
 
 #endif // ndef __GINAC_TINFOS_H__
index 29eef7c..91d30f7 100644 (file)
@@ -28,7 +28,6 @@
 
 #include <string>
 #include <stdexcept>
-#include <functional>
 #if defined(HAVE_SSTREAM)
 #include <sstream>
 #elif defined(HAVE_STRSTREAM)
@@ -183,8 +182,8 @@ int permutation_sign(It first, It last)
        return sign;
 }
 
-template <class It, class Cmp>
-int permutation_sign(It first, It last, Cmp comp)
+template <class It, class Cmp, class Swap>
+int permutation_sign(It first, It last, Cmp comp, Swap swapit)
 {
        if (first == last)
                return 0;
@@ -200,7 +199,7 @@ int permutation_sign(It first, It last, Cmp comp)
                bool swapped = false;
                while (i != first) {
                        if (comp(*i, *other)) {
-                               std::iter_swap(other, i);
+                               swapit(*other, *i);
                                flag = other;
                                swapped = true;
                                sign = -sign;
@@ -219,7 +218,7 @@ int permutation_sign(It first, It last, Cmp comp)
                swapped = false;
                while (i != last) {
                        if (comp(*other, *i)) {
-                               std::iter_swap(i, other);
+                               swapit(*i, *other);
                                flag = other;
                                swapped = true;
                                sign = -sign;
@@ -237,8 +236,8 @@ int permutation_sign(It first, It last, Cmp comp)
 }
 
 /* Implementation of shaker sort, only compares adjacent elements. */
-template <class It, class Cmp>
-void shaker_sort(It first, It last, Cmp comp)
+template <class It, class Cmp, class Swap>
+void shaker_sort(It first, It last, Cmp comp, Swap swapit)
 {
        if (first == last)
                return;
@@ -253,7 +252,7 @@ void shaker_sort(It first, It last, Cmp comp)
                bool swapped = false;
                while (i != first) {
                        if (comp(*i, *other)) {
-                               std::iter_swap(other, i);
+                               swapit(*other, *i);
                                flag = other;
                                swapped = true;
                        }
@@ -270,7 +269,7 @@ void shaker_sort(It first, It last, Cmp comp)
                swapped = false;
                while (i != last) {
                        if (comp(*other, *i)) {
-                               std::iter_swap(i, other);
+                               swapit(*i, *other);
                                flag = other;
                                swapped = true;
                        }
@@ -284,8 +283,8 @@ void shaker_sort(It first, It last, Cmp comp)
 }
 
 /* In-place cyclic permutation of a container (no copying, only swapping). */
-template <class It>
-void cyclic_permutation(It first, It last, It new_first)
+template <class It, class Swap>
+void cyclic_permutation(It first, It last, It new_first, Swap swapit)
 {
        unsigned num = last - first;
 again:
@@ -296,7 +295,7 @@ again:
        if (num1 >= num2) {
                It a = first, b = new_first;
                while (b != last) {
-                       std::iter_swap(a, b);
+                       swapit(*a, *b);
                        ++a; ++b;
                }
                if (num1 > num2) {
@@ -308,7 +307,7 @@ again:
                It a = new_first, b = last;
                do {
                        --a; --b;
-                       std::iter_swap(a, b);
+                       swapit(*a, *b);
                } while (a != first);
                last -= num1;
                num = num2;
@@ -316,14 +315,6 @@ again:
        }
 }
 
-/* Function objects for STL sort() etc. */
-struct ex_is_less : public std::binary_function<ex, ex, bool> {
-       bool operator() (const ex &lh, const ex &rh) const { return lh.compare(rh) < 0; }
-};
-
-struct ex_is_equal : public std::binary_function<ex, ex, bool> {
-       bool operator() (const ex &lh, const ex &rh) const { return lh.is_equal(rh); }
-};
 
 // Collection of `construct on first use' wrappers for safely avoiding
 // internal object replication without running into the `static
index 2e76e19..2d79cf7 100644 (file)
@@ -39,7 +39,7 @@ public:
        /** Construct wildcard with specified label. */
        wildcard(unsigned label);
 
-       // functions overriding virtual functions from bases classes
+       // functions overriding virtual functions from base classes
 public:
        void print(const print_context & c, unsigned level = 0) const;
        unsigned calchash(void) const;