44427a6ac08e3e7db642bf2d1381e3491ba1612a
[ginac.git] / ginac / idx.cpp
1 /** @file idx.cpp
2  *
3  *  Implementation of GiNaC's indices. */
4
5 /*
6  *  GiNaC Copyright (C) 1999-2001 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 <stdexcept>
24
25 #include "idx.h"
26 #include "symbol.h"
27 #include "lst.h"
28 #include "print.h"
29 #include "archive.h"
30 #include "utils.h"
31 #include "debugmsg.h"
32
33 namespace GiNaC {
34
35 GINAC_IMPLEMENT_REGISTERED_CLASS(idx, basic)
36 GINAC_IMPLEMENT_REGISTERED_CLASS(varidx, idx)
37
38 //////////
39 // default constructor, destructor, copy constructor assignment operator and helpers
40 //////////
41
42 idx::idx() : inherited(TINFO_idx)
43 {
44         debugmsg("idx default constructor", LOGLEVEL_CONSTRUCT);
45 }
46
47 varidx::varidx() : covariant(false)
48 {
49         debugmsg("varidx default constructor", LOGLEVEL_CONSTRUCT);
50         tinfo_key = TINFO_varidx;
51 }
52
53 void idx::copy(const idx & other)
54 {
55         inherited::copy(other);
56         value = other.value;
57         dim = other.dim;
58 }
59
60 void varidx::copy(const varidx & other)
61 {
62         inherited::copy(other);
63         covariant = other.covariant;
64 }
65
66 DEFAULT_DESTROY(idx)
67 DEFAULT_DESTROY(varidx)
68
69 //////////
70 // other constructors
71 //////////
72
73 idx::idx(const ex & v, const ex & d) : inherited(TINFO_idx), value(v), dim(d)
74 {
75         debugmsg("idx constructor from ex,ex", LOGLEVEL_CONSTRUCT);
76         if (is_dim_numeric())
77                 if (!dim.info(info_flags::posint))
78                         throw(std::invalid_argument("dimension of space must be a positive integer"));
79 }
80
81 varidx::varidx(const ex & v, const ex & d, bool cov) : inherited(v, d), covariant(cov)
82 {
83         debugmsg("varidx constructor from ex,ex,bool", LOGLEVEL_CONSTRUCT);
84         tinfo_key = TINFO_varidx;
85 }
86
87 //////////
88 // archiving
89 //////////
90
91 idx::idx(const archive_node &n, const lst &sym_lst) : inherited(n, sym_lst)
92 {
93         debugmsg("idx constructor from archive_node", LOGLEVEL_CONSTRUCT);
94         n.find_ex("value", value, sym_lst);
95         n.find_ex("dim", dim, sym_lst);
96 }
97
98 varidx::varidx(const archive_node &n, const lst &sym_lst) : inherited(n, sym_lst)
99 {
100         debugmsg("varidx constructor from archive_node", LOGLEVEL_CONSTRUCT);
101         n.find_bool("covariant", covariant);
102 }
103
104 void idx::archive(archive_node &n) const
105 {
106         inherited::archive(n);
107         n.add_ex("value", value);
108         n.add_ex("dim", dim);
109 }
110
111 void varidx::archive(archive_node &n) const
112 {
113         inherited::archive(n);
114         n.add_bool("covariant", covariant);
115 }
116
117 DEFAULT_UNARCHIVE(idx)
118 DEFAULT_UNARCHIVE(varidx)
119
120 //////////
121 // functions overriding virtual functions from bases classes
122 //////////
123
124 void idx::print(const print_context & c, unsigned level) const
125 {
126         debugmsg("idx print", LOGLEVEL_PRINT);
127
128         if (is_of_type(c, print_tree)) {
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_of_type(c, print_latex))
140                         c.s << ".";
141                 bool need_parens = !(is_ex_exactly_of_type(value, numeric) || is_ex_of_type(value, symbol));
142                 if (need_parens)
143                         c.s << "(";
144                 value.print(c);
145                 if (need_parens)
146                         c.s << ")";
147         }
148 }
149
150 void varidx::print(const print_context & c, unsigned level) const
151 {
152         debugmsg("varidx print", LOGLEVEL_PRINT);
153
154         if (is_of_type(c, print_tree)) {
155
156                 c.s << std::string(level, ' ') << class_name()
157                     << std::hex << ", hash=0x" << hashvalue << ", flags=0x" << flags << std::dec
158                     << (covariant ? ", covariant" : ", contravariant")
159                     << std::endl;
160                 unsigned delta_indent = static_cast<const print_tree &>(c).delta_indent;
161                 value.print(c, level + delta_indent);
162                 dim.print(c, level + delta_indent);
163
164         } else {
165
166                 if (!is_of_type(c, print_latex)) {
167                         if (covariant)
168                                 c.s << ".";
169                         else
170                                 c.s << "~";
171                 }
172                 bool need_parens = !(is_ex_exactly_of_type(value, numeric) || is_ex_of_type(value, symbol));
173                 if (need_parens)
174                         c.s << "(";
175                 value.print(c);
176                 if (need_parens)
177                         c.s << ")";
178         }
179 }
180
181 bool idx::info(unsigned inf) const
182 {
183         if (inf == info_flags::idx)
184                 return true;
185         return inherited::info(inf);
186 }
187
188 unsigned idx::nops() const
189 {
190         // don't count the dimension as that is not really a sub-expression
191         return 1;
192 }
193
194 ex & idx::let_op(int i)
195 {
196         GINAC_ASSERT(i == 0);
197         return value;
198 }
199
200 /** Returns order relation between two indices of the same type. The order
201  *  must be such that dummy indices lie next to each other. */
202 int idx::compare_same_type(const basic & other) const
203 {
204         GINAC_ASSERT(is_of_type(other, idx));
205         const idx &o = static_cast<const idx &>(other);
206
207         int cmpval = value.compare(o.value);
208         if (cmpval)
209                 return cmpval;
210         return dim.compare(o.dim);
211 }
212
213 int varidx::compare_same_type(const basic & other) const
214 {
215         GINAC_ASSERT(is_of_type(other, varidx));
216         const varidx &o = static_cast<const varidx &>(other);
217
218         int cmpval = inherited::compare_same_type(other);
219         if (cmpval)
220                 return cmpval;
221
222         // Check variance last so dummy indices will end up next to each other
223         if (covariant != o.covariant)
224                 return covariant ? -1 : 1;
225         return 0;
226 }
227
228 ex idx::subs(const lst & ls, const lst & lr) const
229 {
230         GINAC_ASSERT(ls.nops() == lr.nops());
231
232         // First look for index substitutions
233         for (unsigned i=0; i<ls.nops(); i++) {
234                 if (is_equal(*(ls.op(i)).bp)) {
235
236                         // Substitution index->index
237                         if (is_ex_of_type(lr.op(i), idx))
238                                 return lr.op(i);
239
240                         // Otherwise substitute value
241                         idx *i_copy = static_cast<idx *>(duplicate());
242                         i_copy->value = lr.op(i);
243                         i_copy->clearflag(status_flags::hash_calculated);
244                         return i_copy->setflag(status_flags::dynallocated);
245                 }
246         }
247
248         // None, substitute objects in value (not in dimension)
249         const ex &subsed_value = value.subs(ls, lr);
250         if (are_ex_trivially_equal(value, subsed_value))
251                 return *this;
252
253         idx *i_copy = static_cast<idx *>(duplicate());
254         i_copy->value = subsed_value;
255         i_copy->clearflag(status_flags::hash_calculated);
256         return i_copy->setflag(status_flags::dynallocated);
257 }
258
259 //////////
260 // new virtual functions
261 //////////
262
263 bool idx::is_dummy_pair_same_type(const basic & other) const
264 {
265         const idx &o = static_cast<const idx &>(other);
266
267         // Only pure symbols form dummy pairs, "2n+1" doesn't
268         if (!is_ex_of_type(value, symbol))
269                 return false;
270
271         // Value must be equal, of course
272         if (!value.is_equal(o.value))
273                 return false;
274
275         // Also the dimension
276         return dim.is_equal(o.dim);
277 }
278
279 bool varidx::is_dummy_pair_same_type(const basic & other) const
280 {
281         const varidx &o = static_cast<const varidx &>(other);
282
283         // Variance must be opposite
284         if (covariant == o.covariant)
285                 return false;
286
287         return inherited::is_dummy_pair_same_type(other);
288 }
289
290 //////////
291 // non-virtual functions
292 //////////
293
294 ex varidx::toggle_variance(void) const
295 {
296         varidx *i_copy = static_cast<varidx *>(duplicate());
297         i_copy->covariant = !i_copy->covariant;
298         i_copy->clearflag(status_flags::hash_calculated);
299         return i_copy->setflag(status_flags::dynallocated);
300 }
301
302 //////////
303 // global functions
304 //////////
305
306 bool is_dummy_pair(const idx & i1, const idx & i2)
307 {
308         // The indices must be of exactly the same type
309         if (i1.tinfo() != i2.tinfo())
310                 return false;
311
312         // Same type, let the indices decide whether they are paired
313         return i1.is_dummy_pair_same_type(i2);
314 }
315
316 bool is_dummy_pair(const ex & e1, const ex & e2)
317 {
318         // The expressions must be indices
319         if (!is_ex_of_type(e1, idx) || !is_ex_of_type(e2, idx))
320                 return false;
321
322         return is_dummy_pair(ex_to_idx(e1), ex_to_idx(e2));
323 }
324
325 /** Bring a vector of indices into a canonic order. Dummy indices will lie
326  *  next to each other after the sorting. */
327 static void sort_index_vector(exvector &v)
328 {
329         // Nothing to sort if less than 2 elements
330         if (v.size() < 2)
331                 return;
332
333         // Simple bubble sort algorithm should be sufficient for the small
334         // number of indices expected
335         exvector::iterator it1 = v.begin(), itend = v.end(), next_to_last_idx = itend - 1;
336         while (it1 != next_to_last_idx) {
337                 exvector::iterator it2 = it1 + 1;
338                 while (it2 != itend) {
339                         if (it1->compare(*it2) > 0)
340                                 it1->swap(*it2);
341                         it2++;
342                 }
343                 it1++;
344         }
345 }
346
347
348 void find_free_and_dummy(exvector::const_iterator it, exvector::const_iterator itend, exvector & out_free, exvector & out_dummy)
349 {
350         out_free.clear();
351         out_dummy.clear();
352
353         // No indices? Then do nothing
354         if (it == itend)
355                 return;
356
357         // Only one index? Then it is a free one if it's not numeric
358         if (itend - it == 1) {
359                 if (ex_to_idx(*it).is_symbolic())
360                         out_free.push_back(*it);
361                 return;
362         }
363
364         // Sort index vector. This will cause dummy indices come to lie next
365         // to each other (because the sort order is defined to guarantee this).
366         exvector v(it, itend);
367         sort_index_vector(v);
368
369         // Find dummy pairs and free indices
370         it = v.begin(); itend = v.end();
371         exvector::const_iterator last = it++;
372         while (it != itend) {
373                 if (is_dummy_pair(*it, *last)) {
374                         out_dummy.push_back(*last);
375                         it++;
376                         if (it == itend)
377                                 return;
378                 } else {
379                         if (!it->is_equal(*last) && ex_to_idx(*last).is_symbolic())
380                                 out_free.push_back(*last);
381                 }
382                 last = it++;
383         }
384         if (ex_to_idx(*last).is_symbolic())
385                 out_free.push_back(*last);
386 }
387
388 exvector index_set_difference(const exvector & set1, const exvector & set2)
389 {
390         exvector ret;
391
392         exvector::const_iterator ait = set1.begin(), aitend = set1.end();
393         while (ait != aitend) {
394                 exvector::const_iterator bit = set2.begin(), bitend = set2.end();
395                 bool found = false;
396                 while (bit != bitend) {
397                         if (ait->is_equal(*bit)) {
398                                 found = true;
399                                 break;
400                         }
401                         bit++;
402                 }
403                 if (!found)
404                         ret.push_back(*ait);
405                 ait++;
406         }
407
408         return ret;
409 }
410
411 } // namespace GiNaC