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