3 * Implementation of GiNaC's symmetry definitions. */
6 * GiNaC Copyright (C) 1999-2001 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., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
27 #define DO_GINAC_ASSERT
28 #include "assertion.h"
32 #include "numeric.h" // for factorial()
40 GINAC_IMPLEMENT_REGISTERED_CLASS(symmetry, basic)
43 // default constructor, destructor, copy constructor assignment operator and helpers
46 symmetry::symmetry() : type(none)
48 debugmsg("symmetry default constructor", LOGLEVEL_CONSTRUCT);
49 tinfo_key = TINFO_symmetry;
52 void symmetry::copy(const symmetry & other)
54 inherited::copy(other);
56 indices = other.indices;
57 children = other.children;
60 DEFAULT_DESTROY(symmetry)
66 symmetry::symmetry(unsigned i) : type(none)
68 debugmsg("symmetry constructor from unsigned", LOGLEVEL_CONSTRUCT);
70 tinfo_key = TINFO_symmetry;
73 symmetry::symmetry(symmetry_type t, const symmetry &c1, const symmetry &c2) : type(t)
75 debugmsg("symmetry constructor from symmetry_type,symmetry &,symmetry &", LOGLEVEL_CONSTRUCT);
77 tinfo_key = TINFO_symmetry;
84 /** Construct object from archive_node. */
85 symmetry::symmetry(const archive_node &n, const lst &sym_lst) : inherited(n, sym_lst)
87 debugmsg("symmetry ctor from archive_node", LOGLEVEL_CONSTRUCT);
90 if (!(n.find_unsigned("type", t)))
91 throw (std::runtime_error("unknown symmetry type in archive"));
92 type = (symmetry_type)t;
97 if (n.find_ex("child", e, sym_lst, i))
98 add(ex_to_symmetry(e));
107 if (n.find_unsigned("index", u, i))
116 /** Archive the object. */
117 void symmetry::archive(archive_node &n) const
119 inherited::archive(n);
121 n.add_unsigned("type", type);
123 if (children.empty()) {
124 std::set<unsigned>::const_iterator i = indices.begin(), iend = indices.end();
126 n.add_unsigned("index", *i);
130 exvector::const_iterator i = children.begin(), iend = children.end();
132 n.add_ex("child", *i);
138 DEFAULT_UNARCHIVE(symmetry)
141 // functions overriding virtual functions from bases classes
144 int symmetry::compare_same_type(const basic & other) const
146 GINAC_ASSERT(is_of_type(other, symmetry));
147 const symmetry &o = static_cast<const symmetry &>(other);
149 // All symmetry trees are equal. They are not supposed to appear in
150 // ordinary expressions anyway...
154 void symmetry::print(const print_context & c, unsigned level = 0) const
156 debugmsg("symmetry print", LOGLEVEL_PRINT);
158 if (children.empty()) {
159 if (indices.size() > 0)
160 c.s << *(indices.begin());
163 case none: c.s << '!'; break;
164 case symmetric: c.s << '+'; break;
165 case antisymmetric: c.s << '-'; break;
166 case cyclic: c.s << '@'; break;
167 default: c.s << '?'; break;
170 for (unsigned i=0; i<children.size(); i++) {
171 children[i].print(c);
172 if (i != children.size() - 1)
180 // non-virtual functions in this class
183 symmetry &symmetry::add(const symmetry &c)
185 // All children must have the same number of indices
186 if (type != none && !children.empty()) {
187 GINAC_ASSERT(is_ex_exactly_of_type(children[0], symmetry));
188 if (ex_to_symmetry(children[0]).indices.size() != c.indices.size())
189 throw (std::logic_error("symmetry:add(): children must have same number of indices"));
192 // Compute union of indices and check whether the two sets are disjoint
193 std::set<unsigned> un;
194 set_union(indices.begin(), indices.end(), c.indices.begin(), c.indices.end(), inserter(un, un.begin()));
195 if (un.size() != indices.size() + c.indices.size())
196 throw (std::logic_error("symmetry::add(): the same index appears in more than one child"));
202 children.push_back(c);
206 void symmetry::validate(unsigned n)
208 if (indices.upper_bound(n - 1) != indices.end())
209 throw (std::range_error("symmetry::verify(): index values are out of range"));
210 if (type != none && indices.empty()) {
211 for (unsigned i=0; i<n; i++)
220 class sy_is_less : public std::binary_function<ex, ex, bool> {
221 exvector::iterator v;
224 sy_is_less(exvector::iterator v_) : v(v_) {}
226 bool operator() (const ex &lh, const ex &rh) const
228 GINAC_ASSERT(is_ex_exactly_of_type(lh, symmetry));
229 GINAC_ASSERT(is_ex_exactly_of_type(rh, symmetry));
230 GINAC_ASSERT(ex_to_symmetry(lh).indices.size() == ex_to_symmetry(rh).indices.size());
231 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();
232 while (ait != aitend) {
233 int cmpval = v[*ait].compare(v[*bit]);
244 class sy_swap : public std::binary_function<ex, ex, void> {
245 exvector::iterator v;
250 sy_swap(exvector::iterator v_, bool &s) : v(v_), swapped(s) {}
252 void operator() (const ex &lh, const ex &rh)
254 GINAC_ASSERT(is_ex_exactly_of_type(lh, symmetry));
255 GINAC_ASSERT(is_ex_exactly_of_type(rh, symmetry));
256 GINAC_ASSERT(ex_to_symmetry(lh).indices.size() == ex_to_symmetry(rh).indices.size());
257 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();
258 while (ait != aitend) {
259 v[*ait].swap(v[*bit]);
266 int canonicalize(exvector::iterator v, const symmetry &symm)
268 // No children? Then do nothing
269 if (symm.children.empty())
272 // Canonicalize children first
273 bool something_changed = false;
275 exvector::const_iterator first = symm.children.begin(), last = symm.children.end();
276 while (first != last) {
277 GINAC_ASSERT(is_ex_exactly_of_type(*first, symmetry));
278 int child_sign = canonicalize(v, ex_to_symmetry(*first));
281 if (child_sign != INT_MAX) {
282 something_changed = true;
288 // Now reorder the children
289 first = symm.children.begin();
291 case symmetry::symmetric:
292 shaker_sort(first, last, sy_is_less(v), sy_swap(v, something_changed));
294 case symmetry::antisymmetric:
295 sign *= permutation_sign(first, last, sy_is_less(v), sy_swap(v, something_changed));
297 case symmetry::cyclic:
298 cyclic_permutation(first, last, min_element(first, last, sy_is_less(v)), sy_swap(v, something_changed));
303 return something_changed ? sign : INT_MAX;
307 // Symmetrize/antisymmetrize over a vector of objects
308 static ex symm(const ex & e, exvector::const_iterator first, exvector::const_iterator last, bool asymmetric)
310 // Need at least 2 objects for this operation
311 int num = last - first;
315 // Transform object vector to a list
317 iv_lst.insert(iv_lst.begin(), first, last);
318 lst orig_lst(iv_lst, true);
320 // Create index vectors for permutation
321 unsigned *iv = new unsigned[num], *iv2;
322 for (unsigned i=0; i<num; i++)
324 iv2 = (asymmetric ? new unsigned[num] : NULL);
326 // Loop over all permutations (the first permutation, which is the
327 // identity, is unrolled)
329 while (std::next_permutation(iv, iv + num)) {
331 for (unsigned i=0; i<num; i++)
332 new_lst.append(orig_lst.op(iv[i]));
333 ex term = e.subs(orig_lst, new_lst);
335 memcpy(iv2, iv, num * sizeof(unsigned));
336 term *= permutation_sign(iv2, iv2 + num);
344 return sum / factorial(numeric(num));
347 ex symmetrize(const ex & e, exvector::const_iterator first, exvector::const_iterator last)
349 return symm(e, first, last, false);
352 ex antisymmetrize(const ex & e, exvector::const_iterator first, exvector::const_iterator last)
354 return symm(e, first, last, true);
357 ex symmetrize_cyclic(const ex & e, exvector::const_iterator first, exvector::const_iterator last)
359 // Need at least 2 objects for this operation
360 int num = last - first;
364 // Transform object vector to a list
366 iv_lst.insert(iv_lst.begin(), first, last);
367 lst orig_lst(iv_lst, true);
368 lst new_lst = orig_lst;
370 // Loop over all cyclic permutations (the first permutation, which is
371 // the identity, is unrolled)
373 for (unsigned i=0; i<num-1; i++) {
374 ex perm = new_lst.op(0);
375 new_lst.remove_first().append(perm);
376 sum += e.subs(orig_lst, new_lst);
381 /** Symmetrize expression over a list of objects (symbols, indices). */
382 ex ex::symmetrize(const lst & l) const
386 for (unsigned i=0; i<l.nops(); i++)
387 v.push_back(l.op(i));
388 return symm(*this, v.begin(), v.end(), false);
391 /** Antisymmetrize expression over a list of objects (symbols, indices). */
392 ex ex::antisymmetrize(const lst & l) const
396 for (unsigned i=0; i<l.nops(); i++)
397 v.push_back(l.op(i));
398 return symm(*this, v.begin(), v.end(), true);
401 /** Symmetrize expression by cyclic permutation over a list of objects
402 * (symbols, indices). */
403 ex ex::symmetrize_cyclic(const lst & l) const
407 for (unsigned i=0; i<l.nops(); i++)
408 v.push_back(l.op(i));
409 return GiNaC::symmetrize_cyclic(*this, v.begin(), v.end());