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