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