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