Add support for Texinfo-5.0.
[ginac.git] / ginac / symmetry.cpp
1 /** @file symmetry.cpp
2  *
3  *  Implementation of GiNaC's symmetry definitions. */
4
5 /*
6  *  GiNaC Copyright (C) 1999-2011 Johannes Gutenberg University Mainz, Germany
7  *
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.
12  *
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.
17  *
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
21  */
22
23 #include "symmetry.h"
24 #include "lst.h"
25 #include "add.h"
26 #include "numeric.h" // for factorial()
27 #include "operators.h"
28 #include "archive.h"
29 #include "utils.h"
30 #include "hash_seed.h"
31
32 #include <functional>
33 #include <iostream>
34 #include <limits>
35 #include <stdexcept>
36
37 namespace GiNaC {
38
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))
42
43 /*
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)
47       constructor.
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.
55 */
56
57 //////////
58 // default constructor
59 //////////
60
61 symmetry::symmetry() :  type(none)
62 {
63         setflag(status_flags::evaluated | status_flags::expanded);
64 }
65
66 //////////
67 // other constructors
68 //////////
69
70 symmetry::symmetry(unsigned i) :  type(none)
71 {
72         indices.insert(i);
73         setflag(status_flags::evaluated | status_flags::expanded);
74 }
75
76 symmetry::symmetry(symmetry_type t, const symmetry &c1, const symmetry &c2) :  type(t)
77 {
78         add(c1); add(c2);
79         setflag(status_flags::evaluated | status_flags::expanded);
80 }
81
82 //////////
83 // archiving
84 //////////
85
86 /** Construct object from archive_node. */
87 void symmetry::read_archive(const archive_node &n, lst &sym_lst)
88 {
89         inherited::read_archive(n, sym_lst);
90         unsigned t;
91         if (!(n.find_unsigned("type", t)))
92                 throw (std::runtime_error("unknown symmetry type in archive"));
93         type = (symmetry_type)t;
94
95         unsigned i = 0;
96         while (true) {
97                 ex e;
98                 if (n.find_ex("child", e, sym_lst, i))
99                         add(ex_to<symmetry>(e));
100                 else
101                         break;
102                 i++;
103         }
104
105         if (i == 0) {
106                 while (true) {
107                         unsigned u;
108                         if (n.find_unsigned("index", u, i))
109                                 indices.insert(u);
110                         else
111                                 break;
112                         i++;
113                 }
114         }
115 }
116 GINAC_BIND_UNARCHIVER(symmetry);
117
118 /** Archive the object. */
119 void symmetry::archive(archive_node &n) const
120 {
121         inherited::archive(n);
122
123         n.add_unsigned("type", type);
124
125         if (children.empty()) {
126                 std::set<unsigned>::const_iterator i = indices.begin(), iend = indices.end();
127                 while (i != iend) {
128                         n.add_unsigned("index", *i);
129                         i++;
130                 }
131         } else {
132                 exvector::const_iterator i = children.begin(), iend = children.end();
133                 while (i != iend) {
134                         n.add_ex("child", *i);
135                         i++;
136                 }
137         }
138 }
139
140 //////////
141 // functions overriding virtual functions from base classes
142 //////////
143
144 int symmetry::compare_same_type(const basic & other) const
145 {
146         GINAC_ASSERT(is_a<symmetry>(other));
147
148         // For archiving purposes we need to have an ordering of symmetries.
149         const symmetry &othersymm = ex_to<symmetry>(other);
150
151         // Compare type.
152         if (type > othersymm.type)
153                 return 1;
154         if (type < othersymm.type)
155                 return -1;
156
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)
161                 return 1;
162         if (this_size < that_size)
163                 return -1;
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) {
167                 if(*i < *j)
168                         return 1;
169                 if(*i > *j)
170                         return -1;
171         }
172
173         // Compare the children.
174         if (children.size() > othersymm.children.size())
175                 return 1;
176         if (children.size() < othersymm.children.size())
177                 return -1;
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]));
181                 if (cmpval)
182                         return cmpval;
183         }
184
185         return 0;
186 }
187
188 unsigned symmetry::calchash() const
189 {
190         unsigned v = make_hash_seed(typeid(*this));
191
192         if (type == none) {
193                 v = rotate_left(v);
194                 if (!indices.empty())
195                         v ^= *(indices.begin());
196         } else {
197                 for (exvector::const_iterator i=children.begin(); i!=children.end(); ++i)
198                 {
199                         v = rotate_left(v);
200                         v ^= i->gethash();
201                 }
202         }
203
204         if (flags & status_flags::evaluated) {
205                 setflag(status_flags::hash_calculated);
206                 hashvalue = v;
207         }
208
209         return v;
210 }
211
212 void symmetry::do_print(const print_context & c, unsigned level) const
213 {
214         if (children.empty()) {
215                 if (indices.size() > 0)
216                         c.s << *(indices.begin());
217                 else
218                         c.s << "none";
219         } else {
220                 switch (type) {
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;
226                 }
227                 c.s << '(';
228                 size_t num = children.size();
229                 for (size_t i=0; i<num; i++) {
230                         children[i].print(c);
231                         if (i != num - 1)
232                                 c.s << ",";
233                 }
234                 c.s << ')';
235         }
236 }
237
238 void symmetry::do_print_tree(const print_tree & c, unsigned level) const
239 {
240         c.s << std::string(level, ' ') << class_name() << " @" << this
241             << std::hex << ", hash=0x" << hashvalue << ", flags=0x" << flags << std::dec
242             << ", type=";
243
244         switch (type) {
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;
250         }
251
252         c.s << ", indices=(";
253         if (!indices.empty()) {
254                 std::set<unsigned>::const_iterator i = indices.begin(), end = indices.end();
255                 --end;
256                 while (i != end)
257                         c.s << *i++ << ",";
258                 c.s << *i;
259         }
260         c.s << ")\n";
261
262         exvector::const_iterator i = children.begin(), end = children.end();
263         while (i != end) {
264                 i->print(c, level + c.delta_indent);
265                 ++i;
266         }
267 }
268
269 //////////
270 // non-virtual functions in this class
271 //////////
272
273 bool symmetry::has_nonsymmetric() const
274 {
275         if (type == antisymmetric || type == cyclic)
276                 return true;
277
278         for (exvector::const_iterator i=children.begin(); i!=children.end(); ++i)
279                 if (ex_to<symmetry>(*i).has_nonsymmetric())
280                         return true;
281
282         return false;
283 }
284
285 bool symmetry::has_cyclic() const
286 {
287         if (type == cyclic)
288                 return true;
289
290         for (exvector::const_iterator i=children.begin(); i!=children.end(); ++i)
291                 if (ex_to<symmetry>(*i).has_cyclic())
292                         return true;
293
294         return false;
295 }
296
297 symmetry &symmetry::add(const symmetry &c)
298 {
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"));
304         }
305
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"));
311
312         // Set new index set
313         indices.swap(un);
314
315         // Add child node
316         children.push_back(c);
317         return *this;
318 }
319
320 void symmetry::validate(unsigned n)
321 {
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++)
326                         add(i);
327         }
328 }
329
330 //////////
331 // global functions
332 //////////
333
334 static const symmetry & index0()
335 {
336         static ex s = (new symmetry(0))->setflag(status_flags::dynallocated);
337         return ex_to<symmetry>(s);
338 }
339
340 static const symmetry & index1()
341 {
342         static ex s = (new symmetry(1))->setflag(status_flags::dynallocated);
343         return ex_to<symmetry>(s);
344 }
345
346 static const symmetry & index2()
347 {
348         static ex s = (new symmetry(2))->setflag(status_flags::dynallocated);
349         return ex_to<symmetry>(s);
350 }
351
352 static const symmetry & index3()
353 {
354         static ex s = (new symmetry(3))->setflag(status_flags::dynallocated);
355         return ex_to<symmetry>(s);
356 }
357
358 const symmetry & not_symmetric()
359 {
360         static ex s = (new symmetry)->setflag(status_flags::dynallocated);
361         return ex_to<symmetry>(s);
362 }
363
364 const symmetry & symmetric2()
365 {
366         static ex s = (new symmetry(symmetry::symmetric, index0(), index1()))->setflag(status_flags::dynallocated);
367         return ex_to<symmetry>(s);
368 }
369
370 const symmetry & symmetric3()
371 {
372         static ex s = (new symmetry(symmetry::symmetric, index0(), index1()))->add(index2()).setflag(status_flags::dynallocated);
373         return ex_to<symmetry>(s);
374 }
375
376 const symmetry & symmetric4()
377 {
378         static ex s = (new symmetry(symmetry::symmetric, index0(), index1()))->add(index2()).add(index3()).setflag(status_flags::dynallocated);
379         return ex_to<symmetry>(s);
380 }
381
382 const symmetry & antisymmetric2()
383 {
384         static ex s = (new symmetry(symmetry::antisymmetric, index0(), index1()))->setflag(status_flags::dynallocated);
385         return ex_to<symmetry>(s);
386 }
387
388 const symmetry & antisymmetric3()
389 {
390         static ex s = (new symmetry(symmetry::antisymmetric, index0(), index1()))->add(index2()).setflag(status_flags::dynallocated);
391         return ex_to<symmetry>(s);
392 }
393
394 const symmetry & antisymmetric4()
395 {
396         static ex s = (new symmetry(symmetry::antisymmetric, index0(), index1()))->add(index2()).add(index3()).setflag(status_flags::dynallocated);
397         return ex_to<symmetry>(s);
398 }
399
400 class sy_is_less : public std::binary_function<ex, ex, bool> {
401         exvector::iterator v;
402
403 public:
404         sy_is_less(exvector::iterator v_) : v(v_) {}
405
406         bool operator() (const ex &lh, const ex &rh) const
407         {
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]);
414                         if (cmpval < 0)
415                                 return true;
416                         else if (cmpval > 0)
417                                 return false;
418                         ++ait; ++bit;
419                 }
420                 return false;
421         }
422 };
423
424 class sy_swap : public std::binary_function<ex, ex, void> {
425         exvector::iterator v;
426
427 public:
428         bool &swapped;
429
430         sy_swap(exvector::iterator v_, bool &s) : v(v_), swapped(s) {}
431
432         void operator() (const ex &lh, const ex &rh)
433         {
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]);
440                         ++ait; ++bit;
441                 }
442                 swapped = true;
443         }
444 };
445
446 int canonicalize(exvector::iterator v, const symmetry &symm)
447 {
448         // Less than two elements? Then do nothing
449         if (symm.indices.size() < 2)
450                 return std::numeric_limits<int>::max();
451
452         // Canonicalize children first
453         bool something_changed = false;
454         int sign = 1;
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));
459                 if (child_sign == 0)
460                         return 0;
461                 if (child_sign != std::numeric_limits<int>::max()) {
462                         something_changed = true;
463                         sign *= child_sign;
464                 }
465                 first++;
466         }
467
468         // Now reorder the children
469         first = symm.children.begin();
470         switch (symm.type) {
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));
474                         break;
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));
478                         if (sign == 0)
479                                 return 0;
480                         break;
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));
484                         break;
485                 default:
486                         break;
487         }
488         return something_changed ? sign : std::numeric_limits<int>::max();
489 }
490
491
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)
494 {
495         // Need at least 2 objects for this operation
496         unsigned num = last - first;
497         if (num < 2)
498                 return e;
499
500         // Transform object vector to a lst (for subs())
501         lst orig_lst(first, last);
502
503         // Create index vectors for permutation
504         unsigned *iv = new unsigned[num], *iv2;
505         for (unsigned i=0; i<num; i++)
506                 iv[i] = i;
507         iv2 = (asymmetric ? new unsigned[num] : NULL);
508
509         // Loop over all permutations (the first permutation, which is the
510         // identity, is unrolled)
511         exvector sum_v;
512         sum_v.push_back(e);
513         while (std::next_permutation(iv, iv + num)) {
514                 lst new_lst;
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);
518                 if (asymmetric) {
519                         memcpy(iv2, iv, num * sizeof(unsigned));
520                         term *= permutation_sign(iv2, iv2 + num);
521                 }
522                 sum_v.push_back(term);
523         }
524         ex sum = (new add(sum_v))->setflag(status_flags::dynallocated);
525
526         delete[] iv;
527         delete[] iv2;
528
529         return sum / factorial(numeric(num));
530 }
531
532 ex symmetrize(const ex & e, exvector::const_iterator first, exvector::const_iterator last)
533 {
534         return symm(e, first, last, false);
535 }
536
537 ex antisymmetrize(const ex & e, exvector::const_iterator first, exvector::const_iterator last)
538 {
539         return symm(e, first, last, true);
540 }
541
542 ex symmetrize_cyclic(const ex & e, exvector::const_iterator first, exvector::const_iterator last)
543 {
544         // Need at least 2 objects for this operation
545         unsigned num = last - first;
546         if (num < 2)
547                 return e;
548
549         // Transform object vector to a lst (for subs())
550         lst orig_lst(first, last);
551         lst new_lst = orig_lst;
552
553         // Loop over all cyclic permutations (the first permutation, which is
554         // the identity, is unrolled)
555         ex sum = e;
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);
560         }
561         return sum / num;
562 }
563
564 /** Symmetrize expression over a list of objects (symbols, indices). */
565 ex ex::symmetrize(const lst & l) const
566 {
567         exvector v(l.begin(), l.end());
568         return symm(*this, v.begin(), v.end(), false);
569 }
570
571 /** Antisymmetrize expression over a list of objects (symbols, indices). */
572 ex ex::antisymmetrize(const lst & l) const
573 {
574         exvector v(l.begin(), l.end());
575         return symm(*this, v.begin(), v.end(), true);
576 }
577
578 /** Symmetrize expression by cyclic permutation over a list of objects
579  *  (symbols, indices). */
580 ex ex::symmetrize_cyclic(const lst & l) const
581 {
582         exvector v(l.begin(), l.end());
583         return GiNaC::symmetrize_cyclic(*this, v.begin(), v.end());
584 }
585
586 } // namespace GiNaC