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