#ifndef around namespace GiNaC { }
[ginac.git] / ginac / idx.cpp
1 /** @file idx.cpp
2  *
3  *  Implementation of GiNaC's indices. */
4
5 /*
6  *  GiNaC Copyright (C) 1999 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 "ex.h"
27 #include "lst.h"
28 #include "relational.h"
29 #include "utils.h"
30 #include "debugmsg.h"
31
32 #ifndef NO_GINAC_NAMESPACE
33 namespace GiNaC {
34 #endif // ndef NO_GINAC_NAMESPACE
35
36 //////////
37 // default constructor, destructor, copy constructor assignment operator and helpers
38 //////////
39
40 // public
41
42 idx::idx() : basic(TINFO_idx), symbolic(true), covariant(false)
43 {
44     debugmsg("idx default constructor",LOGLEVEL_CONSTRUCT);
45     serial=next_serial++;
46     name="index"+ToString(serial);
47 }
48
49 idx::~idx() 
50 {
51     debugmsg("idx destructor",LOGLEVEL_DESTRUCT);
52     destroy(0);
53 }
54
55 idx::idx(idx const & other)
56 {
57     debugmsg("idx copy constructor",LOGLEVEL_CONSTRUCT);
58     copy(other);
59 }
60
61 idx const & idx::operator=(idx const & other)
62 {
63     debugmsg("idx operator=",LOGLEVEL_ASSIGNMENT);
64     if (this != &other) {
65         destroy(1);
66         copy(other);
67     }
68     return *this;
69 }
70
71 // protected
72
73 void idx::copy(idx const & other)
74 {
75     basic::copy(other);
76     serial=other.serial;
77     symbolic=other.symbolic;
78     name=other.name;
79     value=other.value;
80     covariant=other.covariant;
81 }
82
83 void idx::destroy(bool call_parent)
84 {
85     if (call_parent) basic::destroy(call_parent);
86 }
87
88 //////////
89 // other constructors
90 //////////
91
92 // public
93
94 idx::idx(bool cov) : basic(TINFO_idx), symbolic(true), covariant(cov)
95 {
96     debugmsg("idx constructor from bool",LOGLEVEL_CONSTRUCT);
97     serial=next_serial++;
98     name="index"+ToString(serial);
99 }
100
101 idx::idx(string const & n, bool cov) : basic(TINFO_idx),  
102     symbolic(true), name(n), covariant(cov)
103 {
104     debugmsg("idx constructor from string,bool",LOGLEVEL_CONSTRUCT);
105     serial=next_serial++;
106 }
107
108 idx::idx(char const * n, bool cov) : basic(TINFO_idx),  
109     symbolic(true), name(n), covariant(cov)
110 {
111     debugmsg("idx constructor from char*,bool",LOGLEVEL_CONSTRUCT);
112     serial=next_serial++;
113 }
114
115 idx::idx(unsigned const v, bool cov) : basic(TINFO_idx),
116     symbolic(false), value(v), covariant(cov)
117 {
118     debugmsg("idx constructor from unsigned,bool",LOGLEVEL_CONSTRUCT);
119     serial=0;
120 }
121
122
123 //////////
124 // functions overriding virtual functions from bases classes
125 //////////
126
127 // public
128
129 basic * idx::duplicate() const
130 {
131     debugmsg("idx duplicate",LOGLEVEL_DUPLICATE);
132     return new idx(*this);
133 }
134
135 void idx::printraw(ostream & os) const
136 {
137     debugmsg("idx printraw",LOGLEVEL_PRINT);
138
139     os << "idx(";
140
141     if (symbolic) {
142         os << "symbolic,name=" << name;
143     } else {
144         os << "non symbolic,value=" << value;
145     }
146
147     if (covariant) {
148         os << ",covariant";
149     } else {
150         os << ",contravariant";
151     }
152
153     os << ",serial=" << serial;
154     os << ",hash=" << hashvalue << ",flags=" << flags;
155     os << ")";
156 }
157
158 void idx::printtree(ostream & os, unsigned indent) const
159 {
160     debugmsg("idx printtree",LOGLEVEL_PRINT);
161
162     os << string(indent,' ') << "idx: ";
163
164     if (symbolic) {
165         os << "symbolic,name=" << name;
166     } else {
167         os << "non symbolic,value=" << value;
168     }
169
170     if (covariant) {
171         os << ",covariant";
172     } else {
173         os << ",contravariant";
174     }
175
176     os << ", serial=" << serial
177        << ", hash=" << hashvalue << " (0x" << hex << hashvalue << dec << ")"
178        << ", flags=" << flags << endl;
179 }
180
181 void idx::print(ostream & os, unsigned upper_precedence) const
182 {
183     debugmsg("idx print",LOGLEVEL_PRINT);
184
185     if (covariant) {
186         os << "_";
187     } else {
188         os << "~";
189     }
190     if (symbolic) {
191         os << name;
192     } else {
193         os << value;
194     }
195 }
196
197 bool idx::info(unsigned inf) const
198 {
199     if (inf==info_flags::idx) return true;
200     return basic::info(inf);
201 }
202
203 ex idx::subs(lst const & ls, lst const & lr) const
204 {
205     GINAC_ASSERT(ls.nops()==lr.nops());
206 #ifdef DO_GINAC_ASSERT
207     for (int i=0; i<ls.nops(); i++) {
208         GINAC_ASSERT(is_ex_exactly_of_type(ls.op(i),symbol)||
209                is_ex_of_type(ls.op(i),idx));
210     }
211 #endif // def DO_GINAC_ASSERT
212
213     for (int i=0; i<ls.nops(); i++) {
214         if (is_equal(*(ls.op(i)).bp)) {
215             return lr.op(i);
216         }
217     }
218     return *this;
219 }
220
221 // protected
222
223 int idx::compare_same_type(basic const & other) const
224 {
225     GINAC_ASSERT(is_of_type(other,idx));
226     idx const & o=static_cast<idx const &>
227                              (const_cast<basic &>(other));
228
229     if (covariant!=o.covariant) {
230         // different co/contravariant
231         return covariant ? -1 : 1;
232     }
233     if ((!symbolic) && (!o.symbolic)) {
234         // non-symbolic, of equal type: compare values
235         if (value==o.value) {
236             return 0;
237         }
238         return value<o.value ? -1 : 1;
239     }
240     if (symbolic && o.symbolic) {
241         // both symbolic: compare serials
242         if (serial==o.serial) {
243             return 0;
244         }
245         return serial<o.serial ? -1 : 1;
246     }
247     // one symbolic, one value: value is sorted first
248     return o.symbolic ? -1 : 1;
249 }
250
251 bool idx::is_equal_same_type(basic const & other) const
252 {
253     GINAC_ASSERT(is_of_type(other,idx));
254     idx const & o=static_cast<idx const &>
255                              (const_cast<basic &>(other));
256
257     if (covariant!=o.covariant) return false;
258     if (symbolic!=o.symbolic) return false;
259     if (symbolic && o.symbolic) return serial==o.serial;
260     return value==o.value;
261 }    
262
263 unsigned idx::calchash(void) const
264 {
265     hashvalue=golden_ratio_hash(golden_ratio_hash(tinfo_key ^ serial));
266     setflag(status_flags::hash_calculated);
267     return hashvalue;
268 }
269
270 //////////
271 // new virtual functions which can be overridden by derived classes
272 //////////
273
274 // public
275
276 bool idx::is_co_contra_pair(basic const & other) const
277 {
278     // like is_equal_same_type(), but tests for different covariant status
279     GINAC_ASSERT(is_of_type(other,idx));
280     idx const & o=static_cast<idx const &>
281                              (const_cast<basic &>(other));
282
283     if (covariant==o.covariant) return false;
284     if (symbolic!=o.symbolic) return false;
285     if (symbolic && o.symbolic) return serial==o.serial;
286     return value==o.value;
287 }    
288
289 bool idx::is_symbolic(void) const
290 {
291     return symbolic;
292 }
293
294 unsigned idx::get_value(void) const
295 {
296     return value;
297 }
298
299 bool idx::is_covariant(void) const
300 {
301     return covariant;
302 }
303
304 ex idx::toggle_covariant(void) const
305 {
306     idx * i_copy=static_cast<idx *>(duplicate());
307     i_copy->covariant = !i_copy->covariant;
308     i_copy->clearflag(status_flags::hash_calculated);
309     return i_copy->setflag(status_flags::dynallocated);
310 }
311
312 //////////
313 // non-virtual functions in this class
314 //////////
315
316 // none
317
318 //////////
319 // static member variables
320 //////////
321
322 // protected
323
324 unsigned idx::next_serial=0;
325
326 //////////
327 // global constants
328 //////////
329
330 const idx some_idx;
331 type_info const & typeid_idx=typeid(some_idx);
332
333 //////////
334 // other functions
335 //////////
336
337 int canonicalize_indices(exvector & iv, bool antisymmetric)
338 {
339     if (iv.size()<2) {
340         // nothing do to for 0 or 1 indices
341         return INT_MAX;
342     }
343
344     bool something_changed=false;
345     int sig=1;
346     // simple bubble sort algorithm should be sufficient for the small number of indices needed
347     exvector::const_iterator last_idx=iv.end();
348     exvector::const_iterator next_to_last_idx=iv.end()-1;
349     for (exvector::iterator it1=iv.begin(); it1!=next_to_last_idx; ++it1) {
350         for (exvector::iterator it2=it1+1; it2!=last_idx; ++it2) {
351             int cmpval=(*it1).compare(*it2);
352             if (cmpval==1) {
353                 iter_swap(it1,it2);
354                 something_changed=true;
355                 if (antisymmetric) sig=-sig;
356             } else if ((cmpval==0) && antisymmetric) {
357                 something_changed=true;
358                 sig=0;
359             }
360         }
361     }
362     return something_changed ? sig : INT_MAX;
363 }
364
365 exvector idx_intersect(exvector const & iv1, exvector const & iv2)
366 {
367     // build a vector of symbolic indices contained in iv1 and iv2 simultaneously
368     // assumes (but does not test) that each index occurs at most twice
369     exvector iv_intersect;
370     for (exvector::const_iterator cit1=iv1.begin(); cit1!=iv1.end(); ++cit1) {
371         GINAC_ASSERT(is_ex_of_type(*cit1,idx));
372         if (ex_to_idx(*cit1).is_symbolic()) {
373             for (exvector::const_iterator cit2=iv2.begin(); cit2!=iv2.end(); ++cit2) {
374                 GINAC_ASSERT(is_ex_of_type(*cit2,idx));
375                 if ((*cit1).is_equal(*cit2)) {
376                     iv_intersect.push_back(*cit1);
377                     break;
378                 }
379             }
380         }
381     }
382     return iv_intersect;
383 }
384
385 #define TEST_PERMUTATION(A,B,C,P) \
386     if ((iv3[B].is_equal(iv2[0]))&&(iv3[C].is_equal(iv2[1]))) { \
387         if (antisymmetric) *sig=P; \
388         return iv3[A]; \
389     }
390
391 ex permute_free_index_to_front(exvector const & iv3, exvector const & iv2,
392                                bool antisymmetric, int * sig)
393 {
394     // match (return value,iv2) to iv3 by permuting indices
395     // iv3 is always cyclic
396
397     GINAC_ASSERT(iv3.size()==3);
398     GINAC_ASSERT(iv2.size()==2);
399
400     *sig=1;
401     
402     TEST_PERMUTATION(0,1,2,  1);
403     TEST_PERMUTATION(0,2,1, -1);
404     TEST_PERMUTATION(1,0,2, -1);
405     TEST_PERMUTATION(1,2,0,  1);
406     TEST_PERMUTATION(2,0,1,  1);
407     TEST_PERMUTATION(2,1,0, -1);
408     throw(std::logic_error("permute_free_index_to_front(): no valid permutation found"));
409 }
410     
411 unsigned subs_index_in_exvector(exvector & v, ex const & is, ex const & ir)
412 {
413     exvector::iterator it;
414     unsigned replacements=0;
415     unsigned current_replacements;
416
417     GINAC_ASSERT(is_ex_of_type(is,idx));
418     GINAC_ASSERT(is_ex_of_type(ir,idx));
419    
420     for (it=v.begin(); it!=v.end(); ++it) {
421         current_replacements=count_index(*it,is);
422         if (current_replacements>0) {
423             (*it)=(*it).subs(is==ir);
424         }
425         replacements += current_replacements;
426     }
427     return replacements;
428 }
429
430 unsigned count_index(ex const & e, ex const & i)
431 {
432     exvector idxv=e.get_indices();
433     unsigned count=0;
434     for (exvector::const_iterator cit=idxv.begin(); cit!=idxv.end(); ++cit) {
435         if ((*cit).is_equal(i)) count++;
436     }
437     return count;
438 }
439
440 ex subs_indices(ex const & e, exvector const & idxv_subs,
441                 exvector const & idxv_repl)
442 {
443     GINAC_ASSERT(idxv_subs.size()==idxv_repl.size());
444     ex res=e;
445     for (unsigned i=0; i<idxv_subs.size(); ++i) {
446         res=res.subs(idxv_subs[i]==idxv_repl[i]);
447     }
448     return res;
449 }
450
451 #ifndef NO_GINAC_NAMESPACE
452 } // namespace GiNaC
453 #endif // ndef NO_GINAC_NAMESPACE