3 * Implementation of GiNaC's symmetry definitions. */
6 * GiNaC Copyright (C) 1999-2006 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
29 #include "numeric.h" // for factorial()
30 #include "operators.h"
36 GINAC_IMPLEMENT_REGISTERED_CLASS_OPT(symmetry, basic,
37 print_func<print_context>(&symmetry::do_print).
38 print_func<print_tree>(&symmetry::do_print_tree))
41 Some notes about the structure of a symmetry tree:
42 - The leaf nodes of the tree are of type "none", have one index, and no
43 children (of course). They are constructed by the symmetry(unsigned)
45 - Leaf nodes are the only nodes that only have one index.
46 - Container nodes contain two or more children. The "indices" set member
47 is the set union of the index sets of all children, and the "children"
48 vector stores the children themselves.
49 - The index set of each child of a "symm", "anti" or "cycl" node must
50 have the same size. It follows that the children of such a node are
51 either all leaf nodes, or all container nodes with two or more indices.
55 // default constructor
58 symmetry::symmetry() : inherited(&symmetry::tinfo_static), type(none)
60 setflag(status_flags::evaluated | status_flags::expanded);
67 symmetry::symmetry(unsigned i) : inherited(&symmetry::tinfo_static), type(none)
70 setflag(status_flags::evaluated | status_flags::expanded);
73 symmetry::symmetry(symmetry_type t, const symmetry &c1, const symmetry &c2) : inherited(&symmetry::tinfo_static), type(t)
76 setflag(status_flags::evaluated | status_flags::expanded);
83 /** Construct object from archive_node. */
84 symmetry::symmetry(const archive_node &n, lst &sym_lst) : inherited(n, sym_lst)
87 if (!(n.find_unsigned("type", t)))
88 throw (std::runtime_error("unknown symmetry type in archive"));
89 type = (symmetry_type)t;
94 if (n.find_ex("child", e, sym_lst, i))
95 add(ex_to<symmetry>(e));
104 if (n.find_unsigned("index", u, i))
113 /** Archive the object. */
114 void symmetry::archive(archive_node &n) const
116 inherited::archive(n);
118 n.add_unsigned("type", type);
120 if (children.empty()) {
121 std::set<unsigned>::const_iterator i = indices.begin(), iend = indices.end();
123 n.add_unsigned("index", *i);
127 exvector::const_iterator i = children.begin(), iend = children.end();
129 n.add_ex("child", *i);
135 DEFAULT_UNARCHIVE(symmetry)
138 // functions overriding virtual functions from base classes
141 int symmetry::compare_same_type(const basic & other) const
143 GINAC_ASSERT(is_a<symmetry>(other));
145 // For archiving purposes we need to have an ordering of symmetries.
146 const symmetry &othersymm = ex_to<symmetry>(other);
149 if (type > othersymm.type)
151 if (type < othersymm.type)
154 // Compare the index set.
155 size_t this_size = indices.size();
156 size_t that_size = othersymm.indices.size();
157 if (this_size > that_size)
159 if (this_size < that_size)
161 typedef std::set<unsigned>::iterator set_it;
162 set_it end = indices.end();
163 for (set_it i=indices.begin(),j=othersymm.indices.begin(); i!=end; ++i,++j) {
170 // Compare the children.
171 if (children.size() > othersymm.children.size())
173 if (children.size() < othersymm.children.size())
175 for (size_t i=0; i<children.size(); ++i) {
176 int cmpval = ex_to<symmetry>(children[i])
177 .compare_same_type(ex_to<symmetry>(othersymm.children[i]));
185 unsigned symmetry::calchash() const
187 unsigned v = golden_ratio_hash((p_int)tinfo());
191 v ^= *(indices.begin());
193 for (exvector::const_iterator i=children.begin(); i!=children.end(); ++i)
200 if (flags & status_flags::evaluated) {
201 setflag(status_flags::hash_calculated);
208 void symmetry::do_print(const print_context & c, unsigned level) const
210 if (children.empty()) {
211 if (indices.size() > 0)
212 c.s << *(indices.begin());
217 case none: c.s << '!'; break;
218 case symmetric: c.s << '+'; break;
219 case antisymmetric: c.s << '-'; break;
220 case cyclic: c.s << '@'; break;
221 default: c.s << '?'; break;
224 size_t num = children.size();
225 for (size_t i=0; i<num; i++) {
226 children[i].print(c);
234 void symmetry::do_print_tree(const print_tree & c, unsigned level) const
236 c.s << std::string(level, ' ') << class_name() << " @" << this
237 << std::hex << ", hash=0x" << hashvalue << ", flags=0x" << flags << std::dec
241 case none: c.s << "none"; break;
242 case symmetric: c.s << "symm"; break;
243 case antisymmetric: c.s << "anti"; break;
244 case cyclic: c.s << "cycl"; break;
245 default: c.s << "<unknown>"; break;
248 c.s << ", indices=(";
249 if (!indices.empty()) {
250 std::set<unsigned>::const_iterator i = indices.begin(), end = indices.end();
258 exvector::const_iterator i = children.begin(), end = children.end();
260 i->print(c, level + c.delta_indent);
266 // non-virtual functions in this class
269 bool symmetry::has_cyclic() const
274 for (exvector::const_iterator i=children.begin(); i!=children.end(); ++i)
275 if (ex_to<symmetry>(*i).has_cyclic())
281 symmetry &symmetry::add(const symmetry &c)
283 // All children must have the same number of indices
284 if (type != none && !children.empty()) {
285 GINAC_ASSERT(is_exactly_a<symmetry>(children[0]));
286 if (ex_to<symmetry>(children[0]).indices.size() != c.indices.size())
287 throw (std::logic_error("symmetry:add(): children must have same number of indices"));
290 // Compute union of indices and check whether the two sets are disjoint
291 std::set<unsigned> un;
292 set_union(indices.begin(), indices.end(), c.indices.begin(), c.indices.end(), inserter(un, un.begin()));
293 if (un.size() != indices.size() + c.indices.size())
294 throw (std::logic_error("symmetry::add(): the same index appears in more than one child"));
300 children.push_back(c);
304 void symmetry::validate(unsigned n)
306 if (indices.upper_bound(n - 1) != indices.end())
307 throw (std::range_error("symmetry::verify(): index values are out of range"));
308 if (type != none && indices.empty()) {
309 for (unsigned i=0; i<n; i++)
318 static const symmetry & index0()
320 static ex s = (new symmetry(0))->setflag(status_flags::dynallocated);
321 return ex_to<symmetry>(s);
324 static const symmetry & index1()
326 static ex s = (new symmetry(1))->setflag(status_flags::dynallocated);
327 return ex_to<symmetry>(s);
330 static const symmetry & index2()
332 static ex s = (new symmetry(2))->setflag(status_flags::dynallocated);
333 return ex_to<symmetry>(s);
336 static const symmetry & index3()
338 static ex s = (new symmetry(3))->setflag(status_flags::dynallocated);
339 return ex_to<symmetry>(s);
342 const symmetry & not_symmetric()
344 static ex s = (new symmetry)->setflag(status_flags::dynallocated);
345 return ex_to<symmetry>(s);
348 const symmetry & symmetric2()
350 static ex s = (new symmetry(symmetry::symmetric, index0(), index1()))->setflag(status_flags::dynallocated);
351 return ex_to<symmetry>(s);
354 const symmetry & symmetric3()
356 static ex s = (new symmetry(symmetry::symmetric, index0(), index1()))->add(index2()).setflag(status_flags::dynallocated);
357 return ex_to<symmetry>(s);
360 const symmetry & symmetric4()
362 static ex s = (new symmetry(symmetry::symmetric, index0(), index1()))->add(index2()).add(index3()).setflag(status_flags::dynallocated);
363 return ex_to<symmetry>(s);
366 const symmetry & antisymmetric2()
368 static ex s = (new symmetry(symmetry::antisymmetric, index0(), index1()))->setflag(status_flags::dynallocated);
369 return ex_to<symmetry>(s);
372 const symmetry & antisymmetric3()
374 static ex s = (new symmetry(symmetry::antisymmetric, index0(), index1()))->add(index2()).setflag(status_flags::dynallocated);
375 return ex_to<symmetry>(s);
378 const symmetry & antisymmetric4()
380 static ex s = (new symmetry(symmetry::antisymmetric, index0(), index1()))->add(index2()).add(index3()).setflag(status_flags::dynallocated);
381 return ex_to<symmetry>(s);
384 class sy_is_less : public std::binary_function<ex, ex, bool> {
385 exvector::iterator v;
388 sy_is_less(exvector::iterator v_) : v(v_) {}
390 bool operator() (const ex &lh, const ex &rh) const
392 GINAC_ASSERT(is_exactly_a<symmetry>(lh));
393 GINAC_ASSERT(is_exactly_a<symmetry>(rh));
394 GINAC_ASSERT(ex_to<symmetry>(lh).indices.size() == ex_to<symmetry>(rh).indices.size());
395 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();
396 while (ait != aitend) {
397 int cmpval = v[*ait].compare(v[*bit]);
408 class sy_swap : public std::binary_function<ex, ex, void> {
409 exvector::iterator v;
414 sy_swap(exvector::iterator v_, bool &s) : v(v_), swapped(s) {}
416 void operator() (const ex &lh, const ex &rh)
418 GINAC_ASSERT(is_exactly_a<symmetry>(lh));
419 GINAC_ASSERT(is_exactly_a<symmetry>(rh));
420 GINAC_ASSERT(ex_to<symmetry>(lh).indices.size() == ex_to<symmetry>(rh).indices.size());
421 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();
422 while (ait != aitend) {
423 v[*ait].swap(v[*bit]);
430 int canonicalize(exvector::iterator v, const symmetry &symm)
432 // Less than two elements? Then do nothing
433 if (symm.indices.size() < 2)
436 // Canonicalize children first
437 bool something_changed = false;
439 exvector::const_iterator first = symm.children.begin(), last = symm.children.end();
440 while (first != last) {
441 GINAC_ASSERT(is_exactly_a<symmetry>(*first));
442 int child_sign = canonicalize(v, ex_to<symmetry>(*first));
445 if (child_sign != INT_MAX) {
446 something_changed = true;
452 // Now reorder the children
453 first = symm.children.begin();
455 case symmetry::symmetric:
456 // Sort the children in ascending order
457 shaker_sort(first, last, sy_is_less(v), sy_swap(v, something_changed));
459 case symmetry::antisymmetric:
460 // Sort the children in ascending order, keeping track of the signum
461 sign *= permutation_sign(first, last, sy_is_less(v), sy_swap(v, something_changed));
465 case symmetry::cyclic:
466 // Permute the smallest child to the front
467 cyclic_permutation(first, last, min_element(first, last, sy_is_less(v)), sy_swap(v, something_changed));
472 return something_changed ? sign : INT_MAX;
476 // Symmetrize/antisymmetrize over a vector of objects
477 static ex symm(const ex & e, exvector::const_iterator first, exvector::const_iterator last, bool asymmetric)
479 // Need at least 2 objects for this operation
480 unsigned num = last - first;
484 // Transform object vector to a lst (for subs())
485 lst orig_lst(first, last);
487 // Create index vectors for permutation
488 unsigned *iv = new unsigned[num], *iv2;
489 for (unsigned i=0; i<num; i++)
491 iv2 = (asymmetric ? new unsigned[num] : NULL);
493 // Loop over all permutations (the first permutation, which is the
494 // identity, is unrolled)
496 while (std::next_permutation(iv, iv + num)) {
498 for (unsigned i=0; i<num; i++)
499 new_lst.append(orig_lst.op(iv[i]));
500 ex term = e.subs(orig_lst, new_lst, subs_options::no_pattern|subs_options::no_index_renaming);
502 memcpy(iv2, iv, num * sizeof(unsigned));
503 term *= permutation_sign(iv2, iv2 + num);
511 return sum / factorial(numeric(num));
514 ex symmetrize(const ex & e, exvector::const_iterator first, exvector::const_iterator last)
516 return symm(e, first, last, false);
519 ex antisymmetrize(const ex & e, exvector::const_iterator first, exvector::const_iterator last)
521 return symm(e, first, last, true);
524 ex symmetrize_cyclic(const ex & e, exvector::const_iterator first, exvector::const_iterator last)
526 // Need at least 2 objects for this operation
527 unsigned num = last - first;
531 // Transform object vector to a lst (for subs())
532 lst orig_lst(first, last);
533 lst new_lst = orig_lst;
535 // Loop over all cyclic permutations (the first permutation, which is
536 // the identity, is unrolled)
538 for (unsigned i=0; i<num-1; i++) {
539 ex perm = new_lst.op(0);
540 new_lst.remove_first().append(perm);
541 sum += e.subs(orig_lst, new_lst, subs_options::no_pattern|subs_options::no_index_renaming);
546 /** Symmetrize expression over a list of objects (symbols, indices). */
547 ex ex::symmetrize(const lst & l) const
549 exvector v(l.begin(), l.end());
550 return symm(*this, v.begin(), v.end(), false);
553 /** Antisymmetrize expression over a list of objects (symbols, indices). */
554 ex ex::antisymmetrize(const lst & l) const
556 exvector v(l.begin(), l.end());
557 return symm(*this, v.begin(), v.end(), true);
560 /** Symmetrize expression by cyclic permutation over a list of objects
561 * (symbols, indices). */
562 ex ex::symmetrize_cyclic(const lst & l) const
564 exvector v(l.begin(), l.end());
565 return GiNaC::symmetrize_cyclic(*this, v.begin(), v.end());