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