Added complex conjugation methods and GiNaC function "conjugate".
[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 "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_latex inherited from idx
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(TINFO_idx) {}
57
58 varidx::varidx() : covariant(false)
59 {
60         tinfo_key = TINFO_varidx;
61 }
62
63 spinidx::spinidx() : dotted(false)
64 {
65         tinfo_key = TINFO_spinidx;
66 }
67
68 //////////
69 // other constructors
70 //////////
71
72 idx::idx(const ex & v, const ex & d) : inherited(TINFO_idx), 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 = TINFO_varidx;
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 = TINFO_spinidx;
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         idx *copy = duplicate();
246         copy->setflag(status_flags::dynallocated);
247         copy->clearflag(status_flags::hash_calculated);
248         copy->value = f(value);
249         return *copy;
250 }
251
252 /** Returns order relation between two indices of the same type. The order
253  *  must be such that dummy indices lie next to each other. */
254 int idx::compare_same_type(const basic & other) const
255 {
256         GINAC_ASSERT(is_a<idx>(other));
257         const idx &o = static_cast<const idx &>(other);
258
259         int cmpval = value.compare(o.value);
260         if (cmpval)
261                 return cmpval;
262         return dim.compare(o.dim);
263 }
264
265 bool idx::match_same_type(const basic & other) const
266 {
267         GINAC_ASSERT(is_a<idx>(other));
268         const idx &o = static_cast<const idx &>(other);
269
270         return dim.is_equal(o.dim);
271 }
272
273 int varidx::compare_same_type(const basic & other) const
274 {
275         GINAC_ASSERT(is_a<varidx>(other));
276         const varidx &o = static_cast<const varidx &>(other);
277
278         int cmpval = inherited::compare_same_type(other);
279         if (cmpval)
280                 return cmpval;
281
282         // Check variance last so dummy indices will end up next to each other
283         if (covariant != o.covariant)
284                 return covariant ? -1 : 1;
285
286         return 0;
287 }
288
289 bool varidx::match_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         if (covariant != o.covariant)
295                 return false;
296
297         return inherited::match_same_type(other);
298 }
299
300 int spinidx::compare_same_type(const basic & other) const
301 {
302         GINAC_ASSERT(is_a<spinidx>(other));
303         const spinidx &o = static_cast<const spinidx &>(other);
304
305         // Check dottedness first so dummy indices will end up next to each other
306         if (dotted != o.dotted)
307                 return dotted ? -1 : 1;
308
309         int cmpval = inherited::compare_same_type(other);
310         if (cmpval)
311                 return cmpval;
312
313         return 0;
314 }
315
316 bool spinidx::match_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         if (dotted != o.dotted)
322                 return false;
323         return inherited::match_same_type(other);
324 }
325
326 unsigned idx::calchash() const
327 {
328         // NOTE: The code in simplify_indexed() assumes that canonically
329         // ordered sequences of indices have the two members of dummy index
330         // pairs lying next to each other. The hash values for indices must
331         // be devised accordingly. The easiest (only?) way to guarantee the
332         // desired ordering is to make indices with the same value have equal
333         // hash keys. That is, the hash values must not depend on the index
334         // dimensions or other attributes (variance etc.).
335         // The compare_same_type() methods will take care of the rest.
336         unsigned v = golden_ratio_hash(tinfo());
337         v = rotate_left(v);
338         v ^= value.gethash();
339
340         // Store calculated hash value only if object is already evaluated
341         if (flags & status_flags::evaluated) {
342                 setflag(status_flags::hash_calculated);
343                 hashvalue = v;
344         }
345
346         return v;
347 }
348
349 /** By default, basic::evalf would evaluate the index value but we don't want
350  *  a.1 to become a.(1.0). */
351 ex idx::evalf(int level) const
352 {
353         return *this;
354 }
355
356 ex idx::subs(const exmap & m, unsigned options) const
357 {
358         // First look for index substitutions
359         exmap::const_iterator it = m.find(*this);
360         if (it != m.end()) {
361
362                 // Substitution index->index
363                 if (is_a<idx>(it->second))
364                         return it->second;
365
366                 // Otherwise substitute value
367                 idx *i_copy = duplicate();
368                 i_copy->value = it->second;
369                 i_copy->clearflag(status_flags::hash_calculated);
370                 return i_copy->setflag(status_flags::dynallocated);
371         }
372
373         // None, substitute objects in value (not in dimension)
374         const ex &subsed_value = value.subs(m, options);
375         if (are_ex_trivially_equal(value, subsed_value))
376                 return *this;
377
378         idx *i_copy = duplicate();
379         i_copy->value = subsed_value;
380         i_copy->clearflag(status_flags::hash_calculated);
381         return i_copy->setflag(status_flags::dynallocated);
382 }
383
384 /** Implementation of ex::diff() for an index always returns 0.
385  *
386  *  @see ex::diff */
387 ex idx::derivative(const symbol & s) const
388 {
389         return _ex0;
390 }
391
392 //////////
393 // new virtual functions
394 //////////
395
396 bool idx::is_dummy_pair_same_type(const basic & other) const
397 {
398         const idx &o = static_cast<const idx &>(other);
399
400         // Only pure symbols form dummy pairs, "2n+1" doesn'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 = 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() const
457 {
458         varidx *i_copy = 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() const
465 {
466         spinidx *i_copy = 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() const
473 {
474         spinidx *i_copy = 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