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