3 * Implementation of GiNaC's symmetry definitions. */
6 * GiNaC Copyright (C) 1999-2016 Johannes Gutenberg University Mainz, Germany
8 * This program is free software; you can redistribute it and/or modify
9 * it under the terms of the GNU General Public License as published by
10 * the Free Software Foundation; either version 2 of the License, or
11 * (at your option) any later version.
13 * This program is distributed in the hope that it will be useful,
14 * but WITHOUT ANY WARRANTY; without even the implied warranty of
15 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16 * GNU General Public License for more details.
18 * You should have received a copy of the GNU General Public License
19 * along with this program; if not, write to the Free Software
20 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
26 #include "numeric.h" // for factorial()
27 #include "operators.h"
30 #include "hash_seed.h"
39 GINAC_IMPLEMENT_REGISTERED_CLASS_OPT(symmetry, basic,
40 print_func<print_context>(&symmetry::do_print).
41 print_func<print_tree>(&symmetry::do_print_tree))
44 Some notes about the structure of a symmetry tree:
45 - The leaf nodes of the tree are of type "none", have one index, and no
46 children (of course). They are constructed by the symmetry(unsigned)
48 - Leaf nodes are the only nodes that only have one index.
49 - Container nodes contain two or more children. The "indices" set member
50 is the set union of the index sets of all children, and the "children"
51 vector stores the children themselves.
52 - The index set of each child of a "symm", "anti" or "cycl" node must
53 have the same size. It follows that the children of such a node are
54 either all leaf nodes, or all container nodes with two or more indices.
58 // default constructor
61 symmetry::symmetry() : type(none)
63 setflag(status_flags::evaluated | status_flags::expanded);
70 symmetry::symmetry(unsigned i) : type(none)
73 setflag(status_flags::evaluated | status_flags::expanded);
76 symmetry::symmetry(symmetry_type t, const symmetry &c1, const symmetry &c2) : type(t)
79 setflag(status_flags::evaluated | status_flags::expanded);
86 /** Construct object from archive_node. */
87 void symmetry::read_archive(const archive_node &n, lst &sym_lst)
89 inherited::read_archive(n, sym_lst);
91 if (!(n.find_unsigned("type", t)))
92 throw (std::runtime_error("unknown symmetry type in archive"));
93 type = (symmetry_type)t;
98 if (n.find_ex("child", e, sym_lst, i))
99 add(ex_to<symmetry>(e));
108 if (n.find_unsigned("index", u, i))
116 GINAC_BIND_UNARCHIVER(symmetry);
118 /** Archive the object. */
119 void symmetry::archive(archive_node &n) const
121 inherited::archive(n);
123 n.add_unsigned("type", type);
125 if (children.empty()) {
126 std::set<unsigned>::const_iterator i = indices.begin(), iend = indices.end();
128 n.add_unsigned("index", *i);
132 exvector::const_iterator i = children.begin(), iend = children.end();
134 n.add_ex("child", *i);
141 // functions overriding virtual functions from base classes
144 int symmetry::compare_same_type(const basic & other) const
146 GINAC_ASSERT(is_a<symmetry>(other));
148 // For archiving purposes we need to have an ordering of symmetries.
149 const symmetry &othersymm = ex_to<symmetry>(other);
152 if (type > othersymm.type)
154 if (type < othersymm.type)
157 // Compare the index set.
158 size_t this_size = indices.size();
159 size_t that_size = othersymm.indices.size();
160 if (this_size > that_size)
162 if (this_size < that_size)
164 typedef std::set<unsigned>::const_iterator set_it;
165 set_it end = indices.end();
166 for (set_it i=indices.begin(),j=othersymm.indices.begin(); i!=end; ++i,++j) {
173 // Compare the children.
174 if (children.size() > othersymm.children.size())
176 if (children.size() < othersymm.children.size())
178 for (size_t i=0; i<children.size(); ++i) {
179 int cmpval = ex_to<symmetry>(children[i])
180 .compare_same_type(ex_to<symmetry>(othersymm.children[i]));
188 unsigned symmetry::calchash() const
190 unsigned v = make_hash_seed(typeid(*this));
194 if (!indices.empty())
195 v ^= *(indices.begin());
197 for (exvector::const_iterator i=children.begin(); i!=children.end(); ++i)
204 if (flags & status_flags::evaluated) {
205 setflag(status_flags::hash_calculated);
212 void symmetry::do_print(const print_context & c, unsigned level) const
214 if (children.empty()) {
215 if (indices.size() > 0)
216 c.s << *(indices.begin());
221 case none: c.s << '!'; break;
222 case symmetric: c.s << '+'; break;
223 case antisymmetric: c.s << '-'; break;
224 case cyclic: c.s << '@'; break;
225 default: c.s << '?'; break;
228 size_t num = children.size();
229 for (size_t i=0; i<num; i++) {
230 children[i].print(c);
238 void symmetry::do_print_tree(const print_tree & c, unsigned level) const
240 c.s << std::string(level, ' ') << class_name() << " @" << this
241 << std::hex << ", hash=0x" << hashvalue << ", flags=0x" << flags << std::dec
245 case none: c.s << "none"; break;
246 case symmetric: c.s << "symm"; break;
247 case antisymmetric: c.s << "anti"; break;
248 case cyclic: c.s << "cycl"; break;
249 default: c.s << "<unknown>"; break;
252 c.s << ", indices=(";
253 if (!indices.empty()) {
254 std::set<unsigned>::const_iterator i = indices.begin(), end = indices.end();
262 exvector::const_iterator i = children.begin(), end = children.end();
264 i->print(c, level + c.delta_indent);
270 // non-virtual functions in this class
273 bool symmetry::has_nonsymmetric() const
275 if (type == antisymmetric || type == cyclic)
278 for (exvector::const_iterator i=children.begin(); i!=children.end(); ++i)
279 if (ex_to<symmetry>(*i).has_nonsymmetric())
285 bool symmetry::has_cyclic() const
290 for (exvector::const_iterator i=children.begin(); i!=children.end(); ++i)
291 if (ex_to<symmetry>(*i).has_cyclic())
297 symmetry &symmetry::add(const symmetry &c)
299 // All children must have the same number of indices
300 if (type != none && !children.empty()) {
301 GINAC_ASSERT(is_exactly_a<symmetry>(children[0]));
302 if (ex_to<symmetry>(children[0]).indices.size() != c.indices.size())
303 throw (std::logic_error("symmetry:add(): children must have same number of indices"));
306 // Compute union of indices and check whether the two sets are disjoint
307 std::set<unsigned> un;
308 set_union(indices.begin(), indices.end(), c.indices.begin(), c.indices.end(), inserter(un, un.begin()));
309 if (un.size() != indices.size() + c.indices.size())
310 throw (std::logic_error("symmetry::add(): the same index appears in more than one child"));
316 children.push_back(c);
320 void symmetry::validate(unsigned n)
322 if (indices.upper_bound(n - 1) != indices.end())
323 throw (std::range_error("symmetry::verify(): index values are out of range"));
324 if (type != none && indices.empty()) {
325 for (unsigned i=0; i<n; i++)
334 static const symmetry & index0()
336 static ex s = (new symmetry(0))->setflag(status_flags::dynallocated);
337 return ex_to<symmetry>(s);
340 static const symmetry & index1()
342 static ex s = (new symmetry(1))->setflag(status_flags::dynallocated);
343 return ex_to<symmetry>(s);
346 static const symmetry & index2()
348 static ex s = (new symmetry(2))->setflag(status_flags::dynallocated);
349 return ex_to<symmetry>(s);
352 static const symmetry & index3()
354 static ex s = (new symmetry(3))->setflag(status_flags::dynallocated);
355 return ex_to<symmetry>(s);
358 const symmetry & not_symmetric()
360 static ex s = (new symmetry)->setflag(status_flags::dynallocated);
361 return ex_to<symmetry>(s);
364 const symmetry & symmetric2()
366 static ex s = (new symmetry(symmetry::symmetric, index0(), index1()))->setflag(status_flags::dynallocated);
367 return ex_to<symmetry>(s);
370 const symmetry & symmetric3()
372 static ex s = (new symmetry(symmetry::symmetric, index0(), index1()))->add(index2()).setflag(status_flags::dynallocated);
373 return ex_to<symmetry>(s);
376 const symmetry & symmetric4()
378 static ex s = (new symmetry(symmetry::symmetric, index0(), index1()))->add(index2()).add(index3()).setflag(status_flags::dynallocated);
379 return ex_to<symmetry>(s);
382 const symmetry & antisymmetric2()
384 static ex s = (new symmetry(symmetry::antisymmetric, index0(), index1()))->setflag(status_flags::dynallocated);
385 return ex_to<symmetry>(s);
388 const symmetry & antisymmetric3()
390 static ex s = (new symmetry(symmetry::antisymmetric, index0(), index1()))->add(index2()).setflag(status_flags::dynallocated);
391 return ex_to<symmetry>(s);
394 const symmetry & antisymmetric4()
396 static ex s = (new symmetry(symmetry::antisymmetric, index0(), index1()))->add(index2()).add(index3()).setflag(status_flags::dynallocated);
397 return ex_to<symmetry>(s);
400 class sy_is_less : public std::binary_function<ex, ex, bool> {
401 exvector::iterator v;
404 sy_is_less(exvector::iterator v_) : v(v_) {}
406 bool operator() (const ex &lh, const ex &rh) const
408 GINAC_ASSERT(is_exactly_a<symmetry>(lh));
409 GINAC_ASSERT(is_exactly_a<symmetry>(rh));
410 GINAC_ASSERT(ex_to<symmetry>(lh).indices.size() == ex_to<symmetry>(rh).indices.size());
411 std::set<unsigned>::const_iterator ait = ex_to<symmetry>(lh).indices.begin(), aitend = ex_to<symmetry>(lh).indices.end(), bit = ex_to<symmetry>(rh).indices.begin();
412 while (ait != aitend) {
413 int cmpval = v[*ait].compare(v[*bit]);
424 class sy_swap : public std::binary_function<ex, ex, void> {
425 exvector::iterator v;
430 sy_swap(exvector::iterator v_, bool &s) : v(v_), swapped(s) {}
432 void operator() (const ex &lh, const ex &rh)
434 GINAC_ASSERT(is_exactly_a<symmetry>(lh));
435 GINAC_ASSERT(is_exactly_a<symmetry>(rh));
436 GINAC_ASSERT(ex_to<symmetry>(lh).indices.size() == ex_to<symmetry>(rh).indices.size());
437 std::set<unsigned>::const_iterator ait = ex_to<symmetry>(lh).indices.begin(), aitend = ex_to<symmetry>(lh).indices.end(), bit = ex_to<symmetry>(rh).indices.begin();
438 while (ait != aitend) {
439 v[*ait].swap(v[*bit]);
446 int canonicalize(exvector::iterator v, const symmetry &symm)
448 // Less than two elements? Then do nothing
449 if (symm.indices.size() < 2)
450 return std::numeric_limits<int>::max();
452 // Canonicalize children first
453 bool something_changed = false;
455 exvector::const_iterator first = symm.children.begin(), last = symm.children.end();
456 while (first != last) {
457 GINAC_ASSERT(is_exactly_a<symmetry>(*first));
458 int child_sign = canonicalize(v, ex_to<symmetry>(*first));
461 if (child_sign != std::numeric_limits<int>::max()) {
462 something_changed = true;
468 // Now reorder the children
469 first = symm.children.begin();
471 case symmetry::symmetric:
472 // Sort the children in ascending order
473 shaker_sort(first, last, sy_is_less(v), sy_swap(v, something_changed));
475 case symmetry::antisymmetric:
476 // Sort the children in ascending order, keeping track of the signum
477 sign *= permutation_sign(first, last, sy_is_less(v), sy_swap(v, something_changed));
481 case symmetry::cyclic:
482 // Permute the smallest child to the front
483 cyclic_permutation(first, last, min_element(first, last, sy_is_less(v)), sy_swap(v, something_changed));
488 return something_changed ? sign : std::numeric_limits<int>::max();
492 // Symmetrize/antisymmetrize over a vector of objects
493 static ex symm(const ex & e, exvector::const_iterator first, exvector::const_iterator last, bool asymmetric)
495 // Need at least 2 objects for this operation
496 unsigned num = last - first;
500 // Transform object vector to a lst (for subs())
501 lst orig_lst(first, last);
503 // Create index vectors for permutation
504 unsigned *iv = new unsigned[num], *iv2;
505 for (unsigned i=0; i<num; i++)
507 iv2 = (asymmetric ? new unsigned[num] : NULL);
509 // Loop over all permutations (the first permutation, which is the
510 // identity, is unrolled)
513 while (std::next_permutation(iv, iv + num)) {
515 for (unsigned i=0; i<num; i++)
516 new_lst.append(orig_lst.op(iv[i]));
517 ex term = e.subs(orig_lst, new_lst, subs_options::no_pattern|subs_options::no_index_renaming);
519 memcpy(iv2, iv, num * sizeof(unsigned));
520 term *= permutation_sign(iv2, iv2 + num);
522 sum_v.push_back(term);
524 ex sum = (new add(sum_v))->setflag(status_flags::dynallocated);
529 return sum / factorial(numeric(num));
532 ex symmetrize(const ex & e, exvector::const_iterator first, exvector::const_iterator last)
534 return symm(e, first, last, false);
537 ex antisymmetrize(const ex & e, exvector::const_iterator first, exvector::const_iterator last)
539 return symm(e, first, last, true);
542 ex symmetrize_cyclic(const ex & e, exvector::const_iterator first, exvector::const_iterator last)
544 // Need at least 2 objects for this operation
545 unsigned num = last - first;
549 // Transform object vector to a lst (for subs())
550 lst orig_lst(first, last);
551 lst new_lst = orig_lst;
553 // Loop over all cyclic permutations (the first permutation, which is
554 // the identity, is unrolled)
556 for (unsigned i=0; i<num-1; i++) {
557 ex perm = new_lst.op(0);
558 new_lst.remove_first().append(perm);
559 sum += e.subs(orig_lst, new_lst, subs_options::no_pattern|subs_options::no_index_renaming);
564 /** Symmetrize expression over a list of objects (symbols, indices). */
565 ex ex::symmetrize(const lst & l) const
567 exvector v(l.begin(), l.end());
568 return symm(*this, v.begin(), v.end(), false);
571 /** Antisymmetrize expression over a list of objects (symbols, indices). */
572 ex ex::antisymmetrize(const lst & l) const
574 exvector v(l.begin(), l.end());
575 return symm(*this, v.begin(), v.end(), true);
578 /** Symmetrize expression by cyclic permutation over a list of objects
579 * (symbols, indices). */
580 ex ex::symmetrize_cyclic(const lst & l) const
582 exvector v(l.begin(), l.end());
583 return GiNaC::symmetrize_cyclic(*this, v.begin(), v.end());