06eef8ee6b03b5058784291c62565ac66abfc78f
[ginac.git] / ginac / simp_lor.cpp
1 /** @file simp_lor.cpp
2  *
3  *  Implementation of GiNaC's simp_lor objects.
4  *  No real implementation yet, to be done.     */
5
6 /*
7  *  GiNaC Copyright (C) 1999-2000 Johannes Gutenberg University Mainz, Germany
8  *
9  *  This program is free software; you can redistribute it and/or modify
10  *  it under the terms of the GNU General Public License as published by
11  *  the Free Software Foundation; either version 2 of the License, or
12  *  (at your option) any later version.
13  *
14  *  This program is distributed in the hope that it will be useful,
15  *  but WITHOUT ANY WARRANTY; without even the implied warranty of
16  *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
17  *  GNU General Public License for more details.
18  *
19  *  You should have received a copy of the GNU General Public License
20  *  along with this program; if not, write to the Free Software
21  *  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
22  */
23
24 #include <string>
25 #include <list>
26 #include <algorithm>
27 #include <iostream>
28 #include <stdexcept>
29 #include <map>
30
31 #include "simp_lor.h"
32 #include "ex.h"
33 #include "mul.h"
34 #include "symbol.h"
35 #include "debugmsg.h"
36 #include "utils.h"
37
38 #ifndef NO_NAMESPACE_GINAC
39 namespace GiNaC {
40 #endif // ndef NO_NAMESPACE_GINAC
41
42 //////////
43 // default constructor, destructor, copy constructor assignment operator and helpers
44 //////////
45
46 // public
47
48 simp_lor::simp_lor() : type(invalid)
49 {
50         debugmsg("simp_lor default constructor",LOGLEVEL_CONSTRUCT);
51         tinfo_key=TINFO_simp_lor;
52 }
53
54 simp_lor::~simp_lor()
55 {
56         debugmsg("simp_lor destructor",LOGLEVEL_DESTRUCT);
57         destroy(false);
58 }
59
60 simp_lor::simp_lor(const simp_lor & other)
61 {
62         debugmsg("simp_lor copy constructor",LOGLEVEL_CONSTRUCT);
63         copy (other);
64 }
65
66 const simp_lor & simp_lor::operator=(const simp_lor & other)
67 {
68         debugmsg("simp_lor operator=",LOGLEVEL_ASSIGNMENT);
69         if (this != &other) {
70                 destroy(true);
71                 copy(other);
72         }
73         return *this;
74 }
75
76 // protected
77
78 void simp_lor::copy(const simp_lor & other)
79 {
80         indexed::copy(other);
81         type=other.type;
82         name=other.name;
83 }
84
85 void simp_lor::destroy(bool call_parent)
86 {
87         if (call_parent) {
88                 indexed::destroy(call_parent);
89         }
90 }
91
92 //////////
93 // other constructors
94 //////////
95
96 // protected
97
98 simp_lor::simp_lor(simp_lor_types const t) : type(t)
99 {
100         debugmsg("simp_lor constructor from simp_lor_types",LOGLEVEL_CONSTRUCT);
101         tinfo_key=TINFO_simp_lor;
102 }
103
104 simp_lor::simp_lor(simp_lor_types const t, const ex & i1, const ex & i2)
105   : indexed(i1,i2), type(t)
106 {
107         debugmsg("simp_lor constructor from simp_lor_types,ex,ex",LOGLEVEL_CONSTRUCT);
108         tinfo_key=TINFO_simp_lor;
109         GINAC_ASSERT(all_of_type_lorentzidx());
110 }
111
112 simp_lor::simp_lor(simp_lor_types const t, const std::string & n, const ex & i1)
113   : indexed(i1), type(t), name(n)
114 {
115         debugmsg("simp_lor constructor from simp_lor_types,string,ex",LOGLEVEL_CONSTRUCT);
116         tinfo_key=TINFO_simp_lor;
117         GINAC_ASSERT(all_of_type_lorentzidx());
118 }
119
120 simp_lor::simp_lor(simp_lor_types const t, const std::string & n, const exvector & iv)
121   : indexed(iv), type(t), name(n)
122 {
123         debugmsg("simp_lor constructor from simp_lor_types,string,exvector",LOGLEVEL_CONSTRUCT);
124         tinfo_key=TINFO_simp_lor;
125         GINAC_ASSERT(all_of_type_lorentzidx());
126 }
127
128 simp_lor::simp_lor(simp_lor_types const t, const std::string & n, exvector * ivp)
129   : indexed(ivp), type(t), name(n)
130 {
131         debugmsg("simp_lor constructor from simp_lor_types,string,exvector*",LOGLEVEL_CONSTRUCT);
132         tinfo_key=TINFO_simp_lor;
133         GINAC_ASSERT(all_of_type_lorentzidx());
134 }
135
136 //////////
137 // functions overriding virtual functions from bases classes
138 //////////
139
140 // public
141
142 basic * simp_lor::duplicate() const
143 {
144         debugmsg("simp_lor duplicate",LOGLEVEL_DUPLICATE);
145         return new simp_lor(*this);
146 }
147
148 void simp_lor::printraw(std::ostream & os) const
149 {
150         debugmsg("simp_lor printraw",LOGLEVEL_PRINT);
151         os << "simp_lor(type=" << (unsigned)type
152            << ",name=" << name << ",indices=";
153         printrawindices(os);
154         os << ",hash=" << hashvalue << ",flags=" << flags << ")";
155 }
156
157 void simp_lor::printtree(std::ostream & os, unsigned indent) const
158 {
159         debugmsg("simp_lor printtree",LOGLEVEL_PRINT);
160         os << std::string(indent,' ') << "simp_lor object: "
161            << "type=" << (unsigned)type
162            << ", name=" << name << ", ";
163         os << seq.size() << " indices" << std::endl;
164         printtreeindices(os,indent);
165         os << std::string(indent,' ') << "hash=" << hashvalue
166            << " (0x" << std::hex << hashvalue << std::dec << ")"
167            << ", flags=" << flags << std::endl;
168 }
169
170 void simp_lor::print(std::ostream & os, unsigned upper_precedence) const
171 {
172         debugmsg("simp_lor print",LOGLEVEL_PRINT);
173         switch (type) {
174         case simp_lor_g:
175                 os << "g";
176                 break;
177         case simp_lor_vec:
178                 os << name;
179                 break;
180         case invalid:
181         default:
182                 os << "INVALID_SIMP_LOR_OBJECT";
183                 break;
184         }
185         printindices(os);
186 }
187
188 void simp_lor::printcsrc(std::ostream & os, unsigned type, unsigned upper_precedence) const
189 {
190         debugmsg("simp_lor print csrc",LOGLEVEL_PRINT);
191         print(os,upper_precedence);
192 }
193
194 bool simp_lor::info(unsigned inf) const
195 {
196         return indexed::info(inf);
197 }
198
199 ex simp_lor::eval(int level) const
200 {
201         if (type==simp_lor_g) {
202                 // canonicalize indices
203                 exvector iv=seq;
204                 int sig=canonicalize_indices(iv,false); // symmetric
205                 if (sig!=INT_MAX) {
206                         // something has changed while sorting indices, more evaluations later
207                         if (sig==0) return _ex0();
208                         return ex(sig)*simp_lor(type,name,iv);
209                 }
210                 const lorentzidx & idx1=ex_to_lorentzidx(seq[0]);
211                 const lorentzidx & idx2=ex_to_lorentzidx(seq[1]);
212                 if ((!idx1.is_symbolic())&&(!idx2.is_symbolic())) {
213                         // both indices are numeric
214                         if ((idx1.get_value()==idx2.get_value())) {
215                                 // both on diagonal
216                                 if (idx1.get_value()==0) {
217                                         // (0,0)
218                                         return _ex1();
219                                 } else {
220                                         if (idx1.is_covariant()!=idx2.is_covariant()) {
221                                                 // (_i,~i) or (~i,_i), i=1..3
222                                                 return _ex1();
223                                         } else {
224                                                 // (_i,_i) or (~i,~i), i=1..3
225                                                 return _ex_1();
226                                         }
227                                 }
228                         } else {
229                                 // at least one off-diagonal
230                                 return _ex0();
231                         }
232                 } else if (idx1.is_symbolic() &&
233                            idx1.is_co_contra_pair(idx2)) {
234                         return Dim()-idx1.get_dim_parallel_space();
235                 }
236         }
237
238         return this->hold();
239 }
240         
241 // protected
242
243 int simp_lor::compare_same_type(const basic & other) const
244 {
245         GINAC_ASSERT(other.tinfo() == TINFO_simp_lor);
246         const simp_lor *o = static_cast<const simp_lor *>(&other);
247         if (type==o->type) {
248                 if (name==o->name) {
249                         return indexed::compare_same_type(other);
250                 }
251                 return name.compare(o->name);
252         }
253         return type < o->type ? -1 : 1;
254 }
255
256 bool simp_lor::is_equal_same_type(const basic & other) const
257 {
258         GINAC_ASSERT(other.tinfo() == TINFO_simp_lor);
259         const simp_lor *o = static_cast<const simp_lor *>(&other);
260         if (type!=o->type) return false;
261         if (name!=o->name) return false;
262         return indexed::is_equal_same_type(other);
263 }
264
265 unsigned simp_lor::return_type(void) const
266 {
267         return return_types::commutative;
268 }
269    
270 unsigned simp_lor::return_type_tinfo(void) const
271 {
272         return tinfo_key;
273 }
274
275 ex simp_lor::thisexprseq(const exvector & v) const
276 {
277         return simp_lor(type,name,v);
278 }
279
280 ex simp_lor::thisexprseq(exvector * vp) const
281 {
282         return simp_lor(type,name,vp);
283 }
284
285 //////////
286 // virtual functions which can be overridden by derived classes
287 //////////
288
289 // none
290
291 //////////
292 // non-virtual functions in this class
293 //////////
294
295 // protected
296
297 bool simp_lor::all_of_type_lorentzidx(void) const
298 {
299         // used only inside of ASSERTs
300         for (exvector::const_iterator cit=seq.begin(); cit!=seq.end(); ++cit) {
301                 if (!is_ex_of_type(*cit,lorentzidx)) return false;
302         }
303         return true;
304 }
305
306 //////////
307 // static member variables
308 //////////
309
310 // none
311
312 //////////
313 // global constants
314 //////////
315
316 const simp_lor some_simp_lor;
317 const type_info & typeid_simp_lor = typeid(some_simp_lor);
318
319 //////////
320 // friend functions
321 //////////
322
323 simp_lor lor_g(const ex & mu, const ex & nu)
324 {
325         return simp_lor(simp_lor::simp_lor_g,mu,nu);
326 }
327
328 simp_lor lor_vec(const std::string & n, const ex & mu)
329 {
330         return simp_lor(simp_lor::simp_lor_vec,n,mu);
331 }
332
333 ex simplify_simp_lor_mul(const ex & m, const scalar_products & sp)
334 {
335         GINAC_ASSERT(is_ex_exactly_of_type(m,mul));
336         exvector v_contracted;
337
338         // collect factors in an exvector, store squares twice
339         unsigned n=m.nops();
340         v_contracted.reserve(2*n);
341         for (unsigned i=0; i<n; ++i) {
342                 ex f=m.op(i);
343                 if (is_ex_exactly_of_type(f,power)&&f.op(1).is_equal(_ex2())) {
344                         v_contracted.push_back(f.op(0));
345                         v_contracted.push_back(f.op(0));
346                 } else {
347                         v_contracted.push_back(f);
348                 }
349         }
350
351         unsigned replacements;
352         bool something_changed=false;
353
354         exvector::iterator it=v_contracted.begin();
355         while (it!=v_contracted.end()) {
356                 // process only lor_g objects
357                 if (is_ex_exactly_of_type(*it,simp_lor) &&
358                         (ex_to_simp_lor(*it).type==simp_lor::simp_lor_g)) {
359                         const simp_lor & g=ex_to_simp_lor(*it);
360                         GINAC_ASSERT(g.seq.size()==2);
361                         const idx & first_idx=ex_to_lorentzidx(g.seq[0]);
362                         const idx & second_idx=ex_to_lorentzidx(g.seq[1]);
363                         // g_{mu,mu} should have been contracted in simp_lor::eval()
364                         GINAC_ASSERT(!first_idx.is_equal(second_idx));
365                         ex saved_g=*it; // save to restore it later
366
367                         // try to contract first index
368                         replacements=0;
369                         if (first_idx.is_symbolic()) {
370                                 replacements = subs_index_in_exvector(v_contracted, first_idx.toggle_covariant(),second_idx);
371                                 if (replacements==0) {
372                                         // not contracted, restore g object
373                                         *it=saved_g;
374                                 } else {
375                                         // a contracted index should occur exactly once
376                                         GINAC_ASSERT(replacements==1);
377                                         *it=_ex1();
378                                         something_changed=true;
379                                 }
380                         }
381
382                         // try second index only if first was not contracted
383                         if ((replacements==0)&&(second_idx.is_symbolic())) {
384                                 // first index not contracted, *it is again the original g object
385                                 replacements = subs_index_in_exvector(v_contracted, second_idx.toggle_covariant(),first_idx);
386                                 if (replacements==0) {
387                                         // not contracted except in itself, restore g object
388                                         *it=saved_g;
389                                 } else {
390                                         // a contracted index should occur exactly once
391                                         GINAC_ASSERT(replacements==1);
392                                         *it=_ex1();
393                                         something_changed=true;
394                                 }
395                         }
396                 }
397                 ++it;
398         }
399
400         // process only lor_vec objects
401         bool jump_to_next=false;
402         exvector::iterator it1=v_contracted.begin();
403         while (it1!=v_contracted.end()-1) {
404                 if (is_ex_exactly_of_type(*it1,simp_lor) && 
405                         (ex_to_simp_lor(*it1).type==simp_lor::simp_lor_vec)) {
406                         exvector::iterator it2=it1+1;
407                         while ((it2!=v_contracted.end())&&!jump_to_next) {
408                                 if (is_ex_exactly_of_type(*it2,simp_lor) && 
409                                         (ex_to_simp_lor(*it2).type==simp_lor::simp_lor_vec)) {
410                                         const simp_lor & vec1=ex_to_simp_lor(*it1);
411                                         const simp_lor & vec2=ex_to_simp_lor(*it2);
412                                         GINAC_ASSERT(vec1.seq.size()==1);
413                                         GINAC_ASSERT(vec2.seq.size()==1);
414                                         const lorentzidx & idx1=ex_to_lorentzidx(vec1.seq[0]);
415                                         const lorentzidx & idx2=ex_to_lorentzidx(vec2.seq[0]);
416                                         if (idx1.is_symbolic() &&
417                                                 idx1.is_co_contra_pair(idx2) &&
418                                                 sp.is_defined(vec1,vec2)) {
419                                                 *it1=sp.evaluate(vec1,vec2);
420                                                 *it2=_ex1();
421                                                 something_changed=true;
422                                                 jump_to_next=true;
423                                         }
424                                 }
425                                 ++it2;
426                         }
427                         jump_to_next=false;
428                 }
429                 ++it1;
430         }
431         if (something_changed) {
432                 return mul(v_contracted);
433         }
434         return m;
435 }
436
437 ex simplify_simp_lor(const ex & e, const scalar_products & sp)
438 {
439         // all simplification is done on expanded objects
440         ex e_expanded=e.expand();
441
442         // simplification of sum=sum of simplifications
443         if (is_ex_exactly_of_type(e_expanded,add)) {
444                 ex sum=_ex0();
445                 for (unsigned i=0; i<e_expanded.nops(); ++i)
446                         sum += simplify_simp_lor(e_expanded.op(i),sp);
447
448                 return sum;
449         }
450
451         // simplification of commutative product=commutative product of simplifications
452         if (is_ex_exactly_of_type(e_expanded,mul)) {
453                 return simplify_simp_lor_mul(e,sp);
454         }
455
456         // cannot do anything
457         return e_expanded;
458 }
459
460 //ex Dim(void)   // FIXME: what's going on here?
461 //{
462 //    static symbol * d=new symbol("dim");
463 //    return *d;
464 //}
465
466 //////////
467 // helper classes
468 //////////
469
470 void scalar_products::reg(const simp_lor & v1, const simp_lor & v2,
471                           const ex & sp)
472 {
473         if (v1.compare_same_type(v2)>0) {
474                 reg(v2,v1,sp);
475                 return;
476         }
477         spm[make_key(v1,v2)]=sp;
478 }
479
480 bool scalar_products::is_defined(const simp_lor & v1, const simp_lor & v2) const
481 {
482         if (v1.compare_same_type(v2)>0) {
483                 return is_defined(v2,v1);
484         }
485         return spm.find(make_key(v1,v2))!=spm.end();
486 }
487
488 ex scalar_products::evaluate(const simp_lor & v1, const simp_lor & v2) const
489 {
490         if (v1.compare_same_type(v2)>0) {
491                 return evaluate(v2,v1);
492         }
493         return spm.find(make_key(v1,v2))->second;
494 }
495
496 void scalar_products::debugprint(void) const
497 {
498         std::cerr << "map size=" << spm.size() << std::endl;
499         for (spmap::const_iterator cit=spm.begin(); cit!=spm.end(); ++cit) {
500                 const spmapkey & k=(*cit).first;
501                 std::cerr << "item key=((" << k.first.first
502                           << "," << k.first.second << "),";
503                 k.second.printraw(cerr);
504                 cerr << ") value=" << (*cit).second << std::endl;
505         }
506 }
507
508 spmapkey scalar_products::make_key(const simp_lor & v1, const simp_lor & v2)
509 {
510         GINAC_ASSERT(v1.type==simp_lor::simp_lor_vec);
511         GINAC_ASSERT(v2.type==simp_lor::simp_lor_vec);
512         lorentzidx anon=ex_to_lorentzidx(v1.seq[0]).create_anonymous_representative();
513         GINAC_ASSERT(anon.is_equal_same_type(ex_to_lorentzidx(v2.seq[0]).create_anonymous_representative()));
514         return spmapkey(strstrpair(v1.name,v2.name),anon);
515 }
516
517 #ifndef NO_NAMESPACE_GINAC
518 } // namespace GiNaC
519 #endif // ndef NO_NAMESPACE_GINAC