3 * Implementation of GiNaC's indices. */
6 * GiNaC Copyright (C) 1999-2008 Johannes Gutenberg University Mainz, Germany
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.
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.
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
30 #include "relational.h"
31 #include "operators.h"
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))
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))
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))
54 // default constructor
59 varidx::varidx() : covariant(false)
63 spinidx::spinidx() : dotted(false)
71 idx::idx(const ex & v, const ex & d) : value(v), dim(d)
74 if (!dim.info(info_flags::posint))
75 throw(std::invalid_argument("dimension of space must be a positive integer"));
78 varidx::varidx(const ex & v, const ex & d, bool cov) : inherited(v, d), covariant(cov)
82 spinidx::spinidx(const ex & v, const ex & d, bool cov, bool dot) : inherited(v, d, cov), dotted(dot)
90 void idx::read_archive(const archive_node& n, lst& sym_lst)
92 inherited::read_archive(n, sym_lst);
93 n.find_ex("value", value, sym_lst);
94 n.find_ex("dim", dim, sym_lst);
96 GINAC_BIND_UNARCHIVER(idx);
98 void varidx::read_archive(const archive_node& n, lst& sym_lst)
100 inherited::read_archive(n, sym_lst);
101 n.find_bool("covariant", covariant);
103 GINAC_BIND_UNARCHIVER(varidx);
105 void spinidx::read_archive(const archive_node& n, lst& sym_lst)
107 inherited::read_archive(n, sym_lst);
108 n.find_bool("dotted", dotted);
110 GINAC_BIND_UNARCHIVER(spinidx);
112 void idx::archive(archive_node &n) const
114 inherited::archive(n);
115 n.add_ex("value", value);
116 n.add_ex("dim", dim);
119 void varidx::archive(archive_node &n) const
121 inherited::archive(n);
122 n.add_bool("covariant", covariant);
125 void spinidx::archive(archive_node &n) const
127 inherited::archive(n);
128 n.add_bool("dotted", dotted);
132 // functions overriding virtual functions from base classes
135 void idx::print_index(const print_context & c, unsigned level) const
137 bool need_parens = !(is_exactly_a<numeric>(value) || is_a<symbol>(value));
143 if (c.options & print_options::print_index_dimensions) {
150 void idx::do_print(const print_context & c, unsigned level) const
153 print_index(c, level);
156 void idx::do_print_latex(const print_latex & c, unsigned level) const
159 print_index(c, level);
163 void idx::do_print_csrc(const print_csrc & c, unsigned level) const
166 if (value.info(info_flags::integer))
167 c.s << ex_to<numeric>(value).to_int();
173 void idx::do_print_tree(const print_tree & c, unsigned level) const
175 c.s << std::string(level, ' ') << class_name() << " @" << this
176 << std::hex << ", hash=0x" << hashvalue << ", flags=0x" << flags << std::dec
178 value.print(c, level + c.delta_indent);
179 dim.print(c, level + c.delta_indent);
182 void varidx::do_print(const print_context & c, unsigned level) const
188 print_index(c, level);
191 void varidx::do_print_tree(const print_tree & c, unsigned level) const
193 c.s << std::string(level, ' ') << class_name() << " @" << this
194 << std::hex << ", hash=0x" << hashvalue << ", flags=0x" << flags << std::dec
195 << (covariant ? ", covariant" : ", contravariant")
197 value.print(c, level + c.delta_indent);
198 dim.print(c, level + c.delta_indent);
201 void spinidx::do_print(const print_context & c, unsigned level) const
209 print_index(c, level);
212 void spinidx::do_print_latex(const print_latex & c, unsigned level) const
218 print_index(c, level);
222 void spinidx::do_print_tree(const print_tree & c, unsigned level) const
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")
229 value.print(c, level + c.delta_indent);
230 dim.print(c, level + c.delta_indent);
233 bool idx::info(unsigned inf) const
236 case info_flags::idx:
237 case info_flags::has_indices:
240 return inherited::info(inf);
243 size_t idx::nops() const
245 // don't count the dimension as that is not really a sub-expression
249 ex idx::op(size_t i) const
251 GINAC_ASSERT(i == 0);
255 ex idx::map(map_function & f) const
257 const ex &mapped_value = f(value);
258 if (are_ex_trivially_equal(value, mapped_value))
261 idx *copy = duplicate();
262 copy->setflag(status_flags::dynallocated);
263 copy->clearflag(status_flags::hash_calculated);
264 copy->value = mapped_value;
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
273 GINAC_ASSERT(is_a<idx>(other));
274 const idx &o = static_cast<const idx &>(other);
276 int cmpval = value.compare(o.value);
279 return dim.compare(o.dim);
282 bool idx::match_same_type(const basic & other) const
284 GINAC_ASSERT(is_a<idx>(other));
285 const idx &o = static_cast<const idx &>(other);
287 return dim.is_equal(o.dim);
290 int varidx::compare_same_type(const basic & other) const
292 GINAC_ASSERT(is_a<varidx>(other));
293 const varidx &o = static_cast<const varidx &>(other);
295 int cmpval = inherited::compare_same_type(other);
299 // Check variance last so dummy indices will end up next to each other
300 if (covariant != o.covariant)
301 return covariant ? -1 : 1;
306 bool varidx::match_same_type(const basic & other) const
308 GINAC_ASSERT(is_a<varidx>(other));
309 const varidx &o = static_cast<const varidx &>(other);
311 if (covariant != o.covariant)
314 return inherited::match_same_type(other);
317 int spinidx::compare_same_type(const basic & other) const
319 GINAC_ASSERT(is_a<spinidx>(other));
320 const spinidx &o = static_cast<const spinidx &>(other);
322 // Check dottedness first so dummy indices will end up next to each other
323 if (dotted != o.dotted)
324 return dotted ? -1 : 1;
326 int cmpval = inherited::compare_same_type(other);
333 bool spinidx::match_same_type(const basic & other) const
335 GINAC_ASSERT(is_a<spinidx>(other));
336 const spinidx &o = static_cast<const spinidx &>(other);
338 if (dotted != o.dotted)
340 return inherited::match_same_type(other);
343 unsigned idx::calchash() const
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);
356 v ^= value.gethash();
358 // Store calculated hash value only if object is already evaluated
359 if (flags & status_flags::evaluated) {
360 setflag(status_flags::hash_calculated);
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
374 ex idx::subs(const exmap & m, unsigned options) const
376 // First look for index substitutions
377 exmap::const_iterator it = m.find(*this);
380 // Substitution index->index
381 if (is_a<idx>(it->second) || (options & subs_options::really_subs_idx))
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);
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))
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);
402 /** Implementation of ex::diff() for an index always returns 0.
405 ex idx::derivative(const symbol & s) const
411 // new virtual functions
414 bool idx::is_dummy_pair_same_type(const basic & other) const
416 const idx &o = static_cast<const idx &>(other);
418 // Only pure symbols form dummy pairs, "2n+1" doesn't
419 if (!is_a<symbol>(value))
422 // Value must be equal, of course
423 if (!value.is_equal(o.value))
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))
431 return is_exactly_a<numeric>(dim) || is_exactly_a<numeric>(o.dim);
434 bool varidx::is_dummy_pair_same_type(const basic & other) const
436 const varidx &o = static_cast<const varidx &>(other);
438 // Variance must be opposite
439 if (covariant == o.covariant)
442 return inherited::is_dummy_pair_same_type(other);
445 bool spinidx::is_dummy_pair_same_type(const basic & other) const
447 const spinidx &o = static_cast<const spinidx &>(other);
449 // Dottedness must be the same
450 if (dotted != o.dotted)
453 return inherited::is_dummy_pair_same_type(other);
458 // non-virtual functions
461 ex idx::replace_dim(const ex & new_dim) const
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);
469 ex idx::minimal_dim(const idx & other) const
471 return GiNaC::minimal_dim(dim, other.dim);
474 ex varidx::toggle_variance() const
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);
482 ex spinidx::toggle_dot() const
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);
490 ex spinidx::toggle_variance_dot() const
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);
503 bool is_dummy_pair(const idx & i1, const idx & i2)
505 // The indices must be of exactly the same type
506 if (typeid(i1) != typeid(i2))
509 // Same type, let the indices decide whether they are paired
510 return i1.is_dummy_pair_same_type(i2);
513 bool is_dummy_pair(const ex & e1, const ex & e2)
515 // The expressions must be indices
516 if (!is_a<idx>(e1) || !is_a<idx>(e2))
519 return is_dummy_pair(ex_to<idx>(e1), ex_to<idx>(e2));
522 void find_free_and_dummy(exvector::const_iterator it, exvector::const_iterator itend, exvector & out_free, exvector & out_dummy)
527 // No indices? Then do nothing
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);
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());
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);
553 if (!it->is_equal(*last) && ex_to<idx>(*last).is_symbolic())
554 out_free.push_back(*last);
558 if (ex_to<idx>(*last).is_symbolic())
559 out_free.push_back(*last);
562 ex minimal_dim(const ex & dim1, const ex & dim2)
564 if (dim1.is_equal(dim2) || dim1 < dim2 || (is_exactly_a<numeric>(dim1) && !is_a<numeric>(dim2)))
566 else if (dim1 > dim2 || (!is_a<numeric>(dim1) && is_exactly_a<numeric>(dim2)))
569 std::ostringstream s;
570 s << "minimal_dim(): index dimensions " << dim1 << " and " << dim2 << " cannot be ordered";
571 throw (std::runtime_error(s.str()));