G_do_hoelder: fix case with real x values which are not of type cl_R.
[ginac.git] / ginac / idx.cpp
1 /** @file idx.cpp
2  *
3  *  Implementation of GiNaC's indices. */
4
5 /*
6  *  GiNaC Copyright (C) 1999-2018 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 "idx.h"
24 #include "symbol.h"
25 #include "lst.h"
26 #include "relational.h"
27 #include "operators.h"
28 #include "archive.h"
29 #include "utils.h"
30 #include "hash_seed.h"
31
32 #include <iostream>
33 #include <sstream>
34 #include <stdexcept>
35
36 namespace GiNaC {
37
38 GINAC_IMPLEMENT_REGISTERED_CLASS_OPT(idx, basic,
39   print_func<print_context>(&idx::do_print).
40   print_func<print_latex>(&idx::do_print_latex).
41   print_func<print_csrc>(&idx::do_print_csrc).
42   print_func<print_tree>(&idx::do_print_tree))
43
44 GINAC_IMPLEMENT_REGISTERED_CLASS_OPT(varidx, idx,
45   print_func<print_context>(&varidx::do_print).
46   print_func<print_latex>(&varidx::do_print_latex).
47   print_func<print_tree>(&varidx::do_print_tree))
48
49 GINAC_IMPLEMENT_REGISTERED_CLASS_OPT(spinidx, varidx,
50   print_func<print_context>(&spinidx::do_print).
51   print_func<print_latex>(&spinidx::do_print_latex).
52   print_func<print_tree>(&spinidx::do_print_tree))
53
54 //////////
55 // default constructor
56 //////////
57
58 idx::idx() {}
59
60 varidx::varidx() : covariant(false)
61 {
62 }
63
64 spinidx::spinidx() : dotted(false)
65 {
66 }
67
68 //////////
69 // other constructors
70 //////////
71
72 idx::idx(const ex & v, const ex & d) :  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 }
82
83 spinidx::spinidx(const ex & v, const ex & d, bool cov, bool dot) : inherited(v, d, cov), dotted(dot)
84 {
85 }
86
87 //////////
88 // archiving
89 //////////
90
91 void idx::read_archive(const archive_node& n, lst& sym_lst) 
92 {
93         inherited::read_archive(n, sym_lst);
94         n.find_ex("value", value, sym_lst);
95         n.find_ex("dim", dim, sym_lst);
96 }
97 GINAC_BIND_UNARCHIVER(idx);
98
99 void varidx::read_archive(const archive_node& n, lst& sym_lst)
100 {
101         inherited::read_archive(n, sym_lst);
102         n.find_bool("covariant", covariant);
103 }
104 GINAC_BIND_UNARCHIVER(varidx);
105
106 void spinidx::read_archive(const archive_node& n, lst& sym_lst)
107 {
108         inherited::read_archive(n, sym_lst);
109         n.find_bool("dotted", dotted);
110 }
111 GINAC_BIND_UNARCHIVER(spinidx);
112
113 void idx::archive(archive_node &n) const
114 {
115         inherited::archive(n);
116         n.add_ex("value", value);
117         n.add_ex("dim", dim);
118 }
119
120 void varidx::archive(archive_node &n) const
121 {
122         inherited::archive(n);
123         n.add_bool("covariant", covariant);
124 }
125
126 void spinidx::archive(archive_node &n) const
127 {
128         inherited::archive(n);
129         n.add_bool("dotted", dotted);
130 }
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_csrc(const print_csrc & c, unsigned level) const
165 {
166         c.s << "[";
167         if (value.info(info_flags::integer))
168                 c.s << ex_to<numeric>(value).to_int();
169         else
170                 value.print(c);
171         c.s << "]";
172 }
173
174 void idx::do_print_tree(const print_tree & c, unsigned level) const
175 {
176         c.s << std::string(level, ' ') << class_name() << " @" << this
177             << std::hex << ", hash=0x" << hashvalue << ", flags=0x" << flags << std::dec
178             << std::endl;
179         value.print(c, level +  c.delta_indent);
180         dim.print(c, level + c.delta_indent);
181 }
182
183 void varidx::do_print(const print_context & c, unsigned level) const
184 {
185         if (covariant)
186                 c.s << ".";
187         else
188                 c.s << "~";
189         print_index(c, level);
190 }
191
192 void varidx::do_print_tree(const print_tree & c, unsigned level) const
193 {
194         c.s << std::string(level, ' ') << class_name() << " @" << this
195             << std::hex << ", hash=0x" << hashvalue << ", flags=0x" << flags << std::dec
196             << (covariant ? ", covariant" : ", contravariant")
197             << std::endl;
198         value.print(c, level + c.delta_indent);
199         dim.print(c, level + c.delta_indent);
200 }
201
202 void spinidx::do_print(const print_context & c, unsigned level) const
203 {
204         if (covariant)
205                 c.s << ".";
206         else
207                 c.s << "~";
208         if (dotted)
209                 c.s << "*";
210         print_index(c, level);
211 }
212
213 void spinidx::do_print_latex(const print_latex & c, unsigned level) const
214 {
215         if (dotted)
216                 c.s << "\\dot{";
217         else
218                 c.s << "{";
219         print_index(c, level);
220         c.s << "}";
221 }
222
223 void spinidx::do_print_tree(const print_tree & c, unsigned level) const
224 {
225         c.s << std::string(level, ' ') << class_name() << " @" << this
226             << std::hex << ", hash=0x" << hashvalue << ", flags=0x" << flags << std::dec
227             << (covariant ? ", covariant" : ", contravariant")
228             << (dotted ? ", dotted" : ", undotted")
229             << std::endl;
230         value.print(c, level + c.delta_indent);
231         dim.print(c, level + c.delta_indent);
232 }
233
234 bool idx::info(unsigned inf) const
235 {
236         switch(inf) {
237                 case info_flags::idx:
238                 case info_flags::has_indices:
239                         return true;
240         }
241         return inherited::info(inf);
242 }
243
244 size_t idx::nops() const
245 {
246         // don't count the dimension as that is not really a sub-expression
247         return 1;
248 }
249
250 ex idx::op(size_t i) const
251 {
252         GINAC_ASSERT(i == 0);
253         return value;
254 }
255
256 ex idx::map(map_function & f) const
257 {
258         const ex &mapped_value = f(value);
259         if (are_ex_trivially_equal(value, mapped_value))
260                 return *this;
261         else {
262                 idx *copy = duplicate();
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         unsigned v = make_hash_seed(typeid(*this));
354         v = rotate_left(v);
355         v ^= value.gethash();
356
357         // Store calculated hash value only if object is already evaluated
358         if (flags & status_flags::evaluated) {
359                 setflag(status_flags::hash_calculated);
360                 hashvalue = v;
361         }
362
363         return v;
364 }
365
366 /** By default, basic::evalf would evaluate the index value but we don't want
367  *  a.1 to become a.(1.0). */
368 ex idx::evalf() const
369 {
370         return *this;
371 }
372
373 ex idx::subs(const exmap & m, unsigned options) const
374 {
375         // First look for index substitutions
376         auto it = m.find(*this);
377         if (it != m.end()) {
378
379                 // Substitution index->index
380                 if (is_a<idx>(it->second) || (options & subs_options::really_subs_idx))
381                         return it->second;
382
383                 // Otherwise substitute value
384                 idx *i_copy = duplicate();
385                 i_copy->value = it->second;
386                 i_copy->clearflag(status_flags::hash_calculated);
387                 return *i_copy;
388         }
389
390         // None, substitute objects in value (not in dimension)
391         const ex &subsed_value = value.subs(m, options);
392         if (are_ex_trivially_equal(value, subsed_value))
393                 return *this;
394
395         idx *i_copy = duplicate();
396         i_copy->value = subsed_value;
397         i_copy->clearflag(status_flags::hash_calculated);
398         return *i_copy;
399 }
400
401 /** Implementation of ex::diff() for an index always returns 0.
402  *
403  *  @see ex::diff */
404 ex idx::derivative(const symbol & s) const
405 {
406         return _ex0;
407 }
408
409 //////////
410 // new virtual functions
411 //////////
412
413 bool idx::is_dummy_pair_same_type(const basic & other) const
414 {
415         const idx &o = static_cast<const idx &>(other);
416
417         // Only pure symbols form dummy pairs, "2n+1" doesn't
418         if (!is_a<symbol>(value))
419                 return false;
420
421         // Value must be equal, of course
422         if (!value.is_equal(o.value))
423                 return false;
424
425         // Dimensions need not be equal but must be comparable (so we can
426         // determine the minimum dimension of contractions)
427         if (dim.is_equal(o.dim))
428                 return true;
429
430         return is_exactly_a<numeric>(dim) || is_exactly_a<numeric>(o.dim);
431 }
432
433 bool varidx::is_dummy_pair_same_type(const basic & other) const
434 {
435         const varidx &o = static_cast<const varidx &>(other);
436
437         // Variance must be opposite
438         if (covariant == o.covariant)
439                 return false;
440
441         return inherited::is_dummy_pair_same_type(other);
442 }
443
444 bool spinidx::is_dummy_pair_same_type(const basic & other) const
445 {
446         const spinidx &o = static_cast<const spinidx &>(other);
447
448         // Dottedness must be the same
449         if (dotted != o.dotted)
450                 return false;
451
452         return inherited::is_dummy_pair_same_type(other);
453 }
454
455
456 //////////
457 // non-virtual functions
458 //////////
459
460 ex idx::replace_dim(const ex & new_dim) const
461 {
462         idx *i_copy = duplicate();
463         i_copy->dim = new_dim;
464         i_copy->clearflag(status_flags::hash_calculated);
465         return *i_copy;
466 }
467
468 ex idx::minimal_dim(const idx & other) const
469 {
470         return GiNaC::minimal_dim(dim, other.dim);
471 }
472
473 ex varidx::toggle_variance() const
474 {
475         varidx *i_copy = duplicate();
476         i_copy->covariant = !i_copy->covariant;
477         i_copy->clearflag(status_flags::hash_calculated);
478         return *i_copy;
479 }
480
481 ex spinidx::toggle_dot() const
482 {
483         spinidx *i_copy = duplicate();
484         i_copy->dotted = !i_copy->dotted;
485         i_copy->clearflag(status_flags::hash_calculated);
486         return *i_copy;
487 }
488
489 ex spinidx::toggle_variance_dot() const
490 {
491         spinidx *i_copy = duplicate();
492         i_copy->covariant = !i_copy->covariant;
493         i_copy->dotted = !i_copy->dotted;
494         i_copy->clearflag(status_flags::hash_calculated);
495         return *i_copy;
496 }
497
498 //////////
499 // global functions
500 //////////
501
502 bool is_dummy_pair(const idx & i1, const idx & i2)
503 {
504         // The indices must be of exactly the same type
505         if (typeid(i1) != typeid(i2))
506                 return false;
507
508         // Same type, let the indices decide whether they are paired
509         return i1.is_dummy_pair_same_type(i2);
510 }
511
512 bool is_dummy_pair(const ex & e1, const ex & e2)
513 {
514         // The expressions must be indices
515         if (!is_a<idx>(e1) || !is_a<idx>(e2))
516                 return false;
517
518         return is_dummy_pair(ex_to<idx>(e1), ex_to<idx>(e2));
519 }
520
521 void find_free_and_dummy(exvector::const_iterator it, exvector::const_iterator itend, exvector & out_free, exvector & out_dummy)
522 {
523         out_free.clear();
524         out_dummy.clear();
525
526         // No indices? Then do nothing
527         if (it == itend)
528                 return;
529
530         // Only one index? Then it is a free one if it's not numeric
531         if (itend - it == 1) {
532                 if (ex_to<idx>(*it).is_symbolic())
533                         out_free.push_back(*it);
534                 return;
535         }
536
537         // Sort index vector. This will cause dummy indices come to lie next
538         // to each other (because the sort order is defined to guarantee this).
539         exvector v(it, itend);
540         shaker_sort(v.begin(), v.end(), ex_is_less(), ex_swap());
541
542         // Find dummy pairs and free indices
543         it = v.begin(); itend = v.end();
544         auto last = it++;
545         while (it != itend) {
546                 if (is_dummy_pair(*it, *last)) {
547                         out_dummy.push_back(*last);
548                         it++;
549                         if (it == itend)
550                                 return;
551                 } else {
552                         if (!it->is_equal(*last) && ex_to<idx>(*last).is_symbolic())
553                                 out_free.push_back(*last);
554                 }
555                 last = it++;
556         }
557         if (ex_to<idx>(*last).is_symbolic())
558                 out_free.push_back(*last);
559 }
560
561 ex minimal_dim(const ex & dim1, const ex & dim2)
562 {
563         if (dim1.is_equal(dim2) || dim1 < dim2 || (is_exactly_a<numeric>(dim1) && !is_a<numeric>(dim2)))
564                 return dim1;
565         else if (dim1 > dim2 || (!is_a<numeric>(dim1) && is_exactly_a<numeric>(dim2)))
566                 return dim2;
567         else {
568                 std::ostringstream s;
569                 s << "minimal_dim(): index dimensions " << dim1 << " and " << dim2 << " cannot be ordered";
570                 throw (std::runtime_error(s.str()));
571         }
572 }
573
574 } // namespace GiNaC