90ac3a425a35bfd6b2158e60b40eb8e8d36277fb
[ginac.git] / ginac / idx.cpp
1 /** @file idx.cpp
2  *
3  *  Implementation of GiNaC's indices. */
4
5 /*
6  *  GiNaC Copyright (C) 1999-2001 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., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
21  */
22
23 #include <stdexcept>
24
25 #include "idx.h"
26 #include "symbol.h"
27 #include "lst.h"
28 #include "print.h"
29 #include "archive.h"
30 #include "utils.h"
31 #include "debugmsg.h"
32
33 namespace GiNaC {
34
35 GINAC_IMPLEMENT_REGISTERED_CLASS(idx, basic)
36 GINAC_IMPLEMENT_REGISTERED_CLASS(varidx, idx)
37 GINAC_IMPLEMENT_REGISTERED_CLASS(spinidx, varidx)
38
39 //////////
40 // default constructor, destructor, copy constructor assignment operator and helpers
41 //////////
42
43 idx::idx() : inherited(TINFO_idx)
44 {
45         debugmsg("idx default constructor", LOGLEVEL_CONSTRUCT);
46 }
47
48 varidx::varidx() : covariant(false)
49 {
50         debugmsg("varidx default constructor", LOGLEVEL_CONSTRUCT);
51         tinfo_key = TINFO_varidx;
52 }
53
54 spinidx::spinidx() : dotted(false)
55 {
56         debugmsg("spinidx default constructor", LOGLEVEL_CONSTRUCT);
57         tinfo_key = TINFO_spinidx;
58 }
59
60 void idx::copy(const idx & other)
61 {
62         inherited::copy(other);
63         value = other.value;
64         dim = other.dim;
65 }
66
67 void varidx::copy(const varidx & other)
68 {
69         inherited::copy(other);
70         covariant = other.covariant;
71 }
72
73 void spinidx::copy(const spinidx & other)
74 {
75         inherited::copy(other);
76         dotted = other.dotted;
77 }
78
79 DEFAULT_DESTROY(idx)
80 DEFAULT_DESTROY(varidx)
81 DEFAULT_DESTROY(spinidx)
82
83 //////////
84 // other constructors
85 //////////
86
87 idx::idx(const ex & v, const ex & d) : inherited(TINFO_idx), value(v), dim(d)
88 {
89         debugmsg("idx constructor from ex,ex", LOGLEVEL_CONSTRUCT);
90         if (is_dim_numeric())
91                 if (!dim.info(info_flags::posint))
92                         throw(std::invalid_argument("dimension of space must be a positive integer"));
93 }
94
95 varidx::varidx(const ex & v, const ex & d, bool cov) : inherited(v, d), covariant(cov)
96 {
97         debugmsg("varidx constructor from ex,ex,bool", LOGLEVEL_CONSTRUCT);
98         tinfo_key = TINFO_varidx;
99 }
100
101 spinidx::spinidx(const ex & v, const ex & d, bool cov, bool dot) : inherited(v, d, cov), dotted(dot)
102 {
103         debugmsg("spinidx constructor from ex,ex,bool,bool", LOGLEVEL_CONSTRUCT);
104         tinfo_key = TINFO_spinidx;
105 }
106
107 //////////
108 // archiving
109 //////////
110
111 idx::idx(const archive_node &n, const lst &sym_lst) : inherited(n, sym_lst)
112 {
113         debugmsg("idx constructor from archive_node", LOGLEVEL_CONSTRUCT);
114         n.find_ex("value", value, sym_lst);
115         n.find_ex("dim", dim, sym_lst);
116 }
117
118 varidx::varidx(const archive_node &n, const lst &sym_lst) : inherited(n, sym_lst)
119 {
120         debugmsg("varidx constructor from archive_node", LOGLEVEL_CONSTRUCT);
121         n.find_bool("covariant", covariant);
122 }
123
124 spinidx::spinidx(const archive_node &n, const lst &sym_lst) : inherited(n, sym_lst)
125 {
126         debugmsg("spinidx constructor from archive_node", LOGLEVEL_CONSTRUCT);
127         n.find_bool("dotted", dotted);
128 }
129
130 void idx::archive(archive_node &n) const
131 {
132         inherited::archive(n);
133         n.add_ex("value", value);
134         n.add_ex("dim", dim);
135 }
136
137 void varidx::archive(archive_node &n) const
138 {
139         inherited::archive(n);
140         n.add_bool("covariant", covariant);
141 }
142
143 void spinidx::archive(archive_node &n) const
144 {
145         inherited::archive(n);
146         n.add_bool("dotted", dotted);
147 }
148
149 DEFAULT_UNARCHIVE(idx)
150 DEFAULT_UNARCHIVE(varidx)
151 DEFAULT_UNARCHIVE(spinidx)
152
153 //////////
154 // functions overriding virtual functions from bases classes
155 //////////
156
157 void idx::print(const print_context & c, unsigned level) const
158 {
159         debugmsg("idx print", LOGLEVEL_PRINT);
160
161         if (is_of_type(c, print_tree)) {
162
163                 c.s << std::string(level, ' ') << class_name()
164                     << std::hex << ", hash=0x" << hashvalue << ", flags=0x" << flags << std::dec
165                     << std::endl;
166                 unsigned delta_indent = static_cast<const print_tree &>(c).delta_indent;
167                 value.print(c, level + delta_indent);
168                 dim.print(c, level + delta_indent);
169
170         } else {
171
172                 if (!is_of_type(c, print_latex))
173                         c.s << ".";
174                 bool need_parens = !(is_ex_exactly_of_type(value, numeric) || is_ex_of_type(value, symbol));
175                 if (need_parens)
176                         c.s << "(";
177                 value.print(c);
178                 if (need_parens)
179                         c.s << ")";
180         }
181 }
182
183 void varidx::print(const print_context & c, unsigned level) const
184 {
185         debugmsg("varidx print", LOGLEVEL_PRINT);
186
187         if (is_of_type(c, print_tree)) {
188
189                 c.s << std::string(level, ' ') << class_name()
190                     << std::hex << ", hash=0x" << hashvalue << ", flags=0x" << flags << std::dec
191                     << (covariant ? ", covariant" : ", contravariant")
192                     << std::endl;
193                 unsigned delta_indent = static_cast<const print_tree &>(c).delta_indent;
194                 value.print(c, level + delta_indent);
195                 dim.print(c, level + delta_indent);
196
197         } else {
198
199                 if (!is_of_type(c, print_latex)) {
200                         if (covariant)
201                                 c.s << ".";
202                         else
203                                 c.s << "~";
204                 }
205                 bool need_parens = !(is_ex_exactly_of_type(value, numeric) || is_ex_of_type(value, symbol));
206                 if (need_parens)
207                         c.s << "(";
208                 value.print(c);
209                 if (need_parens)
210                         c.s << ")";
211         }
212 }
213
214 void spinidx::print(const print_context & c, unsigned level) const
215 {
216         debugmsg("spinidx print", LOGLEVEL_PRINT);
217
218         if (is_of_type(c, print_tree)) {
219
220                 c.s << std::string(level, ' ') << class_name()
221                     << std::hex << ", hash=0x" << hashvalue << ", flags=0x" << flags << std::dec
222                     << (covariant ? ", covariant" : ", contravariant")
223                     << (dotted ? ", dotted" : ", undotted")
224                     << std::endl;
225                 unsigned delta_indent = static_cast<const print_tree &>(c).delta_indent;
226                 value.print(c, level + delta_indent);
227                 dim.print(c, level + delta_indent);
228
229         } else {
230
231                 bool is_tex = is_of_type(c, print_latex);
232                 if (!is_tex) {
233                         if (covariant)
234                                 c.s << ".";
235                         else
236                                 c.s << "~";
237                 }
238                 if (dotted) {
239                         if (is_tex)
240                                 c.s << "\\dot{";
241                         else
242                                 c.s << "*";
243                 }
244                 bool need_parens = !(is_ex_exactly_of_type(value, numeric) || is_ex_of_type(value, symbol));
245                 if (need_parens)
246                         c.s << "(";
247                 value.print(c);
248                 if (need_parens)
249                         c.s << ")";
250                 if (is_tex && dotted)
251                         c.s << "}";
252         }
253 }
254
255 bool idx::info(unsigned inf) const
256 {
257         if (inf == info_flags::idx)
258                 return true;
259         return inherited::info(inf);
260 }
261
262 unsigned idx::nops() const
263 {
264         // don't count the dimension as that is not really a sub-expression
265         return 1;
266 }
267
268 ex & idx::let_op(int i)
269 {
270         GINAC_ASSERT(i == 0);
271         return value;
272 }
273
274 /** Returns order relation between two indices of the same type. The order
275  *  must be such that dummy indices lie next to each other. */
276 int idx::compare_same_type(const basic & other) const
277 {
278         GINAC_ASSERT(is_of_type(other, idx));
279         const idx &o = static_cast<const idx &>(other);
280
281         int cmpval = value.compare(o.value);
282         if (cmpval)
283                 return cmpval;
284         return dim.compare(o.dim);
285 }
286
287 int varidx::compare_same_type(const basic & other) const
288 {
289         GINAC_ASSERT(is_of_type(other, varidx));
290         const varidx &o = static_cast<const varidx &>(other);
291
292         int cmpval = inherited::compare_same_type(other);
293         if (cmpval)
294                 return cmpval;
295
296         // Check variance last so dummy indices will end up next to each other
297         if (covariant != o.covariant)
298                 return covariant ? -1 : 1;
299         return 0;
300 }
301
302 int spinidx::compare_same_type(const basic & other) const
303 {
304         GINAC_ASSERT(is_of_type(other, spinidx));
305         const spinidx &o = static_cast<const spinidx &>(other);
306
307         int cmpval = inherited::compare_same_type(other);
308         if (cmpval)
309                 return cmpval;
310
311         // Check variance and dottedness last so dummy indices will end up next to each other
312         if (covariant != o.covariant)
313                 return covariant ? -1 : 1;
314         if (dotted != o.dotted)
315                 return dotted ? -1 : 1;
316
317         return 0;
318 }
319
320 ex idx::subs(const lst & ls, const lst & lr) const
321 {
322         GINAC_ASSERT(ls.nops() == lr.nops());
323
324         // First look for index substitutions
325         for (unsigned i=0; i<ls.nops(); i++) {
326                 if (is_equal(*(ls.op(i)).bp)) {
327
328                         // Substitution index->index
329                         if (is_ex_of_type(lr.op(i), idx))
330                                 return lr.op(i);
331
332                         // Otherwise substitute value
333                         idx *i_copy = static_cast<idx *>(duplicate());
334                         i_copy->value = lr.op(i);
335                         i_copy->clearflag(status_flags::hash_calculated);
336                         return i_copy->setflag(status_flags::dynallocated);
337                 }
338         }
339
340         // None, substitute objects in value (not in dimension)
341         const ex &subsed_value = value.subs(ls, lr);
342         if (are_ex_trivially_equal(value, subsed_value))
343                 return *this;
344
345         idx *i_copy = static_cast<idx *>(duplicate());
346         i_copy->value = subsed_value;
347         i_copy->clearflag(status_flags::hash_calculated);
348         return i_copy->setflag(status_flags::dynallocated);
349 }
350
351 //////////
352 // new virtual functions
353 //////////
354
355 bool idx::is_dummy_pair_same_type(const basic & other) const
356 {
357         const idx &o = static_cast<const idx &>(other);
358
359         // Only pure symbols form dummy pairs, "2n+1" doesn't
360         if (!is_ex_of_type(value, symbol))
361                 return false;
362
363         // Value must be equal, of course
364         if (!value.is_equal(o.value))
365                 return false;
366
367         // Also the dimension
368         return dim.is_equal(o.dim);
369 }
370
371 bool varidx::is_dummy_pair_same_type(const basic & other) const
372 {
373         const varidx &o = static_cast<const varidx &>(other);
374
375         // Variance must be opposite
376         if (covariant == o.covariant)
377                 return false;
378
379         return inherited::is_dummy_pair_same_type(other);
380 }
381
382 bool spinidx::is_dummy_pair_same_type(const basic & other) const
383 {
384         const spinidx &o = static_cast<const spinidx &>(other);
385
386         // Dottedness must be the same
387         if (dotted != o.dotted)
388                 return false;
389
390         return inherited::is_dummy_pair_same_type(other);
391 }
392
393
394 //////////
395 // non-virtual functions
396 //////////
397
398 ex varidx::toggle_variance(void) const
399 {
400         varidx *i_copy = static_cast<varidx *>(duplicate());
401         i_copy->covariant = !i_copy->covariant;
402         i_copy->clearflag(status_flags::hash_calculated);
403         return i_copy->setflag(status_flags::dynallocated);
404 }
405
406 ex spinidx::toggle_dot(void) const
407 {
408         spinidx *i_copy = static_cast<spinidx *>(duplicate());
409         i_copy->dotted = !i_copy->dotted;
410         i_copy->clearflag(status_flags::hash_calculated);
411         return i_copy->setflag(status_flags::dynallocated);
412 }
413
414 ex spinidx::toggle_variance_dot(void) const
415 {
416         spinidx *i_copy = static_cast<spinidx *>(duplicate());
417         i_copy->covariant = !i_copy->covariant;
418         i_copy->dotted = !i_copy->dotted;
419         i_copy->clearflag(status_flags::hash_calculated);
420         return i_copy->setflag(status_flags::dynallocated);
421 }
422
423 //////////
424 // global functions
425 //////////
426
427 bool is_dummy_pair(const idx & i1, const idx & i2)
428 {
429         // The indices must be of exactly the same type
430         if (i1.tinfo() != i2.tinfo())
431                 return false;
432
433         // Same type, let the indices decide whether they are paired
434         return i1.is_dummy_pair_same_type(i2);
435 }
436
437 bool is_dummy_pair(const ex & e1, const ex & e2)
438 {
439         // The expressions must be indices
440         if (!is_ex_of_type(e1, idx) || !is_ex_of_type(e2, idx))
441                 return false;
442
443         return is_dummy_pair(ex_to_idx(e1), ex_to_idx(e2));
444 }
445
446 /** Bring a vector of indices into a canonic order. Dummy indices will lie
447  *  next to each other after the sorting. */
448 static void sort_index_vector(exvector &v)
449 {
450         // Nothing to sort if less than 2 elements
451         if (v.size() < 2)
452                 return;
453
454         // Simple bubble sort algorithm should be sufficient for the small
455         // number of indices expected
456         exvector::iterator it1 = v.begin(), itend = v.end(), next_to_last_idx = itend - 1;
457         while (it1 != next_to_last_idx) {
458                 exvector::iterator it2 = it1 + 1;
459                 while (it2 != itend) {
460                         if (it1->compare(*it2) > 0)
461                                 it1->swap(*it2);
462                         it2++;
463                 }
464                 it1++;
465         }
466 }
467
468
469 void find_free_and_dummy(exvector::const_iterator it, exvector::const_iterator itend, exvector & out_free, exvector & out_dummy)
470 {
471         out_free.clear();
472         out_dummy.clear();
473
474         // No indices? Then do nothing
475         if (it == itend)
476                 return;
477
478         // Only one index? Then it is a free one if it's not numeric
479         if (itend - it == 1) {
480                 if (ex_to_idx(*it).is_symbolic())
481                         out_free.push_back(*it);
482                 return;
483         }
484
485         // Sort index vector. This will cause dummy indices come to lie next
486         // to each other (because the sort order is defined to guarantee this).
487         exvector v(it, itend);
488         sort_index_vector(v);
489
490         // Find dummy pairs and free indices
491         it = v.begin(); itend = v.end();
492         exvector::const_iterator last = it++;
493         while (it != itend) {
494                 if (is_dummy_pair(*it, *last)) {
495                         out_dummy.push_back(*last);
496                         it++;
497                         if (it == itend)
498                                 return;
499                 } else {
500                         if (!it->is_equal(*last) && ex_to_idx(*last).is_symbolic())
501                                 out_free.push_back(*last);
502                 }
503                 last = it++;
504         }
505         if (ex_to_idx(*last).is_symbolic())
506                 out_free.push_back(*last);
507 }
508
509 exvector index_set_difference(const exvector & set1, const exvector & set2)
510 {
511         exvector ret;
512
513         exvector::const_iterator ait = set1.begin(), aitend = set1.end();
514         while (ait != aitend) {
515                 exvector::const_iterator bit = set2.begin(), bitend = set2.end();
516                 bool found = false;
517                 while (bit != bitend) {
518                         if (ait->is_equal(*bit)) {
519                                 found = true;
520                                 break;
521                         }
522                         bit++;
523                 }
524                 if (!found)
525                         ret.push_back(*ait);
526                 ait++;
527         }
528
529         return ret;
530 }
531
532 } // namespace GiNaC