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